Organic Carbon (DOC) leaching to the European river

Leaching of dissolved organic carbon (DOC) from soils to the river network is an important 15 component of the land carbon (C) budget. At regional to global scales, its significance has been estimated 16 through simple mass budgets, often using multi-year averages of observed fluvial DOC fluxes as proxy of DOC 17 leaching due to the limited availability of observations of the leaching flux itself. This procedure leads to a 18 systematic underestimation of the leaching flux because of the reactivity of DOC during fluvial transport. 19 Moreover, this procedure does not allow to reveal spatio-temporal variability in DOC leaching from soils, which 20 is needed to better understand the drivers of DOC leaching and its impact on the local soil C budget. In this 21 study, we use the land surface model ORCHILEAK to simulate the terrestrial C budget including leaching of 22 DOC from the soil and its subsequent reactive transport through the river network of Europe. The model 23 performance is not only evaluated against the sparse observations of DOC leaching, but also against the more 24 abundant observations of fluxes and reactivity of DOC in rivers, providing further evidence that our simulated 25 DOC leaching fluxes are realistic. The model is then used to simulate the spatio-temporal patterns of DOC 26 leaching across Europe over the period 1972 to 2012, quantifying both the environmental drivers of these 27 patterns as well as the impact of DOC leaching on the land C budget. Over the simulation period, we find that on 28 average 14.3 TgC yr -1 of DOC is leached from land to European rivers, which is only about 0.6% of the 29 terrestrial net primary production, a fraction about one order of magnitude lower than reported for tropical river 30 networks. Of the DOC leaching, on average 12.3 TgC yr -1 is exported to the coast via the river network, the rest 31 being respired in transit. DOC leaching presents a large seasonal variability, with a maximum occurring in 32 winter and a minimum in summer, except for the Northern most part of Europe where the maximum occurs in 33 spring due to the snow melt. DOC leaching rate is generally lower in warm and dry regions, and higher in cold 34 and wet regions of Europe. Furthermore, runoff, and the ratio between runoff from shallower flow paths vs. deep 35 drainage and groundwater flow, is the main driver of the spatial variation of DOC leaching. Temperature, as a 36 major control of DOC production and decomposition rates in the soils, plays only a secondary role. 37 https://doi.org/10.5194/esd-2021-44 Preprint. Discussion started: 22 June 2021 c © Author(s) 2021. CC BY 4.0 License.

DOC from the soil and its subsequent reactive transport through the river network of Europe. The model 23 performance is not only evaluated against the sparse observations of DOC leaching, but also against the more 24 abundant observations of fluxes and reactivity of DOC in rivers, providing further evidence that our simulated 25 DOC leaching fluxes are realistic. The model is then used to simulate the spatio-temporal patterns of DOC 26 leaching across Europe over the period 1972 to 2012, quantifying both the environmental drivers of these 27 patterns as well as the impact of DOC leaching on the land C budget. Over the simulation period, we find that on 28 average 14.3 TgC yr -1 of DOC is leached from land to European rivers, which is only about 0.6% of the 29

Introduction 38
Land ecosystems are an important carbon (C) sink which absorbs about one fourth of anthropogenic CO 2 39 emissions and stores it in increasing biomass and soil carbon pools . This terrestrial 40 sink mitigates the growth rate of atmospheric CO 2 and thus plays an important role in regulating climate change 41 (Ciais et al., 2013). However, the efficiency of that sink is partly alleviated by the permanent, lateral leaching of 42 C from soils, through the river network down to the ocean (Regnier et al., 2013). An accurate understanding of 43 lateral carbon fluxes through the river network is thus necessary to better understand the global C cycle and to 44 inform policies of climate change mitigation (Le Quéré et al., 2018). 45 The identification of riverine C transfers as a key component of the continental C budget constituted an 46 important paradigm shift in our understanding of the global C cycle (Cole et al., 2007). More recently, riverine C 47 cycling was also shown to be affected by anthropogenic perturbation and thus to be an element of the 48 anthropogenic CO 2 budget (Regnier et al., 2013;Le Quéré et al. 2015). Anthropogenic perturbations of riverine 49 C fluxes are manifold and comprise direct impacts through changing C and nutrient inputs following land-use 50 change and agricultural activities, wastewater discharge, and hydraulic management (e.g. Tian et al., 2015;51 Lauerwald et al., 2020;Hastie et al., 2021;Maavara et al., 2017). There are also indirect impacts following 52 climate change and changes in atmospheric composition. Together, these perturbations have accelerated the 53 turnover of C along the terrestrial-inland water continuum. The terrestrial C sink, which is classically estimated 54 without taking into account the C exports through the river network, is thus generally overestimated (Regnier et 55 al., 2013;Lauerwald et al., 2020). 56 The integration of riverine C transfers into the terrestrial C budget requires the quantification of the amount of C 57 lost from soils to the river network. Due to the scarcity of observational data, this flux is not easy to determine 58 based on empirical methods. At global scale, this flux was estimated through budget closure based on estimates 59 of riverine C exports to the coast and estimates of C losses to the atmosphere and aquatic sediments during 60 transport. The existing global estimates of these soil C exports to the river network, as synthesized by Drake et 61 al. (2018), range from 1.1 to 5.1 PgC yr -1a huge uncertainty range reflecting the limitations of empirical 62 estimation approaches and the paucity of underlying data. Over the past decade, a new generation of land surface 63 models (LSMs) have been developed, which represent the export of C from soils to the river network, and in 64 some cases even the transport and cycling of these terrestrial C loads along the river network down to the coast 65 (Smith et al. 2010;Kicklighter et al. 2013;Tian et al., 2015;Lauerwald et al., 2017;Nakhavali et al., 2018). 66 With the exception of the study by Tian et al. 2015, all these studies focus on the lateral export of dissolved 67 organic C (DOC) which is a product of the incomplete decomposition of plant litter and soil organic carbon 68 (SOC). These mechanistically based models allow to predict the leaching of DOC in unmonitored regions and to 69 assess the spatial and temporal variability which, to date, can only be poorly resolved by empirical methods. 70 Moreover, these approaches link the C exports from soils to the river network to the terrestrial C cycle, and thus 71 allow to directly assess the role of these C exports in the terrestrial C budget, its perturbation through land use, 72 land use change and changes in climate and atmospheric chemistry, and its impact on the terrestrial sink for 73 anthropogenic CO 2 emissions. 74 https://doi.org/10.5194/esd-2021-44 Preprint. Discussion started: 22 June 2021 c Author(s) 2021. CC BY 4.0 License. is represented by four different litter and three different soil organic carbon (SOC) pools with different turnover 136 rates. The four litter pools correspond to metabolic aboveground and belowground litter, structural aboveground 137 and belowground litter (Fig. 1). The SOC is subdivided into active, slow and passive pools. In the CENTURY 138 scheme, C from the decomposed structural litter enters the active and the slow pools with the fraction allocated 139 to each pools depending on lignin content of the litter, while the entire metabolic litter pool and the remaining 140 part of structural litter is allocated to the active SOC pool. The SOC pools then feed into each other with the 141 main C flux going from active to slow and passive to represent microbial decomposition of detrital organic 142 matter, and a small return flux of slow and passive C back to the active pool to represent implicitly the C supply 143 in the form of dead microbial biomass. 144 145 146 where i stands for the ith layer, z for the depth along the discretized soil profile, and D stands for the molecular 164 diffusion coefficient of DOC. 165 The advective export of DOC to the river network is proportional to the top (first five layers, 4.5 cm) and bottom 166 (11th layer) DOC concentrations, corresponding to water loss fluxes associated to runoff (for near surface) and 167 drainage (for deep). We also apply a Fickian-type transport to represent the effect of bioturbation on DOC 168 profiles, as well as on the solid phase SOC and litter profiles. Therefore, in contrast to advection, which only 169 affects DOC entrained with water losses, bioturbation transports both SOC and DOC across the different layers. 170 In addition, DOC is diffusing along the concentration gradient within the soil solution. 171 The right hand-side of Fig. 1 summarizes the set of production/decomposition processes that occur in each layer. 172 During litter decomposition, a fraction of the C is directly emitted back to the atmosphere as CO 2 while the 173 remainder feeds the active and slow SOC pools: 174 175 (4) where k L is the kinetic constant for the litter decomposition (dependent on soil moisture and temperature 176 (Camino et al., 2018)) and ω L the fraction of litter that is channelled into DOC production (as opposed to 177 particulate soil carbon SOC). This approach of relating DOC production directly to the decomposing litter is 178 inspired by Nakhavali et al., 2018 (following from the ECOSSE model (Smith et al., 2010)) and is a major 179 modification compared to the previous version of soil DOC and POC cycling from Camino et al., 2018. In 180 equations (4) and (5), the partitioning between SOC production and respiration is defined by the carbon use 181 efficiency (CUE). 182 In turn, active SOC is degraded into both slow and passive SOC and the respiration fluxes associated with these 183 processes are also controlled by the CUE (equations 6 and 7, with k SOC as the kinetic constant for SOC 184 decomposition) and ω SOC as the fraction of SOC that is transformed into DOC): 185 https://doi.org/10.5194/esd-2021-44 Preprint.  (Smith et al., 2010) and JULES (Nakhavali et al., 2018) that we followed here, where 194 the major exchange of C is between the different litter and SOC pools, and the production of DOC is related to 195 these SOC and litter pools by empirical rate constants, which were fitted to reproduce observed DOC turnover 196 times (Kalbitz et al.,2003, Turgeon, 2008 and DOC concentrations in the soil. The much higher DOC 197 production rates simulated by ORCHIDEE-SOM in its original configuration during preliminary tests over 198 Europe led us to implement the new approach (equations 4-7). While preserving the basic structure of 199 ORCHIDEE-SOM, we thus adapted the model in a way that organic C exchange occurs mainly among the 200 particulate litter and SOC pools, similar to the original Century model, while the production of DOC is 201 dependent on leaching rates as used in ECOSSE. In the modified soil carbon module, we used the parameter ω in 202 Equations (4-7) as a scaling factor that determines how much DOC is produced by the decomposition of litter 203 and SOC. This parameter was calculated after Smith et al. (2010) as the ratio of production of DOC from litter 204 ( ) and the SOC pools ( ) to the decomposition rates of litter ( ) and SOC ( ). The initial values 205 for ω were 0.5 % and 3 % for the litter and SOC pools, respectively. Further optimization with regard to 206 reproducing observed soil DOC concentrations led to ω values set at 0.2% for the litter and 1.2% for the SOC 207 pools. 208 Once produced, DOC can then adsorb on particles, and desorb from the adsorbed pool back into free DOC pool 209 following an equilibrium reaction between the dissolved and adsorbed phases. The partitioning is controlled by 210 K D , the so-called equilibrium partition coefficient (equation 9). 211 212 (9) Finally, the DOC pool is subject to decomposition according to equation (10) and then partly feeds into the SOC 213 pools (eq. 11). 214 215 (10) (11)

Manure as an additional C source 216
In Europe, a large fraction of the landscape is dominated by agricultural and grazing activities and manure 217 application represents a significant additional C -in particular DOC-source to the soil in regions dominated by 218 grasslands and croplands. To constrain the C flux from manure infiltrating into the soil, we used the gridded 0.5° 219 resolution input of manure nitrogen (manure-N) applications produced by Zhang et al. (2017) as forcing file.

220
Following the use of that forcing data in the model branch ORCHIDEE-CNP developed by Sun et al. (2021), we 221 assumed that 90% of the total manure-N is in mineral form (i.e. NH 4 + or NO 3 -) and the remaining 10% is in 222 organic form. To convert the organic manure-N into a manure-C flux, a C:N stochiometric ratio of 13.7 was then 223 applied . Finally, the particulate and dissolved organic manure-C were assumed to feed 224 the litter and DOC pools, respectively (Fig. S1). Consistent with ORCHIDEE-CNP (Goll et al., 2017), the 225 fractions of particulate and dissolved manure-C were set to 0.9 and 0.1, respectively. 226 227

Hydrological processes 228
The representation of hydrological processes is handled in two distinct sub-modules. The first one, the hydrology 229 sub-module, simulates the vertical exchange of water in the atmosphere-vegetation-soil system in each model 230 grid cell, while the second one, the river routing module, simulates the horizontal transfers between grid cells. 231 The hydrology is forced by several meteorological fields such as precipitation and air temperature. In the 232 hydrology module, precipitation is divided into interception and throughfall, the latter being further subdivided 233 into surface runoff and infiltration into the soil. The infiltration rate is controlled by the throughfall rate, the 234 slope of the soil surface and the hydraulic conductivity of the soil which is a limiting factor for infiltration. The 235 distribution of water within the soil is represented by the distribution of soil moisture over the discretized soil 236 profile (de Rosnay et al., 2002, d'Orgeval et al., 2008. The water budget within the soil is thus determined by 237 the infiltration rate, the evaporation and transpiration from the soil, and drainage at the bottom of the soil 238 column. The infiltration rate and percolation through the soil profile are then used to compute the advective flux 239 of DOC (equation 2) 240 The second module deals with river routing and represents the horizontal transfers of water from the soil column 241 to the aquatic system though surface runoff and drainage, and further through the river network and adjacent 242 floodplains (Vorosmarty et al., 2000). Note that ORCHILEAK simulates the occasional inundation in the river's 243 floodplains, where decomposition rates of the different carbon pools (litter, SOC and DOC) differ depending on 244 whether the soil is flooded or not. For a detailed description of its features, please refer to Lauerwald et al. 245 (2017).

Model set-up 249
Model domain, land cover and forcing data. The simulated model domain extends over the area covered by 250 the EU-27 (4.1 10 6 km 2 ) between 35°N and 70°N latitude and 10°W and 30°E longitude (Fig. 2). This domain 251 includes 5600 model grid cells at 0.5x0.5° resolution and encompasses 6 broad climate zones according to the 252 Köppen-Geiger classification from Peel et al., 2017 (Fig. 2a). The dominant PFTs within Europe include 253 croplands (20% mainly C3), grasslands (31% of which 24% are C3), and forests (25% of boreal forest, of which 254 16% are needleleaved evergreen and 9% are broadleaved summer-green, and 14% of temperate forest, of which 255 8% are broadleaved summer-green, 3% are needleleaved and 3 % are broadleaved evergreen) (Fig.2b). The 256 spatial distribution of manure application on grasslands and croplands is shown in Fig.2c. Finally, Fig.2d  257 illustrates the actual river network derived from the HydroSHEDS DEM data (Lehner et al., 2008) and the one 258 corresponding to our river routing scheme at 0.5 degree resolution, highlighting that the representation of the 259 river network is not optimal due to the coarse spatial resolution of our model. This coarse resolution limits the 260 possibility of model validation to the downstream parts of larger river networks. Note further that the mouth of 261 the Rhine is more than 100 km too far east, which further limits model validation for that river.

268
The forcing data applied in our study are listed in table 1. They are the same as those used in Lauerwald et al. 269 (2017) except for the meteorological forcing data and the land cover. The WFDEI meteorological forcing dataset 270 used in this study was derived by applying the methodology originally used to create the WATCH Forcing Data 271 (WFD) from the ERA-Interim reanalysis data (Weedon et al., 2014). The dataset has a 0.5° spatial resolution and 272 a 3-hourly time step from 1978 to 2014. The land cover forcing data set, which gives the areal proportion of the 273 13 PFTs within each 0.5° grid cell, was taken from Peng et al. (2017). 274 275  ground temperatures and thus lower evaporation rates. This adjustment was deemed necessary in order to better 283 capture the observed mean and seasonal variability of the discharge along the European river network. 284 Spin-up. Before the model can be used to simulate C dynamics over the past decades, a spin up is needed to 285 reach an assumed steady state for the C fluxes during the pre-industrial period. This steady state is achieved by 286 spinning up ORCHILEAK for 15000 years. The spin up was realized by recursively looping over 27 years of 287 climate forcing using the WFDEI forcing dataset over the 1980-2006 period and constant land cover and 288 atmospheric CO 2 concentration of 286 ppm (Guimberteau et al., 2018) corresponding to year 1861. After the end 289 of the spin-up, the soil C stock across the entire European continent changed by less than 1% over a century of 290 simulation, which we considered close enough to steady-state. 291 Transient runs. Using the steady-state outputs as initial condition, the first part of the transient simulation 292  was carried out with increasing atmospheric CO 2 concentration, changing land use and land cover 293 and with river routing activated while still looping over the same 27 years of climate forcing used for the spin-up 294 because no climate forcing data prior to 1978 was available. From then onwards, the WFDEI atmospheric 295 forcing data was applied over the entire period covered by this product . 296 297

Model evaluation 298
Firstly, the simulated discharges were compared to times series of daily stream flow recorded at eleven gauging 299 stations from "The Global Runoff Data Center (GRDC), 56068 Koblenz, Germany" dataset. For comparison, 300 both observed and simulated discharges were aggregated at the monthly temporal resolution over the years 1980 301 to 2006. This period was chosen based on the GRDC data coverage. Model performance was then evaluated 302 with respect to several variables of the terrestrial C cycle. Firstly, simulated Net Primary Production (NPP) was 303 compared to two different data products. The first one, the CARbon DAta MOdel fraMework (CARDAMOM; 304 Bloom et al., 2015)  basins for which we also evaluated the simulated river discharge and DOC fluxes, and which taken together, 316 represent 19 % of the model domain ( System (GIS) that contains up-to-date information on world soil properties. In particular, this dataset reports the 336 organic carbon content in the soil as well as the soil bulk density. The bulk density was calculated in two 337 different ways. The first one follows the method described in Saxton (1986) where the bulk density is related to 338 the soil texture, this approach tending to overestimate density in high porosity soils or in OC rich soils. The 339 second method uses the SOTWIS database in which the bulk density is calculated as a function of soil type and 340 https://doi.org/10.5194/esd-2021-44 Preprint. Discussion started: 22 June 2021 c Author(s) 2021. CC BY 4.0 License. depth. In this database, all variables are reported for the topsoil (0-30cm) and the sub-soil (30-100cm) horizons. 341 For comparison purposes, our simulated SOC stocks were thus integrated over the same depth intervals. We also 342 assessed briefly the extent to which our model can reproduce the main features in observed soil DOC profiles. 343 To that end, we compared our simulated DOC profile averaged over the entire European forest biome against the Those stations are located near the mouth of large rivers (Danube, Rhine, Rhone, Elbe and Seine) but also 366 include a few locations further upstream the same rivers or at major tributaries (Fig.3). The comparison is 367 performed for the period 1990-2000, except for the Rhone at Beaucaire and the Danube at Svistov for which the 368 observed stream gauge data cover only a shorter period. Overall, the model reproduces the observations well, 369 both in terms of amplitude and seasonality, except for the Elbe at Neu Darchau for which the temporal variability 370 is well captured but the absolute discharge is overestimated. 371 Note that the simulated catchment area often diverges (by -25% to +30 %) from the observed value due to the 372 coarse resolution (0.5°x0.5°) of ORCHILEAK (Table S3) Furthermore the 0.5° resolution is too coarse to be able to represent perfectly the pathways of the river. Our 377 model tends more often to underestimate the catchment area, while its yearly mean discharge is overestimated, 378 except at the Beaucaire station along the Rhone River. The bias can be significant and cannot be explained by 379 the model resolution alone. 380 To evaluate model performance for discharge, we used the Pearson's coefficient of determination (R 2 ) and the Nash Sutcliffe modeling efficiency (NSE, Nash and Sutcliffe (1970)). The R 2 only accounts for the correlation with regard to the temporal variability. With R 2 values comprised between 0.43 and 0.62 for all stations, we conclude that the observed seasonality of the discharge along large European rivers is reasonably well reproduced by the model. The NSE not only accounts for the correlation between observed and simulated temporal signals, but also for the model's ability to reproduce absolute 390 discharges. The statistics confirm our previous observation that the model generally overestimates discharges (low NSEs) except for stations Elbe in Dresden, Rhone in Beaucaire, Rhine in Basel and Danube in Bratislava where both the mean and temporality are well captured. Two stations have negative NSE values, which means that the error variance estimated by the model is significantly larger than the variance of the observations; in others words, the difference between model and observations is significant. The mean error (%), that is the weighted difference between the average from the model and the 395 one from observation, confirms that low NSEs are mostly due to overestimated discharges, which is further demonstrated by high mean errors. More results for other European catchments can be found in table S3.

NPP, biomass and soil organic C stocks
We briefly compare simulated NPP with the gridded observation-based products GIMMS and CARDAMON (section 2.2.2) 400 as C fixation by the vegetation exerts an important control on DOC stocks in the soil and thus on DOC leaching. We first perform our comparison over a large domain comprised between -10° and 30° in longitude east and 35° and 70° in latitude north -covering the area from Ireland to the Western Black Sea (where the Danube flows into) and from the south of Spain to the north of Scandinavia. Over this area (referred to as "Europe" from here onwards), the modeled yearly averaged NPP amounts to 445 gC m -2 yr -1 , a value in remarkable agreement with both GIMMS and CARDAMOM estimates of 430 gC m -2 405 yr -1 and 460 gC m -2 yr -1 , respectively. Those two datasets are entailed with an uncertainty that we assume similar to that reported for the MODIS dataset, i.e. 20% (Turner et al 2006). The total living biomass in Europe is simulated at 15.5 PgC or 2.3 kgC m -2 . This value is in good agreement with the recent estimate by Avitabile and Camia 2018, which report a biomass stock at around 16 PgC. We estimate that the total soil carbon stock amounts to 58 PgC. Averaged over the first meter of the soil horizon, this corresponds to a value of 9.5 kgC m -² which is comparable to that of HWSD (6 kgC m -²) when using the 410 SOTWIS method to compute the bulk density, but significantly lower when applying Saxton's method (22 kgC m -²), plausibly because the latter overestimates the bulk density in OC-rich soils (Kochy et al., 2015). Table 2 summarizes the yearly average NPP at the scale of the five selected European catchments. Simulated NPP is of the same order of magnitude as both observation based datasets, without any systematic bias towards an underestimation or overestimation. To provide error bounds for the observational products, we calculated the average standard deviation 415 between yearly-mean values. For GIMMS, we also included the standard deviation induced by the use of the five distinct meteorological forcing files to assess the NPP (section 2.2.1). We find that our simulated catchment averaged NPP fall within the error bounds of the observational products for the Rhine and the Rhone while for the Danube, Elbe and Seine, simulated NPP is slightly above the upper error range.   GIMMS (1982GIMMS ( -2006 datasets. The mean of the two datasets, along with an assessment of the uncertainty (based on MODIS) and of the standard 430 deviation are also reported. In addition, the modeled biomass stock and soil organic carbon (SOC) content (first 1m) are compared with values reported in the HWSD database, using two methods (Saxton and SOTWIS) to calculate the soil bulk density. All variables and processes are reported for the large-scale basins of focus in this study (see fig. 3 for location)

Soil DOC stocks
Comparison between observed and modeled DOC stocks and fluxes is more difficult than for biomass and SOC because those have not been assessed at large spatial scales. Nevertheless, representative soil DOC concentration profiles for 440 coniferous and broadleaved forests of Europe have been compiled by Camino et al. (2014). These profiles were used to evaluate our model. Overall, we found that ORCHILEAK slightly overestimates DOC concentrations, especially in the very topsoil horizons with modeled values around 100 mg l -1 against 40-60 mg l -1 in the observations (Fig. 5). We also simulated higher concentrations in broadleaved forests than in coniferous forests while Camino et al. (2014)

DOC leaching fluxes
The model simulates a yearly-mean DOC leaching flux over Europe of 14.3 (±10) TgC yr -1 (Fig 6), the standard deviation 460 being here coarsely approximated by spatial variability. The average area specific flux rates is of 2.6 (±2.5) gC m -2 yr -1 . We compared DOC leaching fluxes with site level observations from Kindler et al. 2011, across 17 local measurements, each sampled fortnightly during the period October 2006 until March 2008. Our modeled average of 2.6 (±2.5) gC m -2 yr -1 is of the same order of magnitude as the observed one (4.2 gC m -2 yr -1 ). Although the modeled mean is about 38 % lower than the one measured, the standard deviation representing the spatial variability in simulated DOC leaching fluxes over all our 465 model grid cells encapsulate the observational mean, highlighting a significant heterogeneity that is difficult to embrace with local measurements alone.

Fluvial DOC decomposition and export fluxes
The export of DOC from the European river network to the coast is arguably the best monitored variable against which our model can be evaluated. Using this flux to build confidence in our estimate of the terrestrial DOC leaching requires an assessment of the DOC degradation within rivers, a process that is controlled by the hydrology and the half-lives of reactive DOC compounds. In the model, the first-order decomposition rates at a given temperature of 28°C are equal to 0.3 d -1 and 490 0.01 d -1 for the labile and refractory DOC pools, respectively. Based on those values and the simulated distribution of labile and refractory DOC, the estimated bulk decomposition rate constant averaged over the entire model domain is equal to 0.05 d -1 , which corresponds to a half-life for DOC of about 14 days (Table 3). This rate constant varies across Europe but always remains within the same order of magnitude, with half-lifes ranging from 6 to 20 days (0.035-0.122 d -1 ). These     figure 9). In addition, Abril et al. (2002) report DOC 510 concentrations measured at nine river mouths discharging along the Atlantic façade and the North Sea, three of which (Rhine (NL), Scheldt (BE) and Gironde (F)) resolve the seasonality while the other six (Elbe (GE), Ems (GE), Thames (UK), Loire (FR), Sado (P), Douro (P)) only rely on a single measurement per year. Both GLORICH and Abril et al. (2002) report DOC concentrations at the mouth of the Rhine and the Elbe but their values diverge because in addition to analytical uncertainties, the sampling period and data density are not the same. Measured values are equal to 4.3 and 2.9 mg C l -1 for the Rhine and 515 4.6 and 6.1 mg C l -1 for the Elbe, respectively highlighting inherent variability in observational data. To complement these local samplings, we also compared our simulated DOC concentrations with those of Mattson et al. (2008) for several groups of catchments in Finland (9 spread over the whole country), Denmark (10 draining into Horsens fjord), the UK (10 draining into the River Conwy) and France (5 draining into the River Tech). All measured DOC concentrations ranged from 2.5 mg C https://doi.org/10.5194/esd-2021-44 Preprint. Discussion started: 22 June 2021 c Author(s) 2021. CC BY 4.0 License. l -1 to 10 mg C l -1 except in two regions in the north (Finland and basins flowing into the Baltic sea) where concentrations 520 exceeded 10 mg C l -1 . For most of the data, the model slightly overestimated the river DOC concentrations. The model results also suggest that the concentrations broadly increase with latitude, with the higher values found in humid continental and subarctic climate and the lower ones in the Mediterranean climate, a result in agreement with the observations from Mattson et al. (2008). Such pattern possibly results from decreasing mean annual air temperature and runoff in Northern Europe that favor incomplete decomposition of litter and DOC, thus favouring DOC production, while at the same time 525 DOC turnover rates in the soils are decreased. Also the increased abundance of forests, and in particular coniferous forests, is a valid explanation for higher DOC leaching (Lauerwald et al. 2012). However, it is important to keep in mind that we're missing the peatlands, suggesting that we could lack part of the DOC leaching in subarctic and tundra regions leading to even higher DOC fluxes further in the North. Finally, the comparison reveals that model performance tends to improve with catchment size, likely reflecting the difficulty to capture the DOC dynamics at the small scale with the current resolution of 530 ORCHILEAK. But overall, our model is capable of reproducing observed yearly mean DOC concentrations for a wide range of river basins spread between Finland and Portugal. The temporal evolution of observed river DOC fluxes is only available at four stations (Rhine, Elbe, Rhône and Seine) where DOC time series have been recorded over multi-annual periods (Fig. 9). The multi-year mean modeled DOC fluxes are estimated for the Rhine, Elbe, Rhone and Seine at 11.9, 7.2, 8.8 and 3.2 kg s -1 , respectively. The observations amount 540 respectively to 7.9, 3.6, 4.6, 1.6 kg s -1 . For all stations, the model thus overestimates slightly fluvial DOC fluxes, which is not surprising since the model tends to overestimate the discharge. At these four stations, ORCHILEAK also slightly overestimates river DOC concentrations except for the Seine where concentrations are largely underestimated and discharge largely overestimated. In terms of temporal correlation, the simulated DOC flux for the Rhone compared to the observed one yields a R 2 of 0.6 and a mean error of 92% (results for the Seine, Elbe and Rhine are reported in supplementary table S5). In Finally, although the model-data comparison points to a slight overestimation of the river DOC export flux, our pan-555 European estimate amounts to 12.3 TgC yr -1 . This estimate is in fact about 35 % lower than the one reported in another model study by Li et al. (2019), based on the TRIPLE-HYDRA, a process-based model for which the DOC export flux reaches 19.3 TgC yr -1 . Li et al. (2019) applied the model at the global scale and simulation results were primarily evaluated against observations in the world-largest rivers and for Europe only included the Volga River. Li et al. (2019) then applied the model for multiple rivers in Europe such as the Danube, the Po, and the Elbe. Despite these different scales of analysis, 560 the export fluxes predicted by both models fall within the same order of magnitude.

Manure implementation
The implementation of manure significantly affects DOC leaching from grasslands and croplands (Fig. 10) which cover more than half of the studied region. The average annual input rate of manure into the soil is around 2.5 gC m -2 yr -1 (Fig. 2c). 565 With manure implementation, the DOC leaching rate increase drastically (average of +72% compared to the DOC leaching without manure), in particular in the oceanic and humid continental climate regions, where the average DOC leaching rate changes from 1.6 to 2.7 gC m -2 yr -1 and 1.7 to 2.5 gC m -2 yr -1 , respectively. In whole Europe, manure implementation leads to an increase of total DOC leaching into the river network from 9.8 to 14.3 TgC yr -1 .

Drivers of DOC leaching 575
Here, we analyze what controls the spatial distribution and temporal variability in the DOC leaching. Figure 11 shows seasonal variability of DOC leaching and total runoff (surface runoff plus drainage) in different climate zones of Europe, revealing a clear and consistent relationship between those two fluxes. The seasonal peak in DOC leaching consistently occurs in winter while minimum values are found during summer. These results suggest that both spatial and temporal variability in leaching are correlated to total runoff. 580 To further explore the environmental controls the DOC leaching, we expressed DOC leaching as fraction of terrestrial NPP 585 ( Figure 12). Doing this, we assume that NPP, which is undoubtedly the first source for DOC production, is as well an important control of the DOC leaching flux. Moreover, normalizing DOC leaching by NPP, we strive to show the possible influence of other controls, allowing for a more in-depth analysis of the effect of hydrology and climate on the DOC leaching flux. Figure 12 reveals that the fraction of terrestrial NPP lost to DOC leaching increases, as expected, with total runoff. Moreover, this fraction increases with the contribution of surface runoff to total water loss from surface runoff plus 590 drainage (Fig. 12b). This can be explained by the general decrease in soil DOC concentrations with depth (Fig. 5), leading to higher DOC concentrations in surface runoff than in drainage. In fact, according to our simulations, 97% of the leached DOC is concentrated in the surface runoff. Note that higher total runoff is often associated with a higher contribution of surface runoff, which leads to a 'flushing effect' where high runoff events contribute a disproportionate high fraction of the longterm DOC leaching (Idir et al. 1999, Raymond andSaiers 2010). Finally, we found higher leaching to NPP ratios at lower 595 temperatures (Fig. 12a), hinting at the fact that lower temperatures lead to longer turnover times of DOC in the soil, and thus

605
To better quantify the effects of all these drivers on DOC leaching, we fitted a multi-linear regression model to predict the ratio of DOC leaching to NPP as a function of surface runoff, drainage and temperature at all grid points and for each month over the simulation period (eq. 13). To compare the importance of each predictor for the spatiotemporal patterns of DOC leaching, we normalized all variables Vs of equation 13 according to equation 14 (where i is the cell index).
(p-value < 2*10 -16 except for temperature where p-value = 2.7*10 -5 ) 610 (14) To rule out any significant multi-collinearity in the regression model, we calculated for each predictor the Variance Inflation Factor (VIF). The VIF evaluates the correlations among all predictors which could impact the robustness of the regression model (James et al., 2017). The closer the VIF is to 1, the more robust is the model. In our regression, VIF's of the runoff, https://doi.org/10.5194/esd-2021-44 Preprint. Discussion started: 22 June 2021 c Author(s) 2021. CC BY 4.0 License. drainage and temperature are respectively 1.13, 1.13 and 1.01, confirming that our prediction is robust and not biased by high multicollinearity. 615 In Fig. 12 (c), the DOC leaching simulated by ORCHILEAK is compared with the one predicted by equation 13. Our simple regression model is able to reproduce the simulations with a residual standard error of 0.68% and a R² of 0.45. The coefficients of our regression model reveal that spatio-temporal variability in DOC leaching is mainly driven by the surface runoff (K R ) and drainage (K D ). Air temperature as third control of DOC leaching is of subordinate importance as reflected by its low predictor's coefficient (K T ). 620 Table 4 summarizes for each climate zone in Europe the DOC leaching fluxes, in total numbers and normalized by NPP, as well as other important components of the terrestrial C budget. Since runoff and temperature were identified as the controlling factors of the DOC leaching flux, normalized DOC leaching fluxes are expected to be significantly different among climate zones. Indeed, the fraction of NPP lost to the river network as DOC is the lowest in the semi-arid region (0.13%) where annual precipitation is low (total runoff around 92 mm per year) and temperatures are high. 625 The highest fraction of NPP exported to rivers as DOC is found in the tundra climate and reach 1.22%. That can be 630 explained by high runoff and drainage (reaching 920 mm per year) in this climate zone, but also by low temperatures lowering the fraction of DOC already decomposed within the soil column. The subarctic climate also presents a similarly high DOC leaching to NPP ratio with a value of 0.84%. The Mediterranean, Oceanic and humid continental climate zones present intermediate DOC leaching to NPP ratios of respectively 0.26%, 0.48% and 0.49%. Averaged over the whole of the EU-27, the DOC leaching flux normalized to the NPP amounts to 0.60 %. 635

Comparison with previous assessments of DOC leaching
In one of the first studies on the terrestrial C budget of Europe (Janssens et al., 2003) an imbalance (missing sink) between atmospheric CO 2 inversions and bottom up C stock change accounting was partly attributed to the loss of carbon from land to rivers in the form of DOC of around 4 gC m -2 yr -1 . Our results, 2.6 ± 2.5 gC m -2 yr -1 , support this hypothesis although we 640 suggest a DOC leaching rate slightly lower than this early study. Our lower value may come from the fact that we did not simulate peatlands and organic soils which are known hotspots of DOC leaching (Leifeld and Menichetti 2018), in particular in areas such as the northern UK and Scandinavia. In terms of temporal variability, we found the highest DOC leaching in winter averaged over the continent (8.9 TgC in total for the six months of winter October to March) and the lowest in summer (5.4 TgC over the period April to September), consistent with the findings of Kindler et al. (2011). In terms of 645 drivers of the DOC leaching fluxes, our results are in line with empirical findings by Gielen et al. (2011) that identified hydrology as the main driver of the inter-and intra-annual variability in DOC leaching. Similar conclusions have also been drawn by other empirical studies (Michalzik et al., 2001, Neff and Asner 2001, Worrall and Burt 2007.
It is also interesting to compare our results with recent global and regional model studies of DOC leaching in tropical and boreal ecosystems. For the Amazon and Congo basins, Hastie et al. (2019Hastie et al. ( , 2021 found that 12 and 4 % of the NPP is 650 exported each year to inland waters in the form of DOC, respectivelymuch higher than the one we report for Europe as a whole (0.6%). Note that for these tropical lowland river basins extensive riparian wetlands are an important source of DOC, which are of lower importance in Europe. For the Lena river basin located in the boreal region, Bowring et al. (2020) found a DOC leaching of NPP ratio of about 1.5%. In our model assessment, this ratio reaches a very similar value of 1.2% for the boreal portion of Europe. For the temperate zone, a ratio of 0.35% for the East Coast of the US can be calculated when 655 dividing the average DOC leaching flux of 2.7 gC m -2 yr -1 simulated by Tian et al., 2015 by the average NPP of 780 gC m -2 yr -1 estimated by Zhao et al., (2005). Further, our value is quite similar to the one extracted from the global study by Nakhavali et al. (2020) that amounts to 0.5 % for the European domain only. Overall, this comparison highlights that in Europe, the fraction of NPP lost as DOC to the river network is significantly smaller than in other regions of the world. The lower value is likely due to the lower connectivity between terrestrial and aquatic systems due to the lack of extensive 660 wetlands, which have been reduced by major regulation of the European river network.

3 Implications for the terrestrial carbon budget of Europe
The terrestrial carbon budget is controlled by NPP, heterotrophic respiration, crop and wood harvesting and land use change.
Here we look at the net ecosystem exchange (NEE) which is the net C exchange between land and atmosphere (Kramer et 665 al., 2002). However, this view neglects the leakiness of terrestrial ecosystems that permanently removes a fraction of the land C and export it to the river network. Moreover, we can argue that DOC leaching represents a fraction of NEE, while the remainder of NEE can be attributed to harvest, land use change and changes in biomass and soil C stocks. From 1979 to 2012, the average NEE in Europe is 860 TgC yr -1 (123 gC m -2 yr -1 ), equaling about 28% of the total NPP (Fig. 13b). The ratio of DOC leaching to NEE shows drastic spatial variation, varying from an average value of 0.4% in the semi-arid 670 regions to a value of 5.7% in the tundra. In whole Europe, the DOC leaching is about 3% of the NEE.  We reconstructed the terrestrial and riverine C fluxes in Europe during period 1979-2012 using the ORCHILEAK LSM. The total C leaching from soil to European rivers is 14.3 TgC yr -1 on average, about 0.6 % of the estimated NPP and 3% of the terrestrial net up-take of atmospheric C. This flux shows large spatial and temporal variations. In specific, DOC leaching overall increases from warm and dry regions to cold and wet regions. However, since the model does not represent peatlands yet, the simulation results for subarctic and tundra regions in northern Europe could be biased. In whole Europe, DOC 680 leaching rate is the highest in winter and lowest during the summer, mainly controlled by the seasonal variation of runoff.
The implementation of manure lead to a significant increase in DOC leaching over the oceanic and humid continental region where croplands and grasslands are dominant. Our results contribute to a better assessment of the land-ocean C fluxes in Europe and to a better understanding of the effects of lateral C transfer on the terrestrial C budget. Combined with recent large-scale studies in tropical and boreal biomes as well as along the east coast of the US, an emergent view regarding the 685 global role of DOC leaching on the terrestrial C balance and its underlying drivers is progressively emerging.

Code and data availability
The model code used in this study is available at DOI : 10.14768/75AC2F47-4691-46AF-9B12-B1A9629CBC56 All forcing data set are listed in table 1. Data of observed discharge used in this study are available from the Global Runoff 690 Data Center (GRDC) at www.bafg.de/GRDC. Data of observed DOC concentrations in France are provided by eau de France at http://www.data.eaufrance.fr/.