Increased Terrestrial Carbon Export and CO2 Evasion From Global Inland Waters Since the Preindustrial Era

Global carbon dioxide (CO2) evasion from inland waters (rivers, lakes, and reservoirs) and carbon (C) export from land to oceans constitute critical terms in the global C budget. However, the magnitudes, spatiotemporal patterns, and underlying mechanisms of these fluxes are poorly constrained. Here, we used a coupled terrestrial–aquatic model to assess how multiple changes in climate, land use, atmospheric CO2 concentration, nitrogen (N) deposition, N fertilizer and manure applications have affected global CO2 evasion and riverine C export along the terrestrial‐aquatic continuum. We estimate that terrestrial C loadings, riverine C export, and CO2 evasion in the preindustrial period (1800s) were 1,820 ± 507 (mean ± standard deviation), 765 ± 132, and 841 ± 190 Tg C yr−1, respectively. During 1800–2019, multifactorial global changes caused an increase of 25% (461 Tg C yr−1) in terrestrial C loadings, reaching 2,281 Tg C yr−1 in the 2010s, with 23% (104 Tg C yr−1) of this increase exported to the ocean and 59% (273 Tg C yr−1) being emitted to the atmosphere. Our results showed that global inland water recycles and exports nearly half of the net land C sink into the atmosphere and oceans, highlighting the important role of inland waters in the global C balance, an amount that should be taken into account in future C budgets. Our analysis supports the view that a major feature of the global C cycle–the transfer from land to ocean–has undergone a dramatic change over the last two centuries as a result of human activities.

• Terrestrial carbon loading since 1800 has increased by 25%, with 23% of this increase exported to the ocean and 59% being emitted to the atmosphere • Atmospheric CO 2 and N inputs dominated C export increase, while the climate and land use change dominated the decadal variations in CO 2 evasion • Global inland water recycles and exports nearly half of the net land C sink into the atmosphere and oceans in the 2010s Supporting Information: Supporting Information may be found in the online version of this article.

Introduction
Inland waters (streams, rivers, lakes, and reservoirs), though occupying only 1% of the Earth surface, represent a vital link connecting two of the largest active carbon (C) pools, terrestrial, and marine ecosystems.Inland waters also serve as an important conduit for the exchange of CO 2 with the atmosphere.Terrestrial C cycle changes can drastically alter the patterns of C loading (C transfer from land to inland waters), which further influences C dynamics in coastal and marine ecosystems (Butman et al., 2016;Cole et al., 2007;Najjar et al., 2018).Through this biogeochemical pathway, inland waters transport, bury and remove C before it reaches coastal oceans, thereby significantly altering the global C budget.Globally, inland waters emit about 2.2 Pg C yr −1 (0.7-4.2 Pg C yr −1 ) to the atmosphere (Lauerwald et al., 2023a;Raymond et al., 2013), offsetting 75% of the contemporary terrestrial CO 2 sink (Friedlingstein et al., 2019(Friedlingstein et al., , 2022;;Le Quéré et al., 2018).However, large uncertainties still exist in the magnitude and variability of C dynamics along the terrestrial-aquatic continuum.In particular, it is not well understood how natural and anthropogenic disturbances have changed terrestrial C loading to inland waters, the burial of C in inland waters, evasion of CO 2 from inland waters, and C delivery to coastal waters over a century-long time scale, with an estimated 2-σ errors in all of these flux changes at the global scale of between 50% and 100% (Regnier et al., 2022).
Inland water C has three main forms: dissolved organic carbon (DOC), dissolved inorganic carbon (DIC), and particulate organic carbon (POC).Each form responds differently to environmental changes.The long-term changes of C exports in dissolved forms (DOC and DIC) have been found to be primarily regulated by air temperature (Laudon et al., 2012;Pastor et al., 2003), while hydrological conditions may explain their short-term variations (Raymond & Oh, 2007).DIC is also influenced by CO 2 evasion, which is also affected by climate and land-use change (Lauerwald et al., 2015).Land conversions can substantially modify the geomorphologic conditions of the land surface, and thus alter soil erosion and the loss of POC (Galy et al., 2015).Additionally, other direct and indirect anthropogenic factors, such as elevated atmospheric CO 2 concentrations, atmospheric nitrogen deposition, and extensive fertilizer application, can cause notable changes in inland water C fluxes (Findlay, 2005;Houghton, 2010).
Previous estimates of inland water C fluxes were mostly derived from statistical relationships between observations and environmental factors (Dai et al., 2012;M. Li et al., 2017;Ludwig et al., 1996) or from bookkeeping methods (Butman et al., 2016).For example, Dai et al. (2012) estimated the magnitude of global riverine DOC export to be 0.21 Pg C yr −1 based on observations from 118 rivers around the world.However, these methods are highly dependent on the quality and quantity of field measurements, which limits their use in watersheds with scarce data and high spatio-temporal heterogeneity.Therefore, hybrid models integrating empirical, statistical and mechanistic components, such as the Global Nutrient Export from WaterSheds (Global NEWS) model, have been developed to estimate riverine exports of C and nutrients (Seitzinger et al., 2005) and greenhouse gas emissions from inland waters (Hu et al., 2016;Kroeze et al., 2005;Seitzinger et al., 2000).Nevertheless, the wide use of empirical equations and over-simplified representations of mechanisms for land and aquatic systems in these models can undermine their predictive capability as the performance of empirical equations may substantially decrease under changing environmental conditions (Girardin et al., 2008;Leach & Moore, 2019).
In contrast, process-based modeling is more general and hence more suitable for both hindcasting and forecasting under varied conditions.Though several models have been applied in the estimation of inland water C fluxes in regional and global studies (Hastie et al., 2019), these models are usually run at a coarse spatial resolution that is unable to capture key water transport processes, such as channel routing in headwater zones, which have been shown to be hotspots of greenhouse gas emissions (Battin et al., 2008;Butman & Raymond, 2011;Butman et al., 2016;M. Li et al., 2021;McClain et al., 2003;Raymond et al., 2013).Furthermore, it is difficult for these models to represent the effect of individual land-to-water delivery factors and their interactions within the actual ecosystem (Robertson & Saad, 2013).Yet, such processes largely determine how C cycling in aquatic ecosystems is affected by natural and anthropogenic disturbances.The comparison of the current C models (M.Li et al., 2019;Marescaux et al., 2020;Mayorga et al., 2010;Nakayama, 2022;Nakhavali et al., 2020;Saccardi & Winnick, 2021;H. Zhang et al., 2022) reveals their limited ability to differentiate small stream processes in simulating riverine C dynamics.The representation of C transport across the river-lake-reservoir continuum remains incomplete in these models.Moreover, a subset of these models lacks the inclusion of terrestrial C processes, thereby hampering the overall coupling of terrestrial and aquatic C dynamics.
Here, we describe refinements of the Terrestrial/Aquatic Continuum module of the Dynamic Land Ecosystem Model (DLEM-TAC) that improve representations of coupling of terrestrial and aquatic processes of the C, TIAN ET AL. 10.1029/2023GB007776 3 of 24 nitrogen and water cycles, and their integration across multiple water types (rivers, headwater, lakes and reservoirs) (Tian et al., 2015b;Yao et al., 2021).The process-based, high-resolution DLEM-TAC was used to quantify the magnitude and spatio-temporal patterns of global inland water CO 2 evasion, C burial in aquatic sediments, and the riverine exports of POC, DOC, and DIC from land to oceans during 1800-2019.Moreover, factorial simulations were implemented to attribute changes in global inland water C fluxes to different environmental factors, including climate, land-use change, atmospheric CO 2 concentration, nitrogen deposition, and fertilizer/manure application.Our findings demonstrate the important roles of both climate-related factors as well as human activities in alterations in land-ocean-atmosphere C exchanges.

Methods
The DLEM-Terrestrial/Aquatic Continuum (DLEM-TAC) is an integrated modeling framework that encapsulates land and aquatic ecosystem processes (Figure 1).The land component describes terrestrial biophysical characteristics, plant physiological processes, and soil biogeochemistry at a daily time step (Tian et al., 2010); the aquatic component uses C and nutrient loadings as inputs and simulates the biogeochemical dynamics along the land-river-ocean continuum, such as CO 2 evasion, burial, and transport (Yao et al., 2021).Detailed information on the land component can be found in our previous studies (Pan et al., 2021;Tian et al., 2015a).Here, we focus on the improvements in the aquatic component.

The Routing Scheme of Rivers, Lakes, and Reservoirs
In DLEM-TAC, water transport within the grid cells is separated into three subgrid processes: hillslope routing, subnetwork routing, and main channel routing.A scale-adaptive and physically based model named Model of Scale-Adaptive River Transport (H.Li et al., 2013) was incorporated into DLEM.The water from surface runoff is routed across hillslopes first.The water received by the subnetwork channels from hillslope flow and groundwater discharge flows into the main channel.Note that the subnetwork channels within a 30 arc-min grid cell represent the streams from first to fifth orders (Fekete et al., 2001).The main channel receives water from upstream grid cells and local subnetworks and discharges it to the downstream grid cell.All three sub-grid routing processes use kinematic wave methods (Chow, 1964), which require several physical parameters (channel length, bank-full depth, channel slope, and channel roughness) derived from a 15 arc-second resolution hydrological data set (H.-Y.Li et al., 2015).To represent the floodplain process, the model assumed that, the channel width increases by five times when the water level is higher than the bankfull channel depth.In the scale-adaptive water transport scheme, the length of the subnetwork flow changes with the model grid resolution (H.Li et al., 2013), so the total effective length of small streams (subnetworks) within a grid cell increases with the grid size.Therefore, the parameters of the new water transport module only require minor re-calibration when the model is applied at different resolutions.This way, the scale-adaptive river routing scheme within DLEM allows a simultaneous representation of carbon and nitrogen fluxes on scales from small streams to large rivers within a grid cell (Yao et al., 2020(Yao et al., , 2021)).
We should note that, dam operations would substantially affect the flow regime, and the associated carbon fluxes of reservoirs behind them as well as their downstream rivers.Here, we further linked lakes and reservoirs with small streams and large rivers in DLEM, forming a stream-river-lake-reservoir corridor (Wollheim et al., 2008).Specifically, lakes with upstream areas smaller than a grid cell are linked to small streams, while those with upstream areas larger than a grid cell are linked to large rivers.The outflow rates of natural lakes were calculated based on the prescribed residence time obtained from Messager et al. (2016).The residence times of reservoirs that linked subnetwork corridor are obtained from the global data databased (Lehner et al., 2011).Reservoirs that are linked to main channels are considered as major reservoirs.The operation rules of major reservoirs were adopted from existing algorithms (Biemans et al., 2011;Haddeland et al., 2006;Hanasaki et al., 2006), which require a reference run with reservoir operation turned off to provide natural water flows as model input.

Aquatic Carbon Dynamics in DLEM-TAC
The aquatic C module in DLEM-TAC consists of lateral C transport, the decomposition of organic C, the burial of POC, DIC uptake through primary production, and CO 2 evasion (Figure 1a).The dynamics of different C forms within a water body (river, lake, or reservoir) follow the mass balance: TIAN ET AL. 10.1029/2023GB007776 4 of 24 where ∆M POC /∆t, ∆M DOC /∆t, and ∆M DIC /∆t (g C d −1 ), are the net changes in the total mass (M) of POC, DOC, and DIC, respectively; ∆t denotes the time step (day); F a represents the net lateral transport (inflow minus outflow) of C species through the linked inland water corridor (g C d −1 ); R DOC and R POC are the decomposition and catabolism rates of the organic (dissolved and particulate) C species (d −1 ); v s is the settling velocity of POC (m d −1 ); C POC is the concentration of POC (g C m −3 ); A s is the area (m 2 ) of the water body surface, P is the primary production through photosynthesis in the aquatic system (g C d −1 ), and  CO 2 is the CO 2 evasion to the atmosphere (g C d −1 ).TIAN ET AL. 10.1029/2023GB007776 5 of 24 Particulate Inorganic Carbon (PIC) was not considered in our simulation because we assumed that PIC is not reactively involved in the C cycle of terrestrial and aquatic ecosystems.
Carbon species (DOC, POC, or DIC) from the land (surface runoff) enter the hillslope and subsurface flows and further contribute to the subnetwork flow (Figure 1b).Biogeochemical processes within the hillslope flow and subsurface flow were not considered in this work.The advective C fluxes through the subnetwork and main-channel are described as where F a,sub is the C flux (DOC, POC, or DIC) (g C d −1 ) of the subnetwork flow; F a,main is the C flux (g C d −1 ) of the main-channel flow; F h/c is the C flux (g C d −1 ) of the hillslope flow; F g/c is the C flux (g C d −1 ) from the groundwater to the subnetwork; Q sub is the out flow rate of subnetworks (m 3 s −1 ); C sub is the concentration (g C m −3 ) of C flux in the subnetworks; Q up and Q main are the outflow rates of upstream grid cells and the main channel within the grid cell (m 3 s −1 ), respectively; and C up and C main are the associated concentrations (g C m −3 ) of C species associated with Q up and Q main , respectively.
The first-order decomposition and catabolism rate coefficients are given by where K DOC and K POC are the decomposition and catabolism rates (d −1 ) at the reference temperature T s (20°C); Q 10 is the multiplicative factor applied to the respiration rates when the water temperature T w (°C) increases by 10°C relative to T s ; and T w is the water temperature (°C), which is calculated based on an empirical relationship with air temperature (Mohseni et al., 1998(Mohseni et al., , 1999)).
The settling velocity of POC is estimated by a simplified Stokes' law (Thomann & Mueller, 1987): where v s is the settling velocity (m d −1 ); α represents the effect of the particle shape on the settling velocity; ρ s and ρ w are the densities of the particle and water (g • cm −3 ), respectively; and d is the particle diameter (μm).
The CO 2 exchange between water bodies and the atmosphere is explicitly estimated as where  CO 2 is a net CO 2 evasion (g C d −1 );  CO 2 is the gas transfer velocity (m d −1 ), which is obtained from Raymond et al. (2012);  CO 2 is the dissolved CO 2 concentration (g C m −3 ), which is computed from DIC and pH (note that pH is a static map interpolated based on the GLORICH data set without temporal variations) at a given water temperature T w (Hartmann et al., 2019;Yao et al., 2021), and  CO 2 eq is the equilibrium surface water CO 2 concentration (g C m −3 ) with respect to atmospheric pCO 2 .
The gas exchange rate (or refers to piston velocity)  CO 2 (m d −1 ) is estimated as: where  SCCO 2 is the Schmidt Number for CO 2 (unitless), which can be calculated as (Raymond et al., 2012): 10.1029/2023GB007776 6 of 24 where K 600 (m d −1 ) is the gas exchange coefficient.
The K 600 for streams and rivers was estimated by Raymond et al. (2012): In the updated version of DLEM-TAC (Yao et al., 2022), we used different K 600 for lakes and reservoirs, which was adopted from the work of Tan et al. (2015): where U 10 is the wind speed of 10-m above the ground surface (m s −1 ).More detailed information about the representation of the aquatic C dynamics can be found in our previous publication (e.g., Tian et al., 2015a;Yao et al., 2021Yao et al., , 2022)).

Model Forcing
A database at a spatial resolution of 0.5°, consisting of climate variables, land-use change, nitrogen deposition, nitrogen fertilizer, and manure application (Figure S1 in Supporting Information S1), atmospheric CO 2 concentration, and hydrological data was developed to drive the DLEM-TAC model.
The climate variables, which consist of daily precipitation, daily mean temperature, daily maximum temperature, daily minimum temperature, daily wind speed and daily shortwave radiation, were compiled from the CRU-NCEP data set for 1901-2019 (Viovy, 2018).Climate data of each year during 1800-1900 were randomly sampled from the years between 1901 and 1920.The land use data were from HYDE version 3.2 (Klein Goldewijk et al., 2017), and the land cover cohort within a grid cell is composed of four natural vegetation types, one cropland type, and other non-vegetation land-use types (Tian et al., 2018).Detailed information on DLEM land-use input data can be found in M. Liu and Tian (2010).
The annual atmospheric CO 2 concentration from 1800 to 2019 was obtained from the NOAA GLOBALVIEW-CO 2 data set (https://www.esrl.noaa.gov).The nitrogen deposition data set was obtained from the Atmospheric Chemistry and Climate Model Intercomparison Project.The N fertilizer data were obtained from Lu and Tian (2017) and the spatial manure N application data were adopted from B. Zhang et al. (2017) (Figure S2 in Supporting Information S1).
The hydrological input data, covering the flow direction, flow distance, and upstream area, were obtained from the Dominant River Tracing (DRT) data set (Wu et al., 2012).The bank-full width and bank-full depth data sets were obtained from the Hydrological Modeling and Analysis Platform (Getirana et al., 2012).The channel density and channel slopes of small streams and rivers were derived from the National Hydrography Dataset plus v2 data (available at: http://www.horizon-systems.com/NHDPlus/index.php).Surface area, upstream area, volume, depth, and averaged residence time for lakes were obtained from the Hydrolakes data set (Messager et al., 2016).The variables for reservoirs, including built year, height, maximum storage, water surface area, residence time, and upstream area, were obtained from the GRanD database (Lehner et al., 2011).

Simulation Protocol
DLEM-TAC simulation covered the globe primarily following three steps (Figure 2): (a) To obtain the initial climatological or steady state pre-industrial conditions, we conducted the equilibrium simulation for each grid cell by holding all the driving forces constant in the year 1800 including climate status (we used 1901 climate because CRU-NCEP data before 1901 are not available), land-use conversions, atmospheric CO 2 concentration, and nitrogen additions.When the local C, nitrogen, and water pools of all the grid cells reached a steady state, the equilibrium run finished (Thornton & Rosenbloom, 2005).(b) Before moving to the year-to-year normal simulation, we conducted a 30-year second spin-up run by randomly selecting climate forcings within the 1800s (Tian et al., 2012).Then we conduct the natural flow simulation with the dam module temporarily disabled, and all the driving variables changing over time.(c) After the natural flow simulation, we set up a management flow simulation, with the dam module turned on.We evaluated the simulated flow discharge and C fluxes with observations and calibrated the model parameters, including terrestrial C loadings and aquatic processes.Following the model calibration, we conducted the base simulation (S1) by changing all the driving forcing over time to investigate the spatial and temporal patterns of inland water C fluxes.
In order to investigate the responses of inland water C fluxes to multiple environmental changes, six factorial experiments were conducted by fixing each of the driving forcings (Figure 2).For each of the simulations S2 through S6, we held climate, land use, atmospheric CO 2 concentration, nitrogen deposition, fertilizer applications, and manure applications constant at the level of 1800, respectively, while allowing the other variables to track their historical trajectories.The effects of climate, land use, atmospheric CO 2 concentration, nitrogen deposition, fertilizer applications, and manure applications on C loading were calculated by subtracting the factorial experiment from the base simulation (i.e., S1-S2, S1-S3, S1-S4, S1-S5, and S1-S6, respectively).

Model Evaluation
To examine the performance of DLEM-TAC, we compared the simulated terrestrial loading and riverine export with observations for each C species (at 56 land sites and 47 world major rivers.The detailed data used for model calibration and validation are shown in Figure S3, Tables S1 and S2 in Supporting Information S1).
The determination coefficient (R 2 ) and Nash-Sutcliffe efficiencies (NSE, Nash and Sutcliffe (1970)) were applied to assess the performance of DLEM-TAC.The NSE ranges between  −∞ and 1.An NSE close to 1 means a good match of the simulated to the observed data; NSE = 0 means that the simulation is as skillful as the average of the observed data (M.Li et al., 2019).Simulated terrestrial loading of DOC, POC, and DIC agreed well with observations (log transformed, R 2 > 0.7; Figure 3a-3c).Aquatic module simulations of the observed annual mean discharges and the C export to the coast were consistent with those recorded by the GEMS-GLORI database (Meybeck & Ragu, 2012), with R 2 values of the log of discharge and export of DOC, POC, and DIC being 0.8, 0.8, 0.6, and 0.7, respectively (Figures 3d-3g), which can be considered as satisfactory (Moriasi et al., 2015).Furthermore, the simulated monthly riverine C exports of the major rivers by  We evaluated the simulated performance of seasonal variability on DOC exports at 10 sites, with 3 rivers located in the tropical region, 5 rivers in the temperate region, and 2 rivers in the boreal region.Comparisons between simulated and database observations had R 2 higher than 0.5 and NSEs higher than zero at most sites.We did not evaluate monthly POC exports due to lack of data availability.Additionally, simulated total organic carbon exports of 5 temperate rivers and 3 boreal rivers were comparable to observed values with R 2 values in most cases higher than 0.4 and NSEs higher than zero.Finally, the DIC exports simulated by DLEM-TAC in 3 tropical rivers and 9 temperate rivers were representative of observations with R 2 values higher than 0.5 and higher than zero at most sites.

Uncertainty Analysis
We evaluated two major sources of uncertainty in estimating global inland water C fluxes, namely terrestrial C loading and water surface area.To evaluate the uncertainty induced by the former, we carried out a literature survey to collect the observed C loading at the site level.We conducted linear regressions to evaluate the model estimated C loadings against the observations (Figure 3).The 95% confidence band of the linear regressions of DOC, POC, and DIC loadings vary about ±20%, ±30%, and ±40%, respectively, with respect to the 1:1 line.We then conducted two model simulations from 1800 to 2019 with the parameters of terrestrial loadings for DOC, POC, and DIC varying ±20%, ±30%, and ±40%, respectively, to represent their minimal and maximal range.For water surface area, we did not investigate the uncertainty originating from the areas of lake and reservoirs and large rivers, because they can be measured using remotely sensed products, leaving the area of headwater streams as the most uncertain one.We therefore implemented an uncertainty analysis for the river shape parameter (r in Text S3 in Supporting Information S1) to represent the global river surface areas varying from 0.89, 0.81, 0.73, 0.65, and 0.56 (10 6 km 2 ), which aligns well with previous estimates (Allen & Pavelsky, 2018;Bastviken et al., 2011;Raymond et al., 2013).We averaged the simulated mean values for two sources of uncertainties and the best-performing simulation, is presented as the final results in this study.The combined standard deviation was calculated as the square root of the sum of the squared individual standard deviations corresponding to each uncertainty source.

The Contemporary Carbon Budget of Global Inland Water Systems
The simulated total terrestrial C loading during 2010-2019 was estimated to be 2281 ± 640 Tg C yr −1 (mean ± 1 standard deviation of the annual average, n = 20, Figure 4), with DIC dominating, followed by DOC and POC.More than one third (38%) of the C loading (869 ± 151 Tg C yr −1 ) was exported to the coastal area, of which DIC accounted for half with DOC and POC occupying a similar share.CO 2 evasion was 1,113 ± 251 Tg C yr −1 , mainly contributed from rivers (777 ± 179 Tg C yr −1 , of which 640 ± 150 Tg C yr −1 from headwater streams and 137 ± 34 Tg C yr −1 from high-order rivers).The C buried in inland waters was 232 ± 41 Tg C yr −1 , much smaller than C export and CO 2 evasion.

Long-Term Trend in Inland Water Carbon Fluxes Since the 1800s
In the 1800s, the terrestrial C loading was estimated to be 1,820 ± 507 Tg C yr −1 , and more than half of the loading was composed of DIC (1,207 ± 82 Tg C yr −1 ).An estimated 42% of the C loading was exported from rivers to the ocean (765 ± 132 Tg C yr −1 ), of which DIC accounted for the largest proportion (404 ± 82 Tg C yr −1 ), followed by DOC (246 ± 31 Tg C yr −1 ) and POC (114 ± 20 Tg C yr −1 ).Nearly half of the C loading leaving inland waters was through CO 2 evasion (841 ± 190 Tg C yr −1 ) (Table 1).
Model simulations of terrestrial C loading increased gradually from the 1800s (1,820 ± 507 Tg C yr −1 ) to the 2010s (2281 ± 640 Tg C yr −1 ), at an average rate of 2 Tg C yr −1 .Compared with the 1800s, the total terrestrial C loading in the 2010s increased by 25%, with the loadings of DIC, DOC, and POC increasing by 22%, 12%, and 62%, respectively (Table 1).Enhanced terrestrial C loading resulted in an increase of 31% in the inland water CO 2 evasion compared with the pre-industrial level, predominantly from headwater streams (Table 1).In addition, due to climate change and human activities, CO 2 evasion from inland waters increased significantly after the 1950s, at a rapid rate of 2 Tg C yr −1 (Figures 4 and 5).As a result, the fractional increase in riverine C export (14%) from the 1800s to the 2010s was smaller than the corresponding increase in C loading (Figure 5).

Spatio-Temporal Patterns of Inland Water Carbon Fluxes
At the global scale, tropical regions and the middle and high latitudes of the Northern Hemisphere were the hot spots of the terrestrial C loading.From the pre-industrial period to the recent decade, the terrestrial DOC loading increased in regions of the tropics and high latitudes, such as the Amazon basin, southeastern Asia, and Europe (Figures 6a-6d).Similarly, POC loading increased in the tropics and the regions around 30-45°N (Figures 6e-6h).At the same time, a small part of the Amazon basin was characterized by a decrease in POC loading.In contrast to organic C loading, the largest increase in DIC loading occurred in regions between 30°N and 60°N, such as eastern North America and Europe (Figures 6i-6l).
Consistent with the pattern of C loading, the hot spots of riverine CO 2 evasion were primarily located in the middle and high latitudes of the Northern Hemisphere and in the tropics (Figures 7a-7d).Compared with the pre-industrial period, the contemporary CO 2 evasion from streams and rivers increased notably around the Equator and in the temperate regions of the Northern Hemisphere (30°N-60°N).Due to the abundance of lakes, the regions around 50°N contributed most of the lacustrine CO 2 evasion.Since the magnitude of riverine C loading was much larger than that of lakes and reservoirs, the total CO 2 evasion from inland waters was dominated by streams and rivers (Figures 7i-7l).
We analyzed the global inland water C fluxes in the 10 main regions defined by the REgional Carbon Cycle Assessment and Processes Phase 2 project (Ciais et al., 2022).In the recent decade, the largest riverine C exports were from South America, Southeast Asia, North America, and Africa, which accounted for 27%, 17%, 11%, and 11% of the global riverine C exports, respectively (Figure 8).Large increases in C exports were found in North America, Southeast Asia, East Asia, and Russia, accounting for 21%, 15%, 17%, and 18% of the increase in global riverine C exports, respectively.Regions of North America, Europe, Russia, South America and East Asia contributed to most of the CO 2 evasion from global inland waters (23%, 18%, 17%, 12%, and 8%, respectively), and experienced large increases in CO 2 evasion since the pre-industrial period (20%, 27%, 27%, 9%, and 5% of the increase in global CO 2 evasion, respectively) (Figure 8).

Environmental Controls Over Inland Water Carbon Fluxes
We quantified the contributions of key environmental drivers to the changes in inland water C fluxes through multiple factorial simulation experiments as described in Figure 3.The results showed that, from the pre-industrial period to the recent decade, increases in atmospheric CO 2 concentration caused increases of 17, 15, and 13 Tg C yr −1 in the exports of DOC, POC, and DIC, respectively (Figure S6 in Supporting Information S1).Nitrogen additions from fertilizer and manure application and nitrogen deposition were associated with increased exports of DOC, POC, and DIC export by 18, 23, and 25 Tg C yr −1 , respectively, over the same time period (Figure S6 in Supporting Information S1).Climate changes primarily influences the exports of DIC (22 Tg C yr −1 ), followed by the exports of DOC (12 Tg C yr −1 ) and POC (7 Tg C yr −1 ) (Figure S6 in Supporting Information S1).Landuse change resulted in an increase in POC exports and DIC exports by 20 and 27 Tg C yr −1 in the 2010s, but reduced DOC exports by 2 Tg C yr −1 (Figure S6 in Supporting Information S1).Changes in climate and land-use were the dominant factors that increased the CO 2 evasion, by about 88 and 108 Tg C yr −1 , respectively, from the 1800s to the 2010s, while other factors (elevated atmospheric CO 2 , nitrogen deposition, and fertilizer and manure application) contributed 81 Tg C yr −1 (Figure 9).Terrestrial C loadings have similar response patterns to environmental change as riverine C exports (Figure 9).Changes in Climate and land use were the primary factors contributing to changes in terrestrial DOC and DIC loadings during the 1800s-2010s, respectively (Figure S7 in Supporting Information S1).The impact of climate change on terrestrial POC loadings was predominant before the 1950s, while anthropogenic disturbances became more influential than climate change in increasing POC loading after that period (Figure S7 in Supporting Information S1).
The spatial pattern of the dominant factors affecting terrestrial C loading and inland water CO 2 evasion also changed from the early 19th century (1800-1819) to the early 21st century (2000-2019) (Figure 10).In the early 19th century, climate variables were largely responsible for the change in C loading and CO 2 evasion.In comparison, the changes in small parts of Europe, South Asia, coastal regions in East Asia, and southeastern regions in North America were dominated by land use (Figures 10a and 10c).As compared to the early 19th century, the area dominated by land use and nitrogen addition significantly expanded in the early 21st century.Increasing atmospheric CO 2 concentration has become an important factor dominating changes in C loading and CO 2 evasion in Africa (Figures 10b and 10d).

Comparison to Previous Studies
The C fluxes estimated by DLEM-TAC fall within the ranges given by previous studies, showing the consistency of this process-based model with other approaches, which are mostly data-driven (Table 2).DLEM-TAC DOC export (averaging 262 ± 33 Tg C yr −1 in the 2010s) is in good agreement with six previous studies, which ranged from 171 to 246 Tg C yr −1 (Table 2).The estimation of global DOC loading of about 280 Tg C yr −1 reported by Dai et al. (2012) and Nakhavali et al. ( 2020) is consistent with our estimate (averaging about 431 Tg C yr −1 in the 2010s) after the C loading from organic soil C is accounted for (around 170 Tg C yr −1 , as reported by Nakhavali et al. (2020)).The significant increasing trend of global DOC export from 1950 to 2019 simulated by DLEM-TAC is consistent with the finding of the Global NEWS model (Seitzinger et al., 2005) but contrasts with the report by M. Li et al. (2019), who claimed a decreasing trend in global DOC export.Nevertheless, DLEM-TAC simulations agree well with the results by M. Li et al. (2019) that DOC export showed a statistically significant increase in tropical regions (see also Lauerwald et al. (2020)).In contrast with the results by M. Li et al. (2019), DLEM-TAC estimated slight increases in arctic DOC export, which is supported by Bowring et al. (2020).A possible reason for such a difference is that the representations of impacts of freeze-thaw cycles are different, which would cause an increase in C loading in recent studies (Rawlins et al., 2021).To address this issue, more observations are needed to constrain the models.Similarly, DLEM-TAC POC export (154 ± 25 Tg C yr −1 in the 2010s) and DIC export (453 ± 94 Tg C yr −1 in the 2010s) are in line with the mean of previous studies (196 Tg C yr −1 for POC and 413 Tg C yr −1 for DIC) (Table 2).The significant increase in DIC export from 1960 to 2019 revealed by DLEM-TAC is supported by a previous data synthesis (M.Li et al., 2017).We should note that, our estimated burial C in rivers is slightly higher than that in lakes and reservoirs, due to the high deposition rate within river floodplains.2).However, DLEM-TAC estimation of total C export is still larger than the one given by Resplandy et al. (2018) because their estimate only accounts for part of the river C export that is outgassed back to the ocean, which does not include the buried fraction (Table 2).The difference may be due to the C loss in estuaries (burial and outgassing), which occurs between the river outlet and the adjacent coastal ocean.Although the total terrestrial C loading given by DLEM-TAC (2281 ± 640 Tg C yr −1 in the 2010s) is much larger than the estimation provided by Cole et al. (2007) and is lower than those by Battin et al. (2023), several studies (Battin et al., 2009;Regnier et al., 2013) support our estimation (Table 2).The performances of the different methods in the estimation of riverine C exports and loadings need to be further evaluated in the future.
Our estimation of riverine CO 2 evasion was in the range of previous studies (Table 3).The riverine CO 2 evasion estimated by DLEM-TAC was around twice as much as the estimate by Cole et al. (2007) (Table 3), probably because their study ignored evasion from small rivers, which are now considered to be an important CO 2 source (M.Li et al., 2021).Our riverine CO 2 evasion is much lower than that given by Raymond et al. (2013) and S. Liu et al. (2022), with the estimated inland water CO 2 emission reached 2.5 Pg C yr −1 .Based on their calculation, the global land C loading can be larger than 3.4 Pg C yr −1 (inland water CO 2 emission plus global riverine export), which is 110% of the estimated C sink by the global terrestrial ecosystem (3.1 Pg C yr −1 ) (Friedlingstein et al., 2022).However, the average fraction of leaching/net ecosystem exchange is about 20% for most of the observation sites in European sites (Kindler et al., 2011), which implies a potentially significant overestimation of inland water CO 2 emissions in the previous studies.The magnitude of lacustrine CO 2 evasion estimated by DLEM-TAC was higher than that of Cole et al. (2007), but lower than that of Raymond et al. (2013).The differences mainly arise from the tropical zone, since most of the observations were collected in the temperate and boreal regions.Hastie et al. (2018) reported a total CO 2 evasion of 272 Tg C yr −1 from lakes at high latitudes, TIAN ET AL.

10.1029/2023GB007776
14 of 24 which is almost equal to our estimate for the CO 2 evasion from global lakes.A reason for such a difference is that the lake area they used was much larger than that used in our study (Messager et al., 2016).The CO 2 evasion from reservoirs presented here is similar to the data synthesis estimate from Deemer et al. (2016), but far lower than that of St. Louis et al. (2000) because of their overestimation of global reservoir surface area (Johnson et al., 2021).
The global CO 2 evasion from inland waters in this study is close to the empirical estimations (Deng et al., 2022).Our simulated CO 2 evasion in the 2010s for most 106 river basins is comparable to their empirical estimates.Nevertheless, our estimates of CO 2 evasion from the Amazon and Congo Basins are lower than those of Deng et al. (2022) and Byrne et al. (2023) (Figure S5 in Supporting Information S1).This is because we do not consider the CO 2 evasion from wetlands and floodplains in our estimation, which are considered as the dominant source of CO 2 evasion from surface water in these two basins (Borges et al., 2015;Richey et al., 2002).Furthermore, inland water CO 2 evasions as well as the riverine C exports from Europe we estimated are much higher than those presented by H. Zhang et al. (2022).It is worth noting, however, that our estimations regarding European riverine C exports aligned closely with the magnitude documented in another global study (M.Li et al., 2017).We showed much lower CO 2 evasions in comparison to those by S. Liu et al. (2022).Our study may underestimate CO 2 evasions from rivers in mountainous regions such as the Rocky Mountains, the Andes, and the Himalayas, considering the high gas transfer velocity documented in S. Liu et al. (2022)'s study.These differences also suggest critical needs to better understand processes controlling C dynamics along the land-aquatic interface as well as better represent these processes in our model.To enhance the credibility of projections concerning inland water C fluxes on a global scale, it is imperative that future projection endeavors to prioritize the harmonization of results obtained from both global and regional studies, particularly within targeted regions.From the pre-industrial period to the recent decade, our increased rates of terrestrial C loading (25%), riverine C export (14%), and CO 2 evasion (32%) are in good agreement with the results reported by Regnier et al. (2022) (26% for terrestrial C loading, 12% for riverine C export, and 28% for CO 2 evasion).

Natural and Anthropogenic Drivers of the Changes in Inland Water Carbon Fluxes
The inter-decadal variations in model-simulated inland water C fluxes are mainly controlled by climate-related variability (Figure 9).Changes in precipitation and air temperature can directly affect ecosystem production and thus the C input to soil (Cleveland et al., 2011;Z. Zhang et al., 2016).Simultaneously, the shifts in precipitation patterns can largely alter runoff and further impact soil leaching and C transport in inland waters (Lal, 2005).
Climate variation together with elevated CO 2 significantly promoted vegetation growth and thus increased C loss from land to inland waters (P.Li et al., 2017).However, our modeling results showed that the contemporary DIC export is comparable to its pre-industrial level, which is supported by evidence derived from observations (Raymond & Hamilton, 2018).This is because stream water temperature is strongly correlated with air temperature (Mohseni et al., 1998(Mohseni et al., , 1999)), the increase in which would increase organic matter decomposition, accelerate gas exchange rates, and substantially increase CO 2 evasions in waters.
The estimated POC export increased substantially from 1900 to 2019, primarily due to land-use change.Most of the increase in POC export occurs in regions where a large area of cropland is converted from forest (Figures S1 and S8 in Supporting Information S1).This is because the erodibility of croplands is much higher than that of forests and grasslands (Williams & Berndt, 1977) leading to rapid loss of POC accumulated previously.The considerable increase in POC export in the arctic region can be mainly attributed to the high sensitivity of the ecosystems to climate change (Hilton et al., 2015) (Figure 9 and Figure S8 in Supporting Information S1).
Overall, nitrogen inputs including nitrogen deposition, nitrogen fertilizer and manure nitrogen applications are positively correlated with riverine C fluxes.We observed a notable increase in nitrogen effect on riverine C fluxes, specifically in recent years, in response to the consistent increasing rates of nitrogen deposition and nitrogen applications.In theory, the nitrogen inputs can alleviate the vegetation nitrogen limitation of global terrestrial ecosystems (Vitousek & Howarth, 1991), which contributes to the increased primary production, soil C pools and the associated C loadings from land to inland waters.

Significance of Inland Water Carbon Fluxes to the Global Carbon Budget
Increased evidence indicates that aquatic C fluxes need to be accounted for in revisions of the overall C balance of the terrestrial ecosystem (Cole et al., 2007;Friedlingstein et al., 2019).Constrained by fossil fuel emissions, atmospheric CO 2 growth rate, and C fluxes within fresh and salt water ecosystems, the global net CO 2 sink by terrestrial vegetation was estimated to be 2.2 Pg C yr −1 in the recent decade (Regnier et al., 2022).Our study indicated that inland water recycles and exports nearly half of the net land C sink into the atmosphere and ocean, highlighting the important role of inland waters in the global land C balance, an amount that should be taken into account in future C budgets.In addition, riverine C export is likely to provide substantial substrate to support microbial activity in the ocean and is an important term in the coastal C budget (Barrón & Duarte, 2015).
Based on our reconstruction of historical C fluxes between atmosphere, land, and water, lateral C fluxes in the pre-industrial period were found to be considerable.The anthropogenic activity has likely perturbed lateral C fluxes in the land-ocean aquatic continuum (LOAC), but this perturbation has been difficult to quantify (Regnier et al., 2022).Our work has quantified the long-term anthropogenic perturbation of lateral C fluxes and found a notable increase (25%) of C transfer from land to inland waters since the pre-industrial period, which intensified inland water CO 2 evasion, and to a lesser extent, export to the ocean.The global C budget of the Global Carbon Project assumed that the C fluxes in the LOAC have not changed since pre-industrial time (Friedlingstein et al., 2021).In contrast, we have highlighted the influence of multiple anthropogenic perturbations on lateral C   (Regnier et al., 2013;Yvon-Durocher et al., 2012).

Uncertainty and Future Research
A major source of uncertainty in the estimation of inland water C fluxes is associated with the input data.For instance, the absence of seasonal variations in the spatial distribution of nitrogen inputs, including atmospheric nitrogen deposition, fertilizer and manure applications, may bias the estimation, as the timing of N inputs has substantial effects on land C and nitrogen cycles (Scharf et al., 2002).Similarly, the atmospheric CO 2 concentration is fixed throughout the year.In fact, atmospheric CO 2 is unevenly distributed across the globe and seasons (Keeling et al., 1989).In this study, we used a static global pH map for calculating CO 2 concentration in aquatic systems, which may introduce uncertainties by ignoring the temporal changes of pH level over time (Stets et al., 2014).In this study, we used five simulated river water surface areas as the uncertainty range, where values (0.73-0.81) close to those of Allen and Pavelsky (2018) can be considered the most robust one.However, the estimates of river water surface area still suffer great uncertainty due to the limit of spatial resolution of remote sensing products.The alteration of surface water area resulting from climate variability and extensive human activities holds considerable implications for inland water CO 2 evasions.For instance, our estimations reveal a contrasting trend in the upward trajectory of inland water CO 2 evasions in East Asia over recent decades compared to a previous study (Ran et al., 2021).This disparity can be attributed to the limitations of using static river data sets, needing to account for the declining riverine area in China during this period.Therefore, a global dynamic inland water network should be a critical task for future research.Additionally, other driving forces, including long-term climate variables and land-use change, are also limited by spatial resolution and data availability.
Uncertainties can also arise from model parameterization and the model structure.Compared to the major component of current Earth system models that can be used to simulate C flux in inland waters, DLEM-TAC is the most   10.1029/2023GB007776 18 of 24 comprehensive model that fully couples aquatic and land processes by incorporating all the inland water types (river, lake and reservoir) with explicit representation of small stream and large rivers.DLEM-TAC has the ability in simulating three carbon species (DOC, DIC and POC) as well as inland water CO 2 evasion (Table S3 in Supporting Information S1).However, the parameters relevant to the aquatic biogeochemical reaction rates in inland waters, such as the organic matter decomposition rate, particle deposition rate, and CO 2 air-water exchange velocity, are largely fixed for all river systems in our simulation.Although our parameter values fall into the ranges of previous studies, the parameters likely vary across river systems.For instance, S. Liu et al. (2022) demonstrated the substantial seasonal variability of gas transfer velocity across different climate zones.Moreover, their study also highlighted the intricate and non-linear association between gas transfer velocity and basin conditions, such as terrain and gas exchange patterns.Employing a very coarse empirical model for the estimation of CO 2 gas transfer velocity may introduce significant errors, particularly in streams characterized by steep terrain.More monitoring sites located within headwater zones are needed to better constrain these parameters.In addition, the current model still lacks a full representation of hydrodynamics and C-associated biogeochemistry although much progress has been made.For instance, a better representation of carbon dynamics in hillslope and groundwater would significantly influence the ratio of inorganic carbon over organic carbon, especially in riparian and hyporheic zones, which need more observations to support the algorithm design and parametrization.Also, we do not consider the vertical stratification of lakes and reservoirs, and the associated carbon reaction and CH 4 flux, which requires systematic observations to quantify the carbon fluxes of multiple interfaces.

Conclusion
In this study, the DLEM-TAC model, which integrates both land and aquatic C processes, was applied to quantify the global inland water C budget from 1800 to 2019.Based on our DLEM-TAC estimates over recent decades, we estimate that around 2.3 Pg C yr −1 entered inland water ecosystems in the 2010s, and a major proportion of C loading was eventually emitted by inland waters to the atmosphere as CO 2 (1.1 Pg C yr −1 ) and exported to the oceans by rivers (0.9 Pg C yr −1 ).Under anthropogenic disturbances, a large increase in terrestrial C loading produced a corresponding increase in CO 2 evasion from inland waters from the 1800s to 2010s.A sustained increase in the exports of DIC, DOC, and POC as well as CO 2 evasion can be attributed to climate change, land-use change, and N application.Although uncertainties exist, our results suggest that inland water C dynamics play a critical role in the global C budget.Our study indicated that inland water recycles and exports nearly half of the net land C sink into the atmosphere and oceans, highlighting the important role of inland waters in the global land C balance, an amount that should be taken into account in future C budget assessment.More observations in headwater zones, arctic regions, and improved process-based modeling tools are needed to better constrain the estimates of inland water C fluxes.

Figure 1 .
Figure 1.Sources and fates of major carbon species in inland water systems (river, lake, and reservoir) represented in the DLEM-TAC modeling framework (a) and the representation of small rivers within the concept model of the scale-adaptive water transport module (b).

Figure 2 .
Figure 2. The flowchart representing the simulation protocol for the lateral carbon fluxes.
performed well in representing the seasonal variability (FigureS4in Supporting Information S1).

Figure 3 .
Figure 3. (a-c) Comparisons of simulated carbon loading, (d) river discharge, and (e-g) carbon export with observations.The original units of the carbon loadings, discharge, and carbon exports were g C m −2 yr −1 , 10 9 m 3 yr −1 , and Gg C yr −1 , respectively.All data are plotted in log10 scale.In the subplots (a-c)c, the error bars of the observations represent the standard deviation; the error bars of the simulation represent the standard deviation of the simulation from 1981 to 2015.The sources of observed data used to validate carbon loading and carbon exports are provided in TablesS1 and S2in Supporting Information S1.The red bands in subplot (a-c) represents the 95% confidence bands.

Figure 4 .
Figure 4.The global inland water carbon budget in the 2010s.Note that carbon burial in rivers includes processes within bankfull channels and floodplains.

Figure 5 .
Figure 5.The long-term trajectories of decadal averages of (a) terrestrial carbon loading, (b) riverine exports, and (c) inland water CO 2 evasion from the 1800s to the 2010s.Note that the vertical scales of the subplots are different.The shaded areas represent decadal variations with ±1 standard deviation.

Figure 6 .
Figure 6.Spatial patterns of terrestrial dissolved organic carbon (top row), particulate organic carbon (middle row), and dissolved inorganic carbon (bottom row) loadings in the preindustrial period (1800s), and the contemporary period (2010s), and the change from the 1800s to the 2010s (third column).The fourth column shows the latitudinal averages of the loadings in the 1800s and 2010s (the shading area indicate uncertainty ranges with ±1 standard deviation).

Figure 7 .
Figure 7. Spatial patterns and latitudinal distribution of riverine CO 2 (top), lacustrine CO 2 (middle), and all inland water CO 2 (bottom) evasion in the preindustrial period (1800s) and the contemporary period (2010s), and the change from the 1800s to the 2010s (third column).The fourth column shows the latitudinal averages of the loadings in the 1800s and 2010s (the shading area indicate uncertainty ranges with ±1 standard deviation).

Figure 8 .
Figure8.The temporal patterns of riverine carbon (dissolved organic carbon, particulate organic carbon, and dissolved inorganic carbon) exports and CO 2 evasion from inland water systems in 10 regions (NA: North America; EU: Europe; WAS: West Asia; RUS: Russia; EAS: East Asia; SEAS: Southeast Asia; OCE: Oceania; SAS: East Asia; AF: Africa; and SA: South America.).Note that the shading areas represent decadal variations with ±1 standard deviation.

Figure 9 .
Figure 9. Decomposition of the factors influencing the long-term changes of (a) terrestrial carbon loading, (b) riverine carbon export, and (c) inland water CO 2 evasion from the 1800s to the 2010s.
should be considered for a robust quantification of the modern global C cycle and climate change mitigation

Figure 10 .
Figure 10.The dominant impact factors on terrestrial carbon loading (a, b) and inland water CO 2 evasion (c, d) in the early 19th century (1800-1819) and the early 21st century (2000-2019).The factor that caused the maximum changes is shown in each grid.Nitrogen addition includes nitrogen deposition and the application of nitrogen fertilizer and manure nitrogen.
The unit of the inland water CO 2 evasion is Tg C •yr −1 .a This study only analyzed evasion from large rivers.b The value includes CO 2 evasion from natural lakes and the lakes with dams.

Table 1
Global Carbon Fluxes Along the Land-Aquatic Interface (Tg C yr −1 , Uncertainty Ranges With Mean ± std)

Table 2
Comparison of Various Estimates on Global Carbon Exports From Rivers to Oceans TIAN ET AL.

Table 3
Comparison of Various Estimates of Global CO 2 Evasion From Inland Water Systems