Estimating the lateral transfer of organic carbon through the European river network using a land surface model

Lateral carbon transport from soils to the ocean through rivers has been acknowledged as a key component of the global carbon cycle, but it is still neglected in most global land surface models (LSMs). Fluvial transport of dissolved organic carbon (DOC) and CO₂ has been implemented in the ORCHIDEE LSM, while erosion-induced delivery of sediment and particulate organic carbon (POC) from land to river was implemented in another version of the model. Based on these two developments, we take the final step towards the full representation of biospheric carbon transport through the land–river continuum. The newly developed model, called ORCHIDEE-C_lateral, simulates the complete lateral transport of water, sediment, POC, DOC, and CO₂ from land to sea through the river network, the deposition of sediment and POC in the river channel and floodplains, and the decomposition of POC and DOC in transit. We parameterized and evaluated ORCHIDEE-C_lateral using observation data in Europe. The model explains 94 %, 75 %, and 83 % of the spatial variations of observed riverine water discharges, bankfull water flows, and riverine sediment discharges in Europe, respectively. The simulated long-term average total organic carbon concentrations and DOC concentrations in river flows are comparable to the observations in major European rivers, although our model generally overestimates the seasonal variation of riverine organic carbon concentrations. Application of ORCHIDEE-C_lateral for Europe reveals that the lateral carbon transfer affects land carbon dynamics in multiple ways, and omission of this process in LSMs may lead to an overestimation of 4.5 % in the simulated annual net terrestrial carbon uptake over Europe. Overall, this study presents a useful tool for simulating large-scale lateral carbon transfer and for predicting the feedbacks between lateral carbon transfer and future climate and land use changes.

1 Introduction

Lateral transfer of organic carbon along the land–river–ocean continuums, involving both spatial redistribution of terrestrial organic carbon and the vertical land–atmosphere carbon exchange, has been acknowledged as a key component of the global carbon cycle (Ciais et al., 2013, 2021; Drake et al., 2018; Regnier et al., 2013, 2022). Erosion of soils and the associated organic carbon, but also leaching of soil dissolved organic carbon (DOC), represent a non-negligible leak in the terrestrial carbon budget and a substantial source of allochthonous organic carbon to inland waters and oceans (Battin et al., 2009; Cole et al., 2007; Raymond et al., 2013; Regnier et al., 2013). As a result of soil aggregate breakdown and desorption, the accelerated mineralization of these eroded and leached soil carbon loads leads to considerable CO₂ emission to the atmosphere (Chappell et al., 2016; Lal, 2003; Van Hemelryck et al., 2011). Meanwhile, the organic carbon that is redeposited and buried in floodplains and lakes might be preserved for a long time, thus creating a CO₂ sink (Stallard, 1998; Van Oost et al., 2007; Wang et al., 2010; Hoffmann, 2022). In addition, lateral redistribution of soil material can alter land–atmosphere CO₂ fluxes indirectly by affecting soil nutrient availability, terrestrial vegetation productivity, and physiochemical properties of inland and coastal waters (Beusen et al., 2005; Vigiak et al., 2017).

Although the important role of lateral carbon transfer in the global carbon cycle has been widely recognized, to date, the estimates of land carbon loss to inland waters, the fate of the terrestrial organic carbon within inland waters, and the net effect of lateral carbon transfer on land–atmosphere CO₂ fluxes remain largely uncertain (Berhe et al., 2007; Doetterl et al., 2016; Lal, 2003; Stallard, 1998; Z. Wang et al., 2014; Zhang et al., 2014). Existing estimates of global carbon loss from soils to inland waters vary from 1.1 to 5.1 Pg (=10¹⁵ g) C per year (yr⁻¹) (Cole et al., 2007; Drake et al., 2018), and the estimated net impact of global lateral carbon redistribution on the land–atmosphere carbon budget ranges from an uptake of atmospheric CO₂ by 1 Pg C yr⁻¹ to a land CO₂ emission of 1 Pg C yr⁻¹ (Lal, 2003; Stallard, 1998; Van Oost et al., 2007; Wang et al., 2017; Regnier et al., 2022). A reliable model which is able to explicitly simulate the lateral carbon flux along the land–river continuum and also the interactions between these lateral fluxes and the comprehensive terrestrial carbon cycle would thus be necessary for projecting changes in the global carbon cycle more accurately.

Global land surface models (LSMs) are important tools to simulate the feedbacks between terrestrial carbon cycle, increasing atmospheric CO₂, and climate and land use change. However, the lateral carbon transfer, especially for particulate organic carbon (POC), is still missing or incompletely represented in existing LSMs (Lauerwald et al., 2017, 2020; Lugato et al., 2016; Naipal et al., 2020; Nakhavali et al., 2021; Tian et al., 2015). It has been hypothesized that the exclusion of lateral carbon transfer in LSMs implies a significant bias in the simulated global land carbon budget (Ciais et al., 2013, 2021; Janssens et al., 2003). For instance, the study of Nakhavali et al. (2021) suggested that about 15 % of the global terrestrial net ecosystem production is exported to inland waters as leached DOC. Lauerwald et al. (2020) showed that the omission of lateral DOC transfer in LSMs might lead to significant underestimation (8.6 %) of the net uptake of atmospheric carbon in the Amazon basin, while terrestrial carbon storage changes in response to the increasing atmospheric CO₂ concentrations were overestimated.

Over the past decade, a number of LSMs have been developed which represent leaching of DOC from soils (Nakhavali et al., 2021; Kicklighter et al., 2013) or the full transport of DOC through the land–river continuum (Lauerwald et al., 2017; Tian et al., 2015). However, the erosion-induced transport of soil POC, which has also been reported to be able to strongly affect the carbon balance of terrestrial ecosystems (Lal, 2003; Van Oost et al., 2007; Tian et al., 2015), is still not represented or is poorly represented in LSMs. The explicit simulation of the complete transport process of POC at large spatial scales is still a major challenge due to the complexity of the processes involved, including erosion-induced sediment and POC delivery to rivers, deposition of sediment and POC in river channels and floodplains, re-detachment of the previously deposited sediments and POC, decomposition and transformation of POC in riverine and flooding waters, and the changes in the soil profile caused by erosion and deposition (Doetterl et al., 2016; Naipal et al., 2020; Zhang et al., 2020).

Several recent model developments have led to the implementation of the lateral transfer of POC in large-scale LSMs. Despite this, there are still some inevitable limitations in these implementations. The Dynamic Land Ecosystem Model (DLEM v2.0, Tian et al., 2015) is able to simulate the erosion-induced POC loss from soil to river and the transport and decomposition of POC in river networks. However, it does not represent the POC deposition in floodplains or the impacts of soil erosion and floodplain deposition on the vertical profiles of soil organic carbon (SOC). The Carbon Erosion DYNAMics model (CE-DYNAM, Naipal et al., 2020) simulates erosion of SOC and its redeposition on the toe slope or floodplains, transport of POC along river channels, and the impact on SOC dynamics at the eroding and deposition sites. However, running at an annual timescale, it mostly addresses the centennial timescale and does not represent deposition and decomposition of POC in river channels. Moreover, CE-DYNAM was only applied over the Rhine catchment and has not been fully coupled into a land surface model, therefore excluding the feedbacks of soil erosion on the fully coupled land and aquatic carbon cycles. There are of course more dedicated hydrology and soil erosion models that explicitly simulate the complete transport, deposition, and decomposition processes of POC in small river basins (e.g., Jetten et al., 2003; Nearing et al., 1989; Neitsch et al., 2011). However, it is difficult to apply these models at large spatial scales (e.g., continental or global scale) due to the limited availability of forcing data (e.g., geometric attributes of river channel), suitable model parameterization, and computational capacity. Moreover, these models have limited capability to represent the full terrestrial C cycle in response to climate change, increasing atmospheric CO₂, and land use change. Therefore, basin-scale models are not an option to assess the impact of soil erosion on the large-scale terrestrial C budget in response to global changes.

Here we describe the development, application, and evaluation of a new branch of the ORCHIDEE LSM (Krinner et al., 2005), hereafter ORCHIDEE-C_lateral, that can be used to simulate the complete lateral transfer processes of water, sediment, POC, and DOC along the land–river–ocean continuum at large spatial scale (e.g., continental and global scale). In previous studies, the leaching and fluvial transfer of DOC and the erosion-induced delivery of sediment and POC from upland soil to river network have been implemented in two different branches of the ORCHIDEE LSM (i.e., ORCHILEAK, Lauerwald et al., 2017, and ORCHIDEE-MUSLE, Zhang et al., 2020). For this new branch, we first merged these two branches and subsequently implemented the fluvial transfer of sediment and POC in the coupled model. ORCHIDEE-C_lateral is calibrated and evaluated using observation data on runoff, bankfull flow, and riverine loads, as well as concentrations of sediment, POC, and DOC across Europe. By applying the calibrated model at the European scale, we estimate the magnitude and spatial distribution of the lateral carbon transfer in European catchments during the period 1901–2014, as well as the potential impacts of lateral carbon transfer on the land carbon balance. Comparing simulation results to those of an alternative simulation run with lateral displacement of C deactivated, we finally quantify the biases in simulated land C budgets that arise ignoring the lateral transfers of C along the land–river continuum.

2 Model development and evaluation

2.1 ORCHIDEE land surface model

The ORCHIDEE LSM comprehensively simulates the cycling of energy, water, and carbon in terrestrial ecosystems (Krinner et al., 2005). The hydrological processes (e.g., rainfall interception, evapotranspiration, and soil water dynamics) and plant photosynthesis in ORCHIDEE are simulated at a time step of 30 min. The carbon cycle processes (e.g., maintenance and growth respiration, carbon allocation, litter decomposition, SOC dynamics, plant phenology and mortality) are simulated at daily time step. In its default configuration, ORCHIDEE represents 13 land cover types, with one for bare soil and 12 for lands covered by vegetation (eight types of forests, two types of grasslands, two types of croplands). Given appropriate land cover maps and parametrization, the number of plant functional types (PFTs) to be represented can, however, be adapted (Zhang et al., 2020).

Our previous implementations of lateral DOC transfer (Lauerwald et al., 2017) and POC delivery from upland to river network (Zhang et al., 2020) were both based on the ORCHIDEE branch ORCHIDEE-SOM (Camino-Serrano et al., 2018, Fig. S1), which provides a depth-dependent description of the water and carbon dynamics in the soil column. Specifically, the vertical soil profile in ORCHIDEE-SOM is described by an 11-layer discretization of a 2 m soil column (Camino-Serrano et al., 2018). Water flows between adjacent soil layers are simulated using the Fokker–Planck equation that resolves water diffusion in non-saturated conditions (Campoy et al., 2013; Guimberteau et al., 2018). Free gravitational drainage occurs in the lowest soil layer when actual soil water content is higher than the residual water content (Campoy et al., 2013). Following the CENTURY model (Parton et al., 1988), ORCHIDEE-SOM represents two litter pools (metabolic and structural) and three SOC pools (active, slow, and passive) that differ in their respective turnover times. The decomposition of each carbon pool is calculated by first-order kinetics based on the corresponding turnover time, soil moisture and temperature as controlling factors, and the priming effects of fresh organic matter (Guenet et al., 2016, 2018). Soil DOC is represented by a labile and a refractory DOC pool, with a high and low turnover rate, respectively. Each DOC pool may be in the soil solution or adsorbed on the mineral matrix. The products of litter and SOC decomposition enter the free DOC pool, which in turn is decomposed following first-order kinetics (Kalbitz et al., 2003) and returns back to SOC. Adsorption and desorption of DOC follow an equilibrium distribution coefficient calculated from soil clay and pH. Free DOC can be transported with the water flux simulated by the soil hydrological module of ORCHIDEE. However, DOC adsorbed to soil minerals can neither be decomposed nor transported (Camino-Serrano et al., 2018). All the described processes occur within each soil layer. At each time step, “the flux of DOC leaving the soil is calculated by multiplying DOC concentrations in soil solution with the runoff (surface layer) and drainage (bottom layer) flux simulated by the hydrological module” (Camino-Serrano et al., 2018, p. 939). More detailed information about the simulation of soil hydrological and biogeochemical processes in ORCHIDEE-SOM can be found in Guenet et al. (2016) and Camino-Serrano et al. (2018).

2.1.1 Lateral transfer of DOC and CO₂

Lateral transfer of DOC and dissolved CO₂ from land to ocean through river network has been implemented in ORCHILEAK (Lauerwald et al., 2017), an ORCHIDEE branch developed from ORCHIDEE-SOM (Fig. S1). The method used in ORCHILEAK to simulate the adsorption, desorption, production, consumption, and transport of DOC within the soil column, as well as DOC export from the soil column with surface runoff and drainage, is similar to that used in ORCHIDEE-SOM. Besides the decomposition of SOC and litter, ORCHILEAK also represents the contribution of wet and dry deposition to soil DOC via throughfall. The direct DOC input from rainfall to aquatic DOC pools is simulated based on the DOC concentration in rainfall and the area fraction of stream and flooding waters in each basin (Table 1). Note that the maximum area fractions of river surface and floodplain in each basin (i.e., each $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cell in this study) are derived from high-resolution topographic data (Table 1). As it is difficult to explicitly represent all real river channels in a global land surface model (due to the limit of computing efficiency of current computers), we assume that there is one virtual river channel in each $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ pixel. The surface area of this virtual river is the sum of all real rivers, and the flow direction of this virtual is assumed to be same to the largest real river (Lauerwald et al., 2015).

Table 1 List of forcing data needed to run ORCHIDEE-C_lateral and the data used to evaluate the simulation results. S_res and T_res are the spatial and temporal resolution of the forcing data, respectively.

^a GFPLAIN250m only covers the regions south of 60^∘ N. We produced a map of floodplain distribution in regions north of the 60^∘ N using the same method for producing GFPLAIN250m (Nardi et al., 2019) based on the ASTER GDEM v3 database (Abrams et al., 2020). ^b The DEM data from HydroSHEDS and GDEM v3 are used to extract the topographic properties (e.g., location, area, and average slope) of headwater basins in regions south and north of 60^∘ N, respectively.^c The Global Runoff Data Centre, 56068 Koblenz, Germany. ^d https://www.eea.europa.eu/data-and-maps/data/sediment-discharges (last access: 12 May 2020). ^e Publications including Van Dijk and Kwaad (1998), Vollmer and Goelz (2006), and reports of the DanubeSediment project (Sediment Management Measures for the Danube, http://www.interreg-danube.eu/approved-projects/danubesediment, last access: 24 June 2020). ^f https://zenodo.org/record/2536040#.YC-QGo9KiUm (last access: 14 October 2020).

Simulation of the lateral transfer of DOC and CO₂ in river networks, i.e., the transfer of DOC and CO₂ from one basin to another based on the streamflow directions obtained from a forcing file (0.5^∘, Table 1), follows the routing scheme of water (Guimberteau et al., 2012). For each basin with a floodplain (defined by forcing data), bankfull flow occurs when stream volume in the river channel exceeds a threshold prescribed by the forcing file (Table 1). DOC and CO₂ in flooding waters can enter into soil DOC and CO₂ pools along with the flooding water infiltrated into soil. The infiltration rate of flooding water depends on soil properties and soil water content, but it does not depend on vegetation cover. On the contrary, DOC and CO₂ originating from the decomposition of submerged litter and SOC in the floodplains is added to the overlying flooding waters. Note that the turnover times of litter and SOC under flooding waters are assumed to be 3 times the litter and SOC turnover times in upland soil (Reddy and Patrick, 1975; Neckles and Neill, 1994; Lauerwald et al., 2017). After removing the infiltrated and evaporated water, the amount of remaining flooding water, as well as the DOC and dissolved CO₂ returning to the river channel at the end of each day, is calculated based on a time constant of flooding water (4.0 d, d'Orgeval et al., 2008) modified by a basin-specific topographic index (f_topo, unitless) (Lauerwald et al., 2017).

Decomposition of DOC in stream and flooding waters is calculated at daily time step based on the prescribed turnover times of labile (2 d) and refractory (80 d) DOC in waters (when temperature is 28 ^∘C) and a temperature factor obtained from Hanson et al. (2011). CO₂ evasion in inland waters is simulated using a much finer integration time step of 6 min. The CO₂ partial pressures (pCO₂) in the water column are first calculated based on the temperature-dependent solubility of CO₂ and the concentration of dissolved CO₂ (Telmer and Veizer, 1999). Then the CO₂ evasion is calculated based on the gas exchange velocity, the water–air gradient in pCO₂, and the surface water area available for gas exchange (Lauerwald et al., 2017). The effect of wind speed on CO₂ evasion is not represented in the current version of ORCHILEAK. In addition, swamp and wetland are represented in the routing scheme of ORCHILEAK. More detailed descriptions can be found in Lauerwald et al. (2017).

2.1.2 Sediment and particulate organic carbon delivery from upland soil to river network

To give an accurate simulation of sediment delivery from uplands to the river network and maintain computational efficiency, an upscaling scheme which integrates information from high-resolution (3") topographic and soil erodibility data into a LSM forcing file at 0.5^∘ spatial resolution has been introduced (see details in Zhang et al., 2020, Fig. 1). With this upscaling scheme, the erosion-induced sediment and POC delivery from upland soils to the river network, as well as the changes in SOC profiles due to soil erosion, had already been implemented in ORCHIDEE-MUSLE (Zhang et al., 2020). The sediment delivery from small headwater basins (which are basins without perennial streams and are extracted from high-resolution, e.g., 3", digital elevation model – DEM – data; Fig. 1a, d) to the river network (i.e., gross upland soil erosion – sediment deposition within headwater basins) is simulated using the Modified Universal Soil Loss Equation model (MUSLE, Williams, 1975). As introduced in Zhang et al. (2020), “the daily sediment delivery rate from each headwater basin (S_{i_ref}, Mg d⁻¹ basin⁻¹) is first calculated for a given set of reference runoff and vegetation cover conditions (Fig. 1e):

$$\begin{array}{}\text{(1)}& S_{i\mathrm{\_}\mathrm{ref}}=a{\left(Q_{i\mathrm{\_}\mathrm{ref}}{q}_{i\mathrm{\_}\mathrm{ref}}\right)}^b{K}_i{\mathrm{LS}}_i{C}_{\mathrm{ref}}{P}_{\mathrm{ref}},\end{array}$$

where Q_{i_ref} is the total water discharge (m³ d⁻¹) at the outlet of headwater basin i for the daily reference runoff condition (R_ref) of 10 mm d⁻¹ (see Table S1 for the definitions of all abbreviations used in this study). In Eq. (1), q_{i_ref} is the daily peak flow rate (m³ s⁻¹) at the headwater basin outlet under the assumed reference runoff condition. Similar to the SWAT model (Soil and Water Assessment Tool, Neitsch et al., 2011), q_{i_ref} was calculated from the reference maximum 30 min runoff (1 mm 30 min⁻¹) depth and drainage area (DA_i, m²) according to the following equation:

$$\begin{array}{}\text{(2)}& q_{i\mathrm{\_}\mathrm{ref}}={\displaystyle \frac{R_{\mathrm{30}\mathrm{\_}\mathrm{ref}}}{\mathrm{30}\times \mathrm{60}}}\left({\mathrm{DA}}_i^{\left(d{{\mathrm{DA}}_i}^c\right)}\right)\mathrm{1000},\end{array}$$

where R_{30_ref} (=1 mm 30 min⁻¹) is the assumed daily maximum 30 min runoff”. The coefficients a and b in Eq. (1) and c and d in Eq. (2) need to be calibrated (see Sect. 2.3 and Table 2). In Eq. (1), the term LS_i is the combined dimensionless slope length and steepness factor calculated based on the DA_i and the average slope steepness (extracted from a DEM) of headwater basin i (Moore and Wilson, 1992). C_ref (0–1, dimensionless) in Eq. (1) represents the cover management factor, which depends on vegetation cover and storage of plant debris (see below). The value of C_ref is set to 0.1 for the reference state. The soil erodibility factor K_i (Mg MJ⁻¹ mm⁻¹) is calculated using the method of the EPIC model (Sharpley and Williams, 1990) based on SOC and soil texture data obtained from the GSDE database (Table 1). The term P_ref (0–1, dimensionless) in Eq. (1) is a factor representing erosion control practices. It was set to 1, as we did not consider the impacts of soil conservation practices in reducing soil erosion rate. Note that it does not matter which value is chosen for R_ref, R_{30_ref}, and C_ref as long as they are used consistently throughout a study.

Figure 1 Upscaling scheme used in ORCHIDEE-MUSLE (Zhang et al., 2020) and ORCHIDEE-C_lateral for calculating the sediment delivery rate from headwater basins to river networks. MUSLE is the Modified Universal Soil Loss Equation; DEM is the digital elevation model (m); K is the soil erodibility factor (Mg MJ⁻¹ mm⁻¹); R_ref is the assumed reference daily runoff depth (10 mm d⁻¹); R_{30_ref} is the assumed reference maximum 30 min runoff depth (1 mm 30 min⁻¹); C_ref (0.1, dimensionless) is the assumed reference cover management factor; R_iday, R_{30_iday}, and C_iday are the simulated daily total runoff depth, daily maximum 30 min runoff depth, and daily cover management factor, respectively. This figure is adapted from Fig. 1 in Zhang et al. (2020).

Table 2 Values of the key parameters used in ORCHIDEE-C_lateral to simulate the lateral transfer of sediment and carbon.

^a For clay, silt, and sand sediment, respectively. ^b For active, slow, and passive POC, respectively.

For the use of these reference sediment delivery estimates in ORCHIDEE-C_lateral, the values were first calculated for each headwater basin derived from high-resolution geodata (Fig. 1e), then aggregated to 0.5^∘ grid cells (Fig. 1f) – the scale used in our simulations and required to maintain computational efficiency (also limited by the availability of climate and land cover forcing data).

This aggregated dataset is then used to force the simulation of the actual daily sediment delivery (S_j, g d⁻¹ grid⁻¹) in ORCHIDEE-C_lateral, simply based on the estimated reference sediment delivery rates of Eq. (1) and on the ratios between actual runoff and land cover conditions as well as the assumed reference conditions used to create that forcing file (Eq. 4, Fig. 1g):

$$\begin{array}{}\text{(3)}& {\displaystyle }{S}_{\mathrm{ref}}={\sum }_{i=\mathrm{1}}^n\left(S_{i\mathrm{\_}\mathrm{ref}}\right)\times {\mathrm{10}}^{\mathrm{6}},\text{(4)}& {\displaystyle }{S}_j=S_{\mathrm{ref}}{\left({\displaystyle \frac{R_j{R}_{\mathrm{30}\mathit{\_}j}}{R_{\mathrm{ref}}{R}_{\mathrm{30}\text{\_ref}}}}\right)}^b{\displaystyle \frac{C_j}{C_{\mathrm{ref}}}},\end{array}$$

where R_j (mm d⁻¹) is the total surface runoff on day j simulated by the hydrological module or ORCHIDEE-MUSLE at 0.5^∘ spatial resolution every 30 min. R_{30_j} (mm 30 min⁻¹) is the maximum value of the 48 half-hour runoffs in each day. C_j (0–1, unitless) is the daily actual cover management factor calculated based on the fraction of surface vegetation cover, the amount of litter stock, and the biomass of living roots in each PFT within each $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cell. R_ref, R_{30_ref}, C_ref, and P_ref are the reference values used to estimate the reference sediment delivery rates as describe above.

Daily POC delivery to the river headstream in each 0.5^∘ grid cell is finally simulated based on the sediment delivery rate and the average SOC concentration of surface soil layers (0–20 cm). We assumed that litter cannot be eroded and transported to the river network; however, it can affect the soil erosion rate through the cover management factor of the MUSLE model (denoted by C_j, Eq. 4). The vertical SOC profile is updated every day based on the average depth of eroded soil for each PFT in each 0.5^∘ grid cell of ORCHIDEE. For a more detailed description of ORCHIDEE-MUSLE, we refer to Zhang et al. (2020).

2.2 Sediment and POC transport in inland water network

By merging the model branches ORCHILEAK and ORCHIDEE-MUSLE, the new branch ORCHIDEE-C_lateral combines the novel features of both sources (DOC and POC) described above. The development of ORCHIDEE-C_lateral is complemented by a representation of the sediment and POC transport through the river network that is completely novel and described below.

2.2.1 Sediment transport

Simulation of sediment transport through the river network basically follows the routing scheme of surface water and DOC of ORCHILEAK (Fig. 2). Along with surface runoff (F_{RO_h2o}, m³ d⁻¹), the sediment delivery (F_{RO_sed}, g d⁻¹) from uplands in each basin (i.e., each 0.5^∘ grid cell in the case of this study) initially feeds an aboveground water reservoir (S_{fast_h2o}, m³) with a so-called fast water residence time. From this fast water reservoir, a delayed outflow feeds into the so-called stream reservoir (S_riv, m³, Fig. 2b). Daily water (F_{Fout_h2o}, m³ d⁻¹) and sediment (F_{Fout_sed}, g d⁻¹) flows from the fast water reservoir to the stream reservoir are calculated from a grid-cell-specific topographic index f_topo (unitless, Vörösmarty et al., 2000) extracted from a forcing file (Table 1) and a reservoir-specific factor τ, which translates f_topo into a water residence time of each reservoir (Eqs. 5, 6). Following Guimberteau et al. (2012), the τ of the fast water reservoir (τ_fast) is set to 3.0 d. As the sediment delivery calculated from MUSLE is the net soil loss from headwater basins (gross soil erosion – soil deposition within headwater basins), we assumed that there is no sediment deposition in the fast reservoir and that all of the sediment in the fast reservoir enters the stream reservoir. In addition, only the surface runoff causes soil erosion. The belowground drainage (F_{DR_h2o}, m³ d⁻¹) only transports DOC and dissolved CO₂ to the stream reservoir (Fig. 2b).

$$\begin{array}{}\text{(5)}& {\displaystyle }{F}_{\text{Fout\_h2o}}={\displaystyle \frac{S_{\text{fast\_h2o}}}{{\mathit{\tau }}_{\mathrm{fast}}{f}_{\mathrm{topo}}}}\text{(6)}& {\displaystyle }{F}_{\text{Fout\_sed}}={\displaystyle \frac{S_{\text{fast\_sed}}}{{\mathit{\tau }}_{\mathrm{fast}}{f}_{\mathrm{topo}}}}\end{array}$$

The budget of the suspended sediment in the stream (S_{riv_sed}, g) is determined by F_{out_sed}, the upstream sediment input (F_{up2riv_sed}, g d⁻¹), the sediment input by flooding water returning to the river (F_{fld2riv_sed}, g d⁻¹), the re-detachment of the previously deposited sediment in the riverbed (F_{rero_sed}, g d⁻¹), the bank erosion (F_{bero_sed}, g d⁻¹), the sediment deposition in the riverbed (F_{rd_sed}, g d⁻¹), and the sediment transported to downstream river stretches (F_{down2riv_sed}, g d⁻¹) and, occasionally, floodplains (F_{riv2fld_sed}, g d⁻¹) (Eq. 7).

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}_{\text{riv\_sed}}}{\mathrm{d}t}}=F_{\text{Fout\_sed}}+F_{\text{up2riv\_sed}}+F_{\text{fld2riv\_sed}}+F_{\text{rero\_sed}}\\ & +F_{\text{bero\_sed}}-F_{\text{rd\_sed}}-F_{\text{down2riv\_sed}}-F_{\text{riv2fld\_sed}}\end{array}\end{array}$$

Sediment transport capacity (TC, g m ⁻³) is defined as the maximum concentration of suspended sediment that a given flow rate can carry. TC and the flow rate determine the amount of sediment that can be transported to the downstream grid cell (e.g., F_{down2riv_sed}, F_{riv2fld_sed}). Suspended sediment loads that are in excess of the maximum possible amount of transported sediment will deposit on the riverbed (F_{rd_sed}). If sediment loads are below that maximum possible amount, erosion of the riverbed (F_{rero_sed}) or riverbank (F_{bero_sed}) takes place (Arnold et al., 1995; Nearing et al., 1989; Neitsch et al., 2011).

Figure 2 Simulated lateral transfer processes of water, sediment, and carbon (POC, DOC, and CO₂) in ORCHIDEE-C_lateral (a) and a schematic plot for the reservoirs and flows of water, sediment, and carbon represented in the routing module of ORCHIDEE-C_lateral (b). S_soil is the soil pool. S_rivbed is the sediment (also POC) deposited on the riverbed. S_fast, S_slow, S_riv, and S_fld are the “fast”, “slow”, stream, and flooding water reservoir, respectively. F_RO and F_DR are the surface runoff and belowground drainage, respectively. F_Fout and F_Sout are the flows from the fast and slow reservoir to the stream reservoir, respectively. F_up2riv and F_down2riv are the upstream inputs and downstream outputs, respectively. F_riv2fld is the output from the river stream to the flooding reservoir. F_fld2riv is the return flow from the flooding reservoir to the stream reservoir. F_bed2fld is the transformation from deposited sediment in a riverbed to floodplain soil. F_bero is bank erosion. F_rd and F_rero are the deposition and re-detachment of sediment and POC in a river channel, respectively. F_sub is the flux of DOC and CO₂ from floodplain soil (originating from the decomposition of submerged litter and soil carbon) to the overlying flooding water. F_fd is the deposition of sediment and POC as well as the infiltration of water and DOC. F_D is the wet and dry deposition of DOC from the atmosphere and plant canopy. DOC_l and DOC_r are the labile and refractory DOC pool, respectively. POC_a, POC_s, and POC_p are the active, slow, and passive POC pool, respectively.

In this study, we used an empirical equation adapted from the WBMsed model, which has been proven effective in simulating the suspended sediment discharges in global large rivers (Cohen et al., 2014), to estimate the TC of suspended sediment concentration in streamflow (g m⁻³):

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathrm{TC}={\displaystyle \frac{\mathit{\omega }{q}_{\mathrm{ave}}^{\mathrm{0.3}}{A}^{\mathrm{0.5}}{\left(\frac{q_{i\mathrm{day}}}{q_{\mathrm{ave}}}\right)}^{e_{\mathrm{1}}}(\mathrm{24}\times \mathrm{60}\times \mathrm{60})}{F_{\text{down2riv\_h20}}}},\text{(9)}& {\displaystyle }{e}_{\mathrm{1}}=\mathrm{1.5}-\max \left(\mathrm{0.8},\mathrm{0.145}{\log }_{\mathrm{10}}\mathrm{DA}\right),\end{array}$$

where ω is the coefficient of proportionality, q_ave (m³ s⁻¹) is the long-term average streamflow rate obtained from a historical simulation by ORCHILEAK (Table 1), q_j (m³ s⁻¹) is streamflow rate on day j, e₁ is an exponent depending on the upstream drainage area (DA, m²), and F_{down2riv_h20} (m³ d⁻¹) is the daily downstream water discharge from the stream reservoir. In the stream reservoir of each basin, net deposition occurs when TC is smaller than the concentration of suspended sediment, and the daily deposited sediment (F_{rd_sed}, g d⁻¹) is calculated based on the surplus of the suspended sediment:

$$\begin{array}{}\text{(10)}& F_{\text{rd\_sed}}=c_{\text{rivdep}}\left(S_{\text{riv\_sed}}-\mathrm{TC}{S}_{\text{riv\_h2o}}\right),\end{array}$$

where c_rivdep (0–1, unitless) is the daily deposited fraction of the sediment surplus. Net erosion of the previously deposited sediment in the riverbed (S_{rivbed_sed}, Fig. 2) or riverbank occurs when TC is larger than the concentration of suspended sediment. We assumed that the erosion of the riverbank occurs only after all of the S_{rivbed_sed} has been eroded. Thus, the daily erosion rate (F_{rero_sed}, g d⁻¹) in a river channel is calculated as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& F_{\text{rero\_sed}}=\\ & \left\{\begin{array}{l}{c}_{\mathrm{ebed}}\left(\mathrm{TC}{S}_{\text{riv\_h2o}}-S_{\text{riv\_sed}}\right),\\ c_{\mathrm{ebed}}\left(\mathrm{TC}{S}_{\text{riv\_h2o}}-S_{\text{riv\_ sed}}\right)\le S_{\text{rivbed\_sed}}\\ S_{\text{rivbed\_sed}}+c_{\mathrm{ebank}}\left(\mathrm{TC}{S}_{\text{riv\_h2o}}-S_{\text{riv\_sed}}-S_{\text{rivbed\_sed}}\right),\\ c_{\mathrm{ebed}}\left(\mathrm{TC}{S}_{\text{riv\_h2o}}-S_{\text{riv\_sed}}\right)>S_{\text{rivbed\_sed}}\end{array}\right.\end{array}\end{array}$$

where c_ebed (0–1, unitless) and c_ebank (0–1, unitless) are the fraction of sediment deficit that can be complemented by erosion of the riverbed and riverbank, respectively. After updating the S_{riv_sed} based on the F_{rd_sed} or F_{rero_sed}, the sediment discharge to the downstream basin (F_{down2riv_sed}, g d⁻¹) is calculated based on the ratio of downstream water discharge to the total stream reservoir.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& F_{\text{down2riv\_sed}}=\\ & \left(S_{\text{riv\_sed}}-F_{\text{rd\_sed}}+F_{\text{rero\_sed}}\right){\displaystyle \frac{F_{\text{down2riv\_h2o}}}{S_{\text{riv\_sh2o}}}}\end{array}\end{array}$$

In each basin, the bankfull flow occurs when S_{riv_h2o} exceeds the maximum water storage of the river channel (S_rivmax, g), which is defined by a forcing file (Table 1). Sediment flow from the stream to the floodplain (F_{riv2fld_sed}, g d⁻¹) follows the flooding water, and it is calculated as

$$\begin{array}{}\text{(13)}& {\displaystyle }{F}_{\text{riv2fld\_sed}}=\left(S_{\text{riv\_sed}}-F_{\text{rd\_sed}}+F_{\text{rero\_sed}}\right){\displaystyle \frac{F_{\text{riv2fld\_h2o}}}{S_{\text{riv\_sh2o}}}}\text{(14)}& {\displaystyle }\begin{array}{rl}{F}_{\text{riv2fld\_h2o}}& =\left(S_{\text{riv\_h2o}}-F_{\text{down2riv\_ h2o}}-S_{\text{rivmax}}\right),\\ & {\displaystyle \frac{f_{\text{A\_fld}}}{f_{\text{A\_fld}}+f_{\text{A\_riv}}}},\end{array}\end{array}$$

where f_{A_fld} (0–1, unitless) and f_{A_riv} (0–1, unitless) represent the fraction of floodplain area and river surface area in each basin, respectively. Following the routing scheme of ORCHILEAK, the bankfull flow of a specific basin is assumed to enter the floodplain in the neighboring downstream basin instead of the basin where it originates.

The sediment balance in a flooding reservoir (S_{fld_sed}, g) is controlled by sediment input from the upstream basins (F_{riv2fld_sed}, g d⁻¹), the sediment flowing back to the stream reservoir (F_{fld2riv_sed}, g d⁻¹), and the sediment deposition (F_{fd_sed}, g d⁻¹) (Fig. 2).

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{d}{S}_{\text{fld\_sed}}}{\mathrm{d}t}}=F_{\text{riv2fld\_sed}}-F_{\text{fld2riv\_ sed}}-F_{\text{fd\_sed}}\end{array}$$

Sediment deposition in a floodplain is calculated as the sum of a natural deposition and the deposition due to evaporation (E_h2o, m³ d⁻¹) and infiltration (I_h2o, m³ d⁻¹) of the flooding waters:

$$\begin{array}{}\text{(16)}& F_{\text{fd\_sed}}=c_{\text{flddep}}{S}_{\text{fld\_sed}}+S_{\text{fld\_sed}}{\displaystyle \frac{E_{\text{h2o}}+I_{\text{h2o}}}{S_{\text{fld\_h2o}}}},\end{array}$$

where c_flddep (0–1, unitless) is the daily deposited fraction of the suspended sediment in flooding waters. After removing the deposited sediment from S_{fld_sed}, F_{fld2riv_sed} is calculated based on the ratio of ratio of F_{fld2riv_h2o} to the total flooding reservoir:

$$\begin{array}{}\text{(17)}& {\displaystyle }{F}_{\text{fld2riv\_sed}}=S_{\text{fld\_sed}}{\displaystyle \frac{F_{\text{fld2riv\_h2o}}}{S_{\text{fld\_h2o}}-E_{\text{h2o}}-I_{\text{h2o}}}},\text{(18)}& {\displaystyle }{F}_{\text{fld2riv\_h2o}}={\displaystyle \frac{S_{\text{fld\_h2o}}-E_{\text{h2o}}-I_{\text{h2o}}}{{\mathit{\tau }}_{\text{flood}}{f}_{\mathrm{topo}}}},\end{array}$$

where τ_flood is a factor which translates f_topo (Table 1) into a water residence time of the flooding reservoir. The same as ORCHILEAK, it is set to 1.4 (d m⁻²) in this study.

Note that as the upland soil in ORCHIDEE is composed of clay, silt, and sand particles, the dynamics of clay, silt, and sand sediment in inland waters are simulated separately. To represent the selective transport of clay, silt, and sand sediment, the model parameters ω (Eq. 8) and c_rivdep (Eq. 10) are set to different values when calculating the sediment transport capacity and the deposition of surplus suspended sediment for different particle sizes (Table 2). Moreover, as our model mainly aims to simulate the lateral transfer of sediment and carbon at the decadal to centennial timescale, rather than covering the past thousands of years or even longer time periods, we did not consider the evolution and diversion of river channels in our study.

2.2.2 POC transport and decomposition

Many studies described the selective transport of POC and sediment of different particle sizes. The enrichment ratio (defined as the ratios of the fraction of any given component in the transported sediment to that in the eroded soils) of POC in the transported sediment generally showed significant positive correlation with the fine sediment particles (e.g., fine silt and clay) but negative correlation with the coarse sediment particles (Galy et al., 2008; Haregeweyn et al., 2008; Nadeu et al., 2011; Nie et al., 2015). In ORCHIDEE-C_lateral, the physical movements of POC in inland water systems are simply assumed to follow the flows of the finest clay sediment (Fig. 2b). For example, the fractions of riverine suspended POC, which is deposited on the riverbed (F_{rd_POC}, g C d⁻¹) or transported to the river channel (F_{down2riv_POC}, g C d⁻¹) or floodplain (F_{riv2fld_POC}, g C d⁻¹), are assumed to be equal to the corresponding fractions of clay sediment (Eqs. 19–21). Also, flows of suspended POC in flooding waters to floodplain soil (F_{fd_POC}, g C d⁻¹) or back to the stream reservoir (F_{fld2riv_POC}, g C d⁻¹), as well as the resuspension of POC from the riverbed (F_{rero_POC}, g C d⁻¹), are scaled to the simulated flows of clay sediment (Eqs. 22–24). Note that, similar to SOC, the POC in aquatic reservoirs is divided into three pools: the active (POC_a), slow (POC_s), and passive pool (POC_p) (Fig. 2a). The eroded active, slow, and passive SOC flows into the corresponding POC pools in the “fast” water reservoir (Fig. 2b).

$$\begin{array}{}\text{(19)}& {\displaystyle }{F}_{\text{rd\_POC}}=S_{\text{riv\_POC}}{\displaystyle \frac{F_{\text{rd\_sed\_clay}}}{S_{\text{riv\_sed\_clay}}}}\text{(20)}& {\displaystyle }{F}_{\text{down2riv\_POC}}=S_{\text{riv\_POC}}{\displaystyle \frac{F_{\text{down2riv\_sed\_clay}}}{S_{\text{riv\_sed\_clay}}}}\text{(21)}& {\displaystyle }{F}_{\text{riv2fld\_POC}}=S_{\text{riv\_POC}}{\displaystyle \frac{F_{\text{riv2fld\_sed\_clay}}}{S_{\text{riv\_sed\_clay}}}}\text{(22)}& {\displaystyle }{F}_{\text{fd\_POC}}=S_{\text{fld\_POC}}{\displaystyle \frac{F_{\text{fd\_sed\_clay}}}{S_{\text{fld\_sed\_clay}}}}\text{(23)}& {\displaystyle }{F}_{\text{fld2riv\_POC}}=S_{\text{fld\_POC}}{\displaystyle \frac{F_{\text{fld2riv\_sed\_clay}}}{S_{\text{fld\_sed\_clay}}}}\text{(24)}& {\displaystyle }{F}_{\text{bed2fld\_POC}}=S_{\text{rivbed\_POC}}{\displaystyle \frac{F_{\text{bed2fld\_sed}}}{S_{\text{rivbed\_sed}}}}\end{array}$$

The representation of POC deposition and transformation in the aquatic reservoirs and bed sediment also involves decomposition, which largely follows the scheme used for SOC (Fig. 2a). However, instead of using the rate modifiers for soil temperature and moisture used in the soil carbon module, daily decomposition rates (F_{POC_i}, g C d⁻¹) of each POC pool (S_{POC_i}, g C) are simulated to vary with water temperature based on the Arrhenius term, which is used to simulate the DOC decomposition in ORCHILEAK (Hanson et al., 2011; Lauerwald et al., 2017):

$$\begin{array}{}\text{(25)}& F_{\text{POC\_i}}=S_{\text{POC\_i}}{\displaystyle \frac{{\mathrm{1.073}}^{\left(T_{\text{water}}-\mathrm{28.0}\right)}}{{\mathit{\tau }}_{\text{poc\_i}}}},\end{array}$$

where T_water (^∘C) is the temperature of water reservoirs and is calculated from local soil temperature using an empirical function (Lauerwald et al., 2017). For the POC stored in bed sediment, the temperature of the stream reservoir is used to calculate the decomposition rate. τ_{POC_i} is the turnover time of the i (active, slow, and passive) POC pool. We assumed that the base turnover times of active (0.3 year) and slow (1.12 years) POC pools are the same as for the corresponding SOC pools. The passive SOC pool is generally regarded as the SOC associated with soil minerals or enclosed in soil aggregates (Parton et al., 1987). During the soil erosion and sediment transport processes, the aggregates break down and the passive POC loses its physical protection from decomposition (Chaplot et al., 2005; Z. Wang et al., 2014; Polyakov and Lal, 2008; X. Wang et al., 2014). To represent the acceleration of passive POC decomposition due to aggregate breakdown, we assume that the turnover time of the passive POC is the same as the active POC (0.3 year) rather than the passive SOC (462 years). Similar to the scheme used to simulate SOC decomposition in ORCHILEAK, the decomposed POC from the active, slow, and passive pools flows to other POC pools, flows to DOC pools, or is released to the atmosphere as CO₂ (Fig. 2). Fractions of the decomposed POC flowing to different POC and DOC pools or to the atmosphere are set to the same values used in ORCHILEAK for simulating the fates of the decomposed SOC pools.

Changes in the vertical SOC profile of floodplain soils following sediment deposition are simulated at the end of every daily modeling time step, after physical transfers and decomposition of POC have been calculated. The sediment deposited on the floodplain becomes part of the surface soil layer, and the active, slow, and passive POC flows into the active, slow, and passive SOC pools in surface soil layer, respectively. SOC in the original surface and subsurface soil layers is sequentially transferred to the adjacent deeper soil layers. As the vertical soil profile in ORCHILEAK is described by an 11-layer discretization of a 2 m soil column, we introduce a deep (>2 m) soil pool (S_deep) to represent the soil and carbon transferred down from the 11th soil layer following ongoing floodplain deposition. Decomposition rates of the organic carbon in this deep soil pool are assumed to be the same as those in the 11th (deepest) soil layer. Note that when the soil erosion rate of the floodplain soil is larger than the sediment deposition rate, sediment and organic carbon in S_deep move up to replenish the stocks of the 11th soil layer.

2.3 Model application and evaluation

In this study, ORCHIDEE-C_lateral was applied over Europe and parts of the Middle East (−30^∘ W–70^∘ E, 34–75^∘ N), where extensive observation datasets are available to calibrate and evaluate our model (Table 1). The return period of daily bankfull flow (P_flooding, year), which represents the average interval between two flooding events, is used in this study to produce the forcing file of S_rivmax from a pre-run of ORCHILEAK. Note that P_floodingis generally shorter than the return period of real flooding events, as the flooding may occur on several continuous days and all the flooding waters occurring on these continuous days are generally regarded as belonging to the same flooding event (Fig. S2 in the Supplement). To our knowledge, existing observational data on P_flooding are still very limited. Therefore, following Schneider et al. (2011), we also use a constant P_flooding to simulate the bankfull flows from European rivers and the observed long-term (1961–2000) average bankfull flow rate (m³ s⁻¹) at 66 sites obtained from Schneider et al. (2011) to calibrate P_flooding (the optimized value is 0.1 year, Table 2). To our knowledge, there are still no large-scale observation data on the sediment delivery rates from land to river networks in Europe. Therefore, following Zhang et al. (2020), the parameters a, b, c, and d in Eqs. (1) and (2) (Table 2) were calibrated for 57 European catchments (Fig. S3d) against the modeled sediment delivery data obtained from the European Soil Data Centre (ESDAC, Borrelli et al., 2018). The sediment delivery data from the ESDAC product were derived from WaTEM/SEDEM model simulations using high-resolution data for topography, soil erodibility, land cover, and rainfall. This model was calibrated and validated using observed sediment fluxes from 24 European catchments (Borrelli et al., 2018).

Parameters controlling sediment transport, deposition, and re-detachment (i.e., ω, c_rivdep, c_flddep, c_ebed, and c_ebank, Table 2) in stream and flooding reservoirs were calibrated against the observed long-term averaged sediment discharge rate (Table 1). We also conducted an analysis to test the sensitivity of the simulated riverine sediment and carbon discharges to these parameters following the method used in Tian et al. (2015). The sensitivity of simulation results was evaluated based on the relative changes in simulated riverine sediment and carbon discharges to a 10 % increase and decrease in each parameter (Table 2). The result of the sensitivity analysis shows that the simulated riverine sediment and POC discharges are most sensitive to c_rivdep in Eq. (10), followed by ω in Eq. (8) (Fig. S4). Compared to c_rivdep and ω, the simulated riverine sediment and POC discharges are less sensitive to c_flddep, c_ebed, and c_ebank. With 10 % changes in c_flddep, c_ebed, or c_ebank, the changes in riverine sediment and POC discharges are generally less than 3 %. In addition, the changes in simulated riverine DOC and CO₂ discharges are mostly less than 1 % with 10 % changes in ω, c_flddep, c_ebed, and c_ebank. Nonetheless, a 10 % change in c_rivdep can lead to a change of about 5 % in the simulated riverine CO₂ discharge (Fig. S4).

After parameter calibration, ORCHIDEE-C_lateral was applied to simulate the lateral transfers of water, sediment. and organic carbon in European rivers over the period 1901–2014. Before this historical simulation, ORCHIDEE-C_lateral was run over 10 000 years (spin-up) until the soil carbon pools reached a steady state. In the “spin-up” simulation, the PFT maps, atmospheric CO₂ concentrations, and meteorological data during 1901–1910 were used repeatedly as forcing data. The final simulated water discharge rates in European rivers were evaluated using observation data at 93 gauging sites (locations see Fig. S3a) from the Global Runoff Database (GRDC, Table 1). The simulated bankfull flows were evaluated against observed long-term (1961–2000) average bankfull flows at 66 sites (Fig. S3b) from Schneider et al. (2011). The simulated riverine sediment discharge rate is evaluated using observation data from the European Environment Agency and existing publications (see Table 1) at 221 gauging sites (Fig. S3c). The riverine total organic carbon (TOC), POC, and DOC concentrations provided by the GLObal RIver Chemistry Database (GLORICH, Hartmann et al., 2019) at 346 sites (Fig. S3d) were used to evaluate the simulated riverine POC and DOC concentrations. Note that observations in the GLORICH database, which are measured at gauging sites with drainage area <1.0×10⁴ km², were excluded from our model evaluation because these small catchments cannot be represented by the coarse river network scheme at 0.5^∘ (ca. 55 km at the Equator). Among the retained 346 gauging sites, TOC concentrations were measured at 188 sites, and DOC was measured at 314 sites. POC was measured at only two sites (Bad Honnef with 51 measurements and Bimmen with 78 measurements) in the Rhine catchment and one site (Rheine, 36 measurements) in the Ems catchment (Fig. S3d).

3 Results and discussion

3.1 Model evaluation

3.1.1 Stream water discharge and bankfull flow

Evaluation of our simulation results using in situ observation data from Europe rivers indicates that ORCHIDEE-C_lateral reproduces the magnitude and interannual variation of water discharge rates well in major European rivers (Figs. 3a and S5). Overall, the simulated riverine water discharge rate explained 94 % (Fig. 3a) of the spatial variation of the observed long-term average water discharge rates across 93 gauging sites in Europe (Fig. S3a). Relative biases (calculated as $\frac{\text{simulation-observation}}{\text{observation}}\times \mathrm{100}$ %, as used through the paper if not otherwise stated) of the simulated average water discharge rates compared to the observations are mostly smaller than 30 % (Fig. 3a). For major European rivers, such as the Rhine, Danube, Elbe, Rhone, and Volga, ORCHIDEE-C_lateral also captures the interannual variation of the water discharge rate (Fig. S5). We recognize that ORCHIDEE-C_lateral may overestimate or underestimate the water discharge rate in some rivers (Fig. 3a), particularly in smaller rivers where discrepancy between the stream routing scheme (delineation of catchment boundaries) extracted from the forcing data at 0.5^∘ resolution and the real river network (Fig. S6) can be substantial. An overestimation or underestimation of the catchment area by the forcing data as respectively found for the Elbe and Rhine will introduce a proportional bias in the average amount of simulated discharge from these catchments. Another problem are stream channel bifurcations, which occur in reality but are not represented in a stream network derived from a digital elevation model. For example, in the Danube river delta, a fraction of the discharge is actually exported to the sea through the Saint George Branch, in addition to the water discharge through the main river channel (Fig. S6b). This explains why the simulated water discharge rate at the outlet of the Danube catchment is larger than the observation at the Ceatal gauging station, Romania (identification number in the GRDC database is 6742900, Fig. S5m), where only the main stream discharge was measured.

Figure 3 Comparison between observed and simulated riverine water discharge rates (a) and bankfull flow rates (b). In panel (a), the error bar denotes the standard deviation of interannual variation. Sources of the observed riverine water discharge rate and bankfull flow rate can be found in Table 1.

With the calibrated return period (0.1 year) of the daily flooding rate (see Sect. 2.3), the simulated bankfull flow rates compare well to observations at the 66 sites for which data were available (Fig. 3b). Overall, the simulation result explained 75 % of the inter-site variation of the observed bankfull flow rates. Relative biases of the simulated bankfull flow rates are generally lower than 30 %, although the relative bias may be larger than 100 % at some sites.

3.1.2 Sediment transport

The simulated area-averaged sediment delivery rates from upland to the river network by ORCHIDEE-C_lateral are overall comparable to those simulated by WaTEM/SEDEM for most catchments in Europe (Figs. 4a and S3d). In the two catchments in the Apennine Peninsula, ORCHIDEE-C_lateral gives a drastically lower estimation of the sediment delivery rates compared to WaTEM/SEDEM. By excluding these two catchments, ORCHIDEE-C_lateral reproduces 72 % of the spatial variation of the sediment delivery rates estimated by WaTEM/SEDEM (Fig. 4a). In addition, the average sediment loss rate over all catchments shown in Fig. S3d is 40.8 g m⁻² yr⁻¹, which is overall comparable to the estimate by WaTEM/SEDEM (42.5 g m⁻² yr⁻¹).

Figure 4 Comparison between the simulated area-averaged sediment delivery rate from uplands to the river network from ORCHIDEE-C_lateral and WaTEM/SEDEM (a), as well as the comparison between observed and simulated annual sediment discharge rates at 221 gauging sites (b). In panel (a), the two hollow dots represent the sediment delivery rates at the two catchments in the Apennine Peninsula (Fig. S3d). The regression function in panel (a) was obtained based on the values of all solid dots, excluding the two hollow dots. In panel (b), the error bar denotes the standard deviation of interannual variation. Sources of the observed annual sediment discharge rate are in Table 1.

ORCHIDEE-C_lateral reproduces 83 % of the inter-site variation of the observed riverine sediment discharge rates across Europe (Fig. 4b). Simulation of the riverine sediment discharge rate at large spatial scale is still a big challenge. It generally needs detailed information on the streamflow, geomorphic properties of river channel, and the particle composition of the suspended sediment (Neitsch et al., 2011). Moreover, the parameters of existing sediment transport models usually require recalibration when they are applied to different catchments (Gassman et al., 2014; Oeurng et al., 2011; Vigiak et al., 2017). In ORCHIDEE-C_lateral, the sediment processes in river networks are simulated using simple empirical functions and parameters based on a routing scheme at a spatial resolution of 0.5^∘ (Sect. 2.2.1). Detailed information about the streamflow (e.g., cross-sectional area) and the geomorphic properties of river channels is not represented. Sediment discharge in all catchments was simulated using a universal parameter set. This may explain why ORCHIDEE-C_lateral fails to capture the observed sediment discharge rates in some specific catchments, especially those with relatively small drainage areas (e.g., <5×10³ km²).

3.1.3 Organic carbon transport

Simulation of the riverine carbon discharge rate at large spatial scale is an even bigger challenge than simulating sediment discharge, as the riverine carbon discharge is controlled by many factors, such as upland topsoil SOC concentrations, soil erosion rate, the transport and deposition rate of the clay fraction in the river channel and on the floodplain, and the decomposition of POC in transit and in aquatic sediments. As described above, the simulated water discharge rate, bankfull flow, and sediment discharge rate are overall comparable to observations (Figs. 3 and 4). The simulated total SOC stock in the top 0–30 cm soil layer in Europe of 107 Pg C is close to the value extracted from the HWSD database (106 Pg C) but significantly lower than the values extracted from some other databases, such as GSDE (249 Pg C), SoilGrids (202 Pg C), S2017 (148 Pg C) and landGIS (226 Pg C) (Fig. S7a). We noticed that the SOC stocks extracted from these observation-based soil databases show considerable difference (vary from 106 to 249 Pg C), as they have been produced using different clusters of site-level SOC measurements and different interpolation methods to produce global gridded SOC stocks from the site-level measurements (Shangguan et al., 2014; Hengl et al., 2014; Sanderman et al., 2017). Distributions of the simulated SOC stock along the latitude gradients (30–75^∘ N) are overall comparable to those extracted from the HWSD and S2017 databases (Fig. S7). But even compared to these two databases, our model still underestimated the SOC stock in southern Europe (30–41^∘ N).

Comparison of the simulated concentrations of riverine organic carbon and the observations obtained from the GLORICH database (Hartmann et al., 2019) indicates that our model can basically capture the TOC and DOC concentrations in European rivers (Figs. 5, 6, S8 and S9). The simulation results explain 34 % and 32 % of the inter-site variation of the observed TOC and DOC concentrations, respectively (Fig. 5). For major European rivers, such as the Rhine, Elbe, Danube, Spree, and Weser, the simulated long-term average TOC and DOC concentrations are overall close to the observations (Figs. 6, S8, and S9). But for the Rhone river in southern France, the DOC concentrations have been systematically overestimated by more than 50 % (Figs. 6 and S9m). In addition, both simulated and observed TOC and DOC concentrations show drastic temporal (both seasonal and interannual) variations (Figs. 5, S8 and S9). Our model seems to have overestimated the temporal variation of TOC and especially DOC concentrations (Figs. S8 and S9). Nonetheless, the simulated temporal variations of TOC and DOC discharge rates are overall comparable to the observations (Figs. S10 and S11), as our model can capture the magnitude and temporal variation of riverine water discharge rates well.

Figure 5 Comparison between the observed and simulated riverine TOC (a, POC + DOC) and DOC (b) concentrations. The dot and error bar denote the mean and standard deviation at each gauging site, respectively. Note that the mean and standard deviation of the simulated concentrations at each site are calculated based on the monthly average value, but the mean and standard deviation of the observed concentrations are based on instantaneous observation.

Figure 6 Comparison between the observed and simulated concentrations of total organic carbon (TOC, a) and dissolved organic carbon (DOC, b) in river flows, as well as the discharge rates of riverine TOC and DOC. The black and pink lines in each box denote the median and mean value, respectively. Box boundaries show the 25th and 75th percentiles, whiskers denote the 10th and 90th percentiles, and the dots below and above each box denote the 5th and 95th percentiles, respectively.

In Europe, the GLORICH database only provides POC concentrations measured at three gauging stations in northwestern Germany (Figs. 7, S3d). The simulated POC concentrations and discharge rates in the Ems river at Rheine are overall comparable to the observations (Fig. 7e, f). However, at the two gauging sites at the river Rhine, the POC concentrations have been significantly underestimated (Fig. 7a–d). We noticed that the stream routing scheme of the Rhine catchment at 0.5^∘ obtained from the forcing data STN-30p (Vörösmarty et al., 2000) differs significantly from the stream routing scheme extracted based on the high-resolution (3") DEM (Fig. S6). Thus, besides the errors in simulated SOC stocks, soil erosion rate, stream discharge rate, and sediment transport and deposition rate, the inaccurate stream routing scheme used in this study might also be an important reason for the underestimation of POC concentration in the Rhine river.

Figure 7 Comparison between observed (instantaneous measurements) and simulated (monthly average values) riverine POC concentrations and POC discharge rates at three gauging sites. The histograms and error bars denote the means and standard deviations of POC concentrations, respectively. Long-term average water discharge rates at Bad Honnef, Bimmen and Rheine during the observation periods are 2023, 2100 and 80 m³ s⁻¹, respectively.

3.2 Lateral carbon transfers in Europe

Based on our simulation results, the average annual sediment delivery from upland to the river network caused by water erosion in Europe (−30^∘ W–70^∘ E, 34–75^∘ N) during 1901–2014 is 2.8±0.4 Pg yr⁻¹ (Fig. 8a). From northern to southern Europe, the sediment delivery rate from upland to the river increases from less than 1.0 g m⁻² yr⁻¹ in the Scandinavia Peninsula, which is covered by mature boreal forests (Fig. S12a), and in the Northern European Plain to more than 600 g m⁻² yr⁻¹ in the mountainous regions of the Apennine Peninsula, Balkan Peninsula, and the Middle East (Figs. 9a, S13a). In total across Europe, 63.2 % (1.8±0.2 Pg yr⁻¹) of the sediment delivered into the river network is deposited in river channels and floodplains, and the remaining 36.8 % (1.0±0.1 Pg yr⁻¹) is exported to the sea (Fig. 8a). Generally, large rivers, like the Danube, Volga, and Ob rivers, carry more sediment to the sea than small rivers (Fig. 9b, c). But several relatively small rivers in the Middle East and the Po river in northern Italy also carry a similarly large amount of sediment to the sea, as the upland soil erosion rates are very high (>200 g m⁻² yr⁻¹) in these catchments (Fig. 9a, c). Spatial distribution of the sediment deposition is controlled by the stream routing scheme and the spatial distribution of floodplains (Fig. 10b). In northern and central Europe, the area-averaged sediment deposition rates (i.e., calculated as the amount of annual sediment deposition in each $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cell divided by the grid cell area) in river channels and floodplains are mostly less than 100.0 g m⁻² yr⁻¹ (Fig. 9d). In the downstream part of the Danube, Po, and several rivers in the Middle East, the sediment deposition rate can exceed 800.0 g m⁻² yr⁻¹. From 1901 to the 1960s, the annual total sediment delivery from uplands to the whole river network of Europe declined significantly (p<0.01, independent sample t test) from about 3.0 Pg yr⁻¹ to about 2.3 Pg yr⁻¹ (Fig. S14a). From 1960 to 2014, the annual sediment delivery rate did not show a significant trend but revealed large interannual variations.

Figure 8 Averaged annual lateral redistribution rate of sediment (a), POC (b), DOC (c), and CO₂ (d) in Europe for the period 1901–2014. F_{sub_DOC} and F_{sub_CO2} are the DOC and CO₂ inputs from floodplain soil (originating from the decomposition of submerged litter and soil carbon) to the overlying flooding water, respectively.

Along with soil erosion and sediment transport, the average annual POC delivery from upland to the river network in the whole of Europe during 1901–2014 is 10.1±1.1 Tg C yr⁻¹ (Fig. 8b). 41.0 % of the POC delivered into the river network is deposited in river channels and floodplains, 2.9 % is decomposed during transport, and the remaining 56.1 % is exported to the sea. Spatial patterns of the area-averaged SOC delivery rate and POC discharge rate basically follow that of sediment (Fig. 10a, c). Although the sediment discharge rates in some rivers in the Middle East can be as high as that in the Danube or Volga river (Fig. 9c), the POC delivery rates in these rivers are much smaller than in the larger ones (Fig. 10c). This is mainly due to the lower SOC stocks in the Middle East compared to those found in the Danube and Volga catchments (Fig. S7). We also note that different from the sediment delivery, the annual total POC delivery from upland to the river network in Europe did not show a significant declining trend from 1901 to the 1960s (Fig. S14b). The increase in SOC stock (Fig. S14c) may have partially offset the decline in sediment delivery rate.

Figure 9 Averaged annual lateral redistribution rate of water and sediment in Europe during 1901–2014. (a) Annual sediment delivery rate from upland to the river network; (b) annual water discharge rate; (c) annual sediment discharge rate; (d) annual net sediment budget in each $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cell. In panel (d), the positive and negative values denote net gain and net loss of sediment, respectively.

Leaching results in an average annual DOC input of 13.5±1.5 Tg C yr⁻¹ from soil to the river network in Europe, and the in situ DOC production caused by wet deposition and the decomposition of riverine POC as well as submerged litter and soil organic carbon under flooding waters amounts to 2.2±0.7 Tg C yr⁻¹ (Fig. 8c). 28.1 % of the total riverine DOC then infiltrates the floodplain soils, 12.9 % is decomposed during riverine transport, and the remaining 59.0 % is exported to the sea. The spatial distribution of the DOC leaching rate is very different from that of POC (Fig. 10b). From northwestern Europe to southeastern Europe and the Middle East, the DOC leaching rates decrease from over 6 g C m⁻² yr⁻¹ to less than 1.0 g C m⁻² yr⁻¹. DOC discharge rates in major European rivers, such as the Rhine, Danube, Volga, Elbe, and Ob, are mostly higher than 100 Tg C yr⁻¹ (Fig. 10d). Comparatively, the DOC discharge rates in southern Europe and the Middle East are significantly lower (<60 Tg C yr⁻¹).

Figure 10 Averaged annual lateral redistribution rate of organic carbon in Europe during 1901–2014. (a) Annual SOC delivery rate from upland to the river network; (b) annual DOC leaching rate; (c) annual POC discharge rate; (d) annual DOC discharge rate.

The average annual leaching rate of CO₂ sourced from the decomposition of upland litter and soil organic carbon (incl. DOC) in the whole of Europe is 14.3±2.2 Tg C yr⁻¹ (Fig. 8a). Decomposition of the submerged litter and organic carbon in floodplains and the decomposition of riverine POC and DOC add an in situ CO₂ production amounting to 7.5±2.7 and 4.1±0.5 Tg C yr⁻¹, respectively. Most of this CO₂ (80.2 %) feeding stream waters is then released back to the atmosphere quickly, in such a way that only 15.8 % of the CO₂ is exported to the sea, and 4.0 % infiltrates the floodplain soils.

3.3 Implications for the terrestrial C budget of Europe

Representing the lateral carbon transport in an LSM is helpful to estimate the terrestrial carbon cycle more accurately. From the year 1901 to 2014, soil erosion and leaching combined resulted in a 5.4 Pg loss of terrestrial carbon to the European river network, with this amount corresponding to about 5 % of the total SOC stock (106 Pg C, Fig. S7a) in the 0–30 cm soil layer. The average annual total delivery of organic carbon (POC + DOC) during the same period is 47.3±6.6 Tg C yr⁻¹ (Fig. 8), which is about 4.7 % of the net ecosystem production (NEP at 993±255 Tg C yr⁻¹, defined as the difference between the vegetation primary production – NPP – and the soil heterotrophic respiration – Rh – due to the decomposition of litter and soil organic matter, i.e., NEP = NPP − Rh) and 19.2 % of the net biome production (NBP at 243±189 Tg C yr⁻¹, defined as the difference between NEP and the land carbon loss – Rd – due to additional disturbances, e.g., harvest, land cover change, and soil erosion, and leaching, i.e., NBP = NEP − Rd − DOC and POC to the river) (Fig. 11b). The annual total export of carbon to the sea surrounding Europe is 19.0±1.4 Tg C yr⁻¹, which amounts to 1.9 % and 8.7 % of the NEE and NBP, respectively.

Figure 11 The simulated average annual carbon budget of the terrestrial ecosystem in Europe during the 1901–2014 when the lateral carbon transport is ignored (a) and considered (b). All fluxes are presented as mean ± standard deviation. NPP is the net primary production. Rh and Rd are the heterotrophic respiration and the respiration due to disturbances like harvest and land cover change, respectively. ΔC_land represents the average annual changes in the total land carbon stock. The percentage following each of these changes in blue represents the average annual relative changes in the corresponding carbon pool. C_land2riv, C_riv2land, and C_riv2sea are the average annual carbon fluxes from land to inland waters, from inland waters to floodplains, and from inland waters to the sea, respectively. SD is the standard deviation.

Besides direct transfers of organic carbon from soil to aquatic systems, the lateral transport of water, sediment, and carbon can also affect the land carbon budget in several indirect ways. First, the lateral redistribution of surface runoff can affect the land carbon budget by altering soil wetness. Our simulation results reveal that the lateral redistribution of runoff can significantly change local soil wetness, especially in floodplains (Fig. S13b), where the increase in soil wetness can be larger than 10 % (Fig. S16b). Soil wetness is a key controlling factor of plant photosynthesis (Knapp et al., 2001; Stocker et al., 2019; Xu et al., 2013). Benefiting from the increase in soil wetness, the NPP in many grid cells with a large area of floodplain increased by more than 5 % (Fig. 11b), although the NPP over the whole of Europe only increased by 1 % (Fig. 11). Changes in soil wetness can further alter soil temperature (Fig. S16a). As soil wetness and temperature are the two most important controlling factors of organic matter decomposition, the lateral redistribution of runoff can affect the local land carbon budget by changing the Rh. Moreover, in ORCHIDEE-C_lateral, the turnover times of litter and SOC under flooding waters (assumed to experience anaerobic conditions) are set to be one-third of the litter and SOC turnover times in upland soil (Reddy and Patrick, 1975; Neckles and Neill, 1994; Lauerwald et al., 2017). Accounting for flooding thus decreases the decomposition rate of litter and SOC stored in floodplain soils.

Second, soil erosion and sediment deposition can affect the land carbon budget by altering the vertical distribution of litter and soil organic carbon. At the net erosion sites of the uplands, the loss of surface soil results in some of the belowground litter and SOC that were originally stored in deeper soil layers emerging to the surface soil layers, and it also results in a fraction of the belowground litter becoming aboveground litter. In the floodplains, the newly deposited sediment becomes part of the surface soil layer, and the belowground litter and SOC in the original surface soil layer are transferred down to the deeper soil layers. As the temperatures and fresh organic matter inputs (sourced from the aboveground litterfall and dead roots), which can impact SOC decomposition rates through the priming effect (Guenet et al., 2010, 2016), in different soil layers are different, changes in the vertical distribution of belowground litter and SOC can therefore lead to changes in the overall decomposition rate of the organic matter in the whole soil column.

Third, soil aggregates mostly break down during soil erosion and sediment transport, and the riverine POC thus loses some of its physically protection from decomposition (Hu and Kuhn, 2016; Lal, 2003). Some modeling studies have assumed that at least 20 % of the eroded SOC would be decomposed during soil erosion and transport processes (Lal, 2003, 2004; Zhang et al., 2014). However, the estimation by Smith et al. (2001) using a conceptual mass balance model suggests that only a tiny fraction of the eroded POC is decomposed and released as CO₂ to the atmosphere. Using laboratory rainfall simulation experiments, van Hemelryck et al. (2011) estimated a 2 %–12 % mineralization of the eroded SOC from a loess soil, and X. Wang et al. (2014) estimated a mineralization of only 1.5 %. In ORCHIDEE-C_lateral, the passive SOC pool is regarded as the SOC associated with soil minerals and protected by soil aggregates. The turnover time of the passive POC in the river stream and flooding waters is assumed to be the same as that of the active POC (0.3 year). Our simulation results suggest that the fraction of total riverine POC that is decomposed during lateral transport from uplands to the sea is 2.9 % in Europe (Fig. 8b), which is larger than the POC decomposition fraction (0.9 %) when the turnover time of the passive POC in rivers is assumed to be the same as that of the passive POC (i.e., no soil aggregates break down). The acceleration of POC decomposition rate due to the breakdown of soil aggregates can thus slightly affect the estimate of the regional land–atmosphere carbon flux. Moreover, the riverine POC and DOC can be transported over a long distance and finally settle or infiltrate floodplains or river channels (especially estuarine deltas) where the local environmental conditions might be quite different from those encountered in the uplands from where these C pools originate. These changes in environmental conditions can affect the decomposition rate of the laterally redistributed organic carbon (Abril et al., 2002).

Comparison between the simulation results from ORCHIDEE-C_lateral with activated and deactivated erosion and river routing modules indicates that ignoring lateral carbon transport processes in LSMs may lead to significant biases in the simulated land carbon budget (Figs. 11 and S14). Although the omission of lateral carbon transport in ORCHIDEE-C_lateral only resulted in a 1 % decrease in simulated average annual total NPP in Europe during 1901–2014 and a 1 % increase in annual total Rh, the annual total NBP (NEP–Rd–DOC and POC to the river) is overestimated by 4.5 %. Over the same period, the lateral carbon transport only induced a 0.09 % decrease in the total SOC and DOC stock in Europe (Fig. S15c), but their spatial distribution was significantly altered (Fig. 12e, f). For instance, in some mountainous regions, the soil erosion induced a reduction of the SOC stock by more than 8 %. On the contrary, the sediment and POC deposition in some floodplains led to an increase in SOC stock by more than 8 % (Fig. 12f).

Figure 12 Changes (first column) and relative changes (second column) in the net primary production (NPP), heterotrophic respiration (Rh), and total soil organic carbon (SOC, 0–2 m) in Europe due to lateral carbon transport during 1901–2014. For each variable, the change is calculated as C_lat−C_nolat, where C_lat and C_nolat are the carbon fluxes or stocks when lateral carbon transport is considered and ignored, respectively. The relative changes is calculated as ($C_{\text{lat}}-C_{\text{nolat}})/C_{\text{nolat}}\times \mathrm{100}$ %.

Consistent with previous studies (Stallard, 1998; Smith et al., 2001; Hoffmann et al., 2013), our simulation results reveal the importance of sediment deposition in floodplains for the overall SOC budget. From 1901 to 2014, erosion and leaching over Europe totally induced a loss of 3.03 Pg organic carbon (POC + DOC) from uplands to the river network, and only 0.65 Pg of this carbon was redeposited onto the floodplains. The total stock of soil organic carbon in Europe thus should have decreased by 2.38 Pg C. However, due to the decrease in decomposition rate of the buried organic carbon (including in situ and ex situ carbon) in floodplain soils, the total stock of soil organic carbon in Europe only decreased by 0.91 Pg C. Floodplains in Europe have protected a total of 2.12 (= 3.03–0.91) Pg of soil organic carbon from being transported to the sea or released to the atmosphere in forms of CO₂. Although the sequestration of organic carbon in floodplains cannot make up all of the soil organic carbon (POC + DOC) loss, the increased organic carbon stock in floodplains (2.12 Pg C) is much higher than the soil POC loss (0.86 Pg C) induced by soil erosion.

3.4 Uncertainties and future work

In the present version of ORCHIDEE-C_lateral, the lateral transfers of sediment and carbon are simulated using a simplified scheme due to the fragmented nature of large-scale forcing (e.g., geomorphic properties of the river channel) and validation data (e.g., continuous sediment and carbon concentration data in river streams and deposition–erosion rates in river channels). We recognize that this simplification induces significant uncertainties in model outputs, especially regarding changes in lateral sediment and particulate carbon transfers under climate change and direct human perturbations. Several physics-based algorithms have been proposed to accurately calculate the TC of streamflows (Arnold et al., 1995; Molinas and Wu, 2001; Nearing et al., 1989). These algorithms mostly require detailed information about the stream power (e.g., flow speed and depth), geomorphic properties of the river channel (e.g., slope and hydraulic radius), and the physical properties of the sediment particles (e.g., median grain size) (Neitsch et al., 2011). They are good predictors to estimate TC in rivers with detailed observation data on local stream, soil, and geomorphic properties. Unfortunately, it is not practical to implement those algorithms in ORCHIDEE-C_lateral due to the lack of appropriate forcing data at large scale as well as the relatively rough representation of streamflow dynamics compared to hydrological models for small basins. For example, runoff and sediment from all headwater basins in one 0.5^∘ grid cell of ORCHIDEE-C_lateral are assumed to flow into one single virtual river channel. Although the total river surface area in each grid cell is represented (obtained from forcing file, Table 1, Lauerwald et al., 2015), the length, width, and depth of the river channel are unknown. Furthermore, in reality, there can be multiple river channels in the area represented by each grid cell, and these channels might flow in different directions.

We also noticed that previous studies have derived empirical functions of upstream drainage area (e.g., Luo et al., 2017) or upstream runoff (e.g., Yamazaki et al., 2011) to calculate the river width and depth, allowing simulation of the water flow in the river channel using physically based algorithms. Unfortunately, to obtain a good fit of the simulated river discharges against observations, the parameters in the empirical functions for calculating river width and depth generally need to be calibrated separately for each catchment (Luo et al., 2017), an approach that is incompatible with large-scale simulations like those performed here. Without such calibration, the simulated geometrical properties of the river channel and runoff are prone to large uncertainties, thus rendering the simulation of sediment transport at continental or global scale using physically based algorithms a more challenging task. Given the difficulty of simulating the detailed hydraulic dynamics of streamflow at large spatial scale, we thus apply a simple approach (Eq. 8) to calculate the sediment transport capacity. Overall, we encourage future studies to produce large-scale databases on the geomorphic properties of global river channels (e.g., river depth and width) and to develop large-scale sediment transport models capable of producing more realistic and accurate simulations of sediment deposition, re-detachment, and transport processes, as well as including the exchanges of water, sediment, and carbon between river streams and floodplains.

The simulation of the soil DOC dynamics and leaching in our model needs to be further improved to better simulate the seasonal variation of riverine DOC and TOC concentrations. The concentration of soil DOC and the DOC decomposition rate during the lateral transport process in the river network are the two key factors controlling DOC concentration in river flow. As only a small fraction (<20 %) of the riverine DOC is decomposed during lateral transport (Fig. 8), the overestimated (Fig. 6) seasonal amplitude in riverine DOC (and TOC) concentrations is likely caused by the uncertainties in the simulated seasonal dynamics of the leached soil DOC. The current scheme used in our model for simulating soil DOC dynamics has been calibrated against observed DOC concentrations at several sites in Europe (Camino-Serrano et al., 2018). Although the calibrated model can overall capture the average concentrations of soil DOC, it is not able to fully capture the temporal dynamics of DOC concentrations (Camino-Serrano et al., 2018). Given this, it is necessary to collect additional observation data on the seasonal dynamics of soil DOC concentration to further calibrate the soil DOC model. In addition, averaged over the various DOC and SOC pools we distinguish in the soils, DOC represents a much more reactive fraction of soil carbon (with a turnover time of several days to a few months) than SOC (with a turnover time of decades to thousands of years). Therefore, soil DOC concentrations experience large seasonal variations, while SOC concentrations generally are much more stable and show very limited seasonal dynamics. Overall, seasonal variations in riverine POC concentrations are mainly controlled by the seasonal dynamics of soil erosion rates rather than by the seasonal SOC dynamics, which explains a partial decoupling in the behavior of POC compared to that of DOC.

Although most processes related to lateral carbon transport have been represented in ORCHIDEE-C_lateral, there are still omitted processes and large uncertainties in our model. For example, many studies suggest that a substantial portion of the eroded sediment and carbon is deposited downhill at adjacent lowlands as colluviums rather than being exported to the river (Berhe et al., 2007; Smith et al., 2001; Hoffmann et al., 2013; Wang et al., 2010). As the deposition of sediment and carbon within headwater basins can also significantly alter the vertical SOC profile and soil micro-environments (e.g., soil moisture, aeration, and density) (Doetterl et al., 2016; Gregorich et al., 1998; Wang et al., 2015; Zhang et al., 2016), omission of this process may result in uncertainties in the simulated vegetation production and SOC decomposition. In addition, the impact of artificial dams and reservoirs on riverine sediment and carbon fluxes is also not represented in our model. Construction of dams generally leads to increased water residence time, nutrient retention, and sediment and carbon trapping in the impounded reservoir (Habersack et al., 2016; Maavara et al., 2017), and it can also affect the downstream flooding regime and frequency (Mei et al., 2016; Timpe and Kaplan, 2017). Estimation by Maavara et al. (2017) suggests that the organic carbon trapped or mineralized in global artificial reservoirs is about 13 % of the total organic carbon carried by global rivers to the oceans. To more accurately simulate the lateral carbon transport, we plan to include the soil and carbon redistribution within headwater basins and the effects of dams and reservoirs on riverine sediment and carbon fluxes in our model in the near future.

The effects of lateral redistribution of water and sediment on vegetation productivity has not been fully represented in our model. As shown above, our model is able to represent the impacts of lateral water redistribution on vegetation productivity though modifying local soil wetness (Figs. 12 and S16). However, in addition to modifying soil wetness, many studies have indicated that soil erosion and sediment deposition can affect vegetation productivity by modifying soil nutrient (e.g., nitrogen – N and phosphorus – P) availability (Bakker et al., 2004; Borrelli et al., 2018; Quine, 2002; Quinton et al., 2010). Recently, terrestrial N and P cycles have already been incorporated into another branch of ORCHIDEE (i.e., the ORCHIDEE-CNP developed by Goll et al., 2017). By coupling our new branch and ORCHIDEE-CNP, it will be possible to develop a more comprehensive LSM that can also simulate the effects of lateral N and P redistribution on vegetation productivity.

Although soils are the major source of riverine organic carbon, domestic, agricultural, and industrial waste, as well as river-borne phytoplankton, can also make significant contributions (Abril et al., 2002; Meybeck, 1993; Hoffmann et al., 2020). Moreover, previous studies have shown that sewage generally contains highly labile POC, while most aquatic production is generally mineralized within a short time (Abril et al., 2002; Caffrey et al., 1998). Omission of organic carbon inputs from manure and sewage could potentially lead to an underestimation of CO₂ evasion from the European river network. Inclusion of these additional carbon sources should thus help improve simulation of aquatic CO₂ evasion.

Uncertainties in our simulation results also stem from the forcing data (Table 1) applied in our model. The routing scheme of water, sediment, and carbon is driven by a map of streamflow direction at 0.5^∘ spatial resolution (Guimberteau et al., 2012). Comparison between this flow direction map and the flow direction map derived based on a high-resolution (3") DEM shows discrepancies between the two river flow networks (Fig. S6). As the flow direction directly determines the area of each catchment and the route of river flows, errors in forcing data for flow direction may thus induce uncertainties in the simulated riverine water, sediment, and carbon discharges. Land cover maps are another source of uncertainty. For instance, croplands generally experience significantly larger soil erosion rates than grasslands and forests (Borrelli et al., 2017; Nunes et al., 2011; Zhang et al., 2020). However, croplands in ORCHIDEE are only represented in a simplified way by segmenting them into C₃ and C₄ crops based on their photosynthesis characteristics. Therefore, our simulations based on land cover data with only two broad groups of crops might not be able to fully capture the seasonal dynamics of planting, canopy growth rate, and harvesting for all crop types. Furthermore, the effects of soil conservation practices, which would decrease erosion rates, are ignored in our model. Panagos et al. (2015) have shown that contour farming as well as stone wall and grass margin techniques have been applied in Europe to reduce the risk of soil erosion. However, these soil conservation practices only reduce the average erosion rate in the European Union by 3 %. Excluding soil conservation practices should thus have a limited impact in our simulation results.

Further model calibration, evaluation, and development are necessary for improving our model. Due to the limitation of observation data, we calibrated the parameters controlling sediment transport, deposition, and re-detachment (i.e., ω, c_rivdep, c_flddep, c_ebed, and c_ebank in Table S1) in streams and flooding reservoirs only against the observed sediment yield. Even though our model can overall capture the lateral transfers of sediment and carbon in many rivers in central and northern Europe, more observation data are crucially needed to further evaluate the performance of our model, in particular in southern Europe. In addition, it is still unknown whether our model can satisfactorily simulate intermediate processes such as sediment deposition in river channels and floodplains, as well as the rate of river channel erosion. It is also unknown whether our model would perform satisfactorily in regions with very different climates than Europe such as the tropical region. Thus, in the future, an important aim will be to further calibrate our model against more detailed observation data (e.g., sediment deposition rate in river channels and floodplains) and extend the model application to regions of contrasting climate, vegetation, and topography. Moreover, the GLORICH database (Hartmann et al., 2019) only provides instantaneous observations of riverine organic carbon concentrations, and it is therefore difficult to evaluate the model's ability to reproduce temporal trends. Therefore, future modeling efforts should be combined with data mining efforts targeting the collection of continuous (e.g., daily) and long-term observational data on organic carbon content and fluxes in streams and rivers.

4 Conclusions

By merging ORCHILEAK (Lauerwald et al., 2017) and an upgraded version of ORCHIDEE-MUSLE (Zhang et al., 2020) for the simulation of DOC and POC from land to sea, respectively, we developed ORCHIDEE-C_lateral, a new branch of the ORCHIDEE LSM. ORCHIDEE-C_lateral simulates the large-scale lateral transport of water, sediment, POC, DOC, and CO₂ from uplands to the sea through river networks, the deposition of sediment and POC in river channels and floodplains, the decomposition POC and DOC during fluvial transport, and the CO₂ evasion to the atmosphere, as well as the changes in soil wetness and vertical SOC profiles due to the lateral redistribution of water, sediment, and carbon.

Evaluation using observation data from European rivers indicates that ORCHIDEE-C_lateral can satisfactorily reproduce the observed riverine discharges of water and sediment, bankfull flows, and organic carbon concentrations in river flows. Application of ORCHIDEE-C_lateral to the entire European river network from 1901 to 2014 reveals that the average annual total carbon delivery to streams and rivers amounts to 47.3±6.6 Tg C yr⁻¹, which corresponds to about 4.7 % of total NEP and 19.2 % of the total NBP of terrestrial ecosystems in Europe. The lateral transfer of water, sediment, and carbon can affect the land carbon dynamics through several different mechanisms. Besides directly inducing a spatial redistribution of organic carbon, it can also affect the regional land carbon budget by altering vertical SOC profiles, as well as the soil wetness and soil temperature, which in turn impact vegetation production and the decomposition of soil organic carbon. Overall, omission of lateral carbon transport in ORCHIDEE potentially results in an underestimation of the annual mean NBP in Europe of 4.5 %. In regions experiencing high soil erosion or high sediment deposition rate, lateral carbon transport also changes the total SOC stock significantly, by more than 8 %.

We recognize that ORCHIDEE-C_lateral still has several limitations and significant uncertainties. To address those, we plan to enhance our model with additional processes, such as sediment deposition at downhill areas or the regulation of lateral transport by dams and reservoirs. We also plan to calibrate and further evaluate our model by extending the observational dataset to regions outside Europe.