Spatiotemporal patterns and drivers of terrestrial dissolved organic carbon (DOC) leaching into the European river network

Leaching of dissolved organic carbon (DOC) from soils into the river network is an important component of the land carbon (C) budget. At regional to global scales, its significance has been estimated through simple mass budgets, often using multi-year averages of observed fluvial DOC fluxes as a proxy of DOC leaching due to the limited availability of observations of the leaching flux itself. This procedure leads to a systematic underestimation of the leaching flux because of the decay of DOC during fluvial transport. Moreover, this procedure does not allow for revealing spatiotemporal variability in DOC leaching from soils, which is vital to better understand the drivers of DOC leaching and its impact on the local soil C budget. In this study, we use the land surface model (LSM) ORCHILEAK to simulate the terrestrial C budget, including leaching of DOC from the soil and its subsequent reactive transport through the river network of Europe. The model performance is evaluated not only against the sparse observations of the soil DOC leaching rate, but also against the more abundant observations of fluxes and reactivity of DOC in rivers, providing further evidence that our simulated DOC fluxes are realistic. The model is then used to simulate the spatiotemporal patterns of DOC leaching across Europe over the period 1972–2012, quantifying both the environmental drivers of these patterns and the impact of DOC leaching on the land C budget. Over the simulation period, we find that, on average, 14.3 Tg C yr⁻¹ of DOC is leached from land into European rivers, which is about 0.6 % of the terrestrial net primary production (NPP), a fraction significantly lower than that reported for tropical river networks. On average, 12.3 Tg C yr⁻¹ of the leached DOC is finally exported to the coast via the river network, and the rest is respired during transit. DOC leaching presents a large seasonal variability, with the maximum occurring in winter and the minimum in summer, except for most parts of northern Europe, where the maximum occurs in spring due to snowmelt. The DOC leaching rate is generally low in warm and dry regions, and high in the cold and wet regions of Europe. Furthermore, runoff and the ratio between runoff from shallower flow paths on one hand and deep drainage and groundwater flow on the other hand are the main drivers of the spatiotemporal variation of DOC leaching. Temperature, as a major control of DOC production and decomposition rates in the soils, plays only a secondary role.

1 Introduction

Terrestrial ecosystems are an important carbon (C) sink as they absorb about one-fourth of anthropogenic CO₂ emissions and store this C in plant biomass and soil carbon pools (Friedlingstein et al., 2020). This terrestrial C sink mitigates the growth rate of atmospheric CO₂ concentration and thus plays an important role in regulating climate change (Ciais et al., 2013). However, the efficiency of that sink is partly alleviated by the permanent, lateral leaching of C from soils, through the river network down into the ocean (Regnier et al., 2013). An accurate understanding of lateral C fluxes through the river network is thus necessary to better understand the global C cycling and to inform policies of climate change mitigation (Le Quéré et al., 2018).

The identification of riverine C transfers as a key component of the continental C budget constituted an important paradigm shift in our understanding of the global C cycle (Cole et al., 2007). More recently, riverine C cycling was also shown to be affected by anthropogenic perturbation and thus to be an element of the anthropogenic CO₂ budget (Regnier et al., 2013; Le Quéré et al., 2015). Anthropogenic perturbations of riverine C fluxes are manifold and comprise direct impacts through changing C and nutrient inputs following land-use change and agricultural activities, wastewater discharge and hydraulic management (e.g., Tian et al., 2015; Lauerwald et al., 2020; Hastie et al., 2021; Maavara et al., 2017). There are also indirect impacts following climate change and changes in atmospheric composition. Together, these perturbations have accelerated the turnover of C along the terrestrial–inland water continuum. The terrestrial C sink, which is classically estimated without taking into account C exports through the river network, is thus generally overestimated (Regnier et al., 2013; Lauerwald et al., 2020).

The integration of riverine C transfers into the terrestrial C budget requires the quantification of the amount of C lost from soils to the river network. Due to the scarcity of observational data, this flux is not easy to estimate using empirical methods. At a global scale, this flux was constrained through budget closure based on estimates of riverine C exports to the coast and estimates of C losses to the atmosphere and aquatic sediments during transport. The existing global estimates of these soil C exports to the river network, as synthesized by Drake et al. (2018), range from 1.1 to 5.1 Pg C yr⁻¹ – a huge uncertainty range reflecting the limitations of empirical estimation approaches and the paucity of underlying data. Over the past decade, a new generation of land surface models (LSMs) have been developed that represent the export of C from soils to the river network, and in some cases even the transport and cycling of these terrestrial C loads along the river network down to the coast (Smith et al. 2010; Kicklighter et al. 2013; Tian et al., 2015; Lauerwald et al., 2017; Nakhavali et al., 2018). With the exception of the study of Tian et al. (2015), these studies all focus on the lateral export of dissolved organic C (DOC), which is a product of the incomplete decomposition of plant litter and soil organic carbon (SOC). These mechanistically based models allow prediction of the leaching of DOC in unmonitored regions and assessment of spatial and temporal variability which, to date, can only be poorly resolved using empirical methods. Moreover, these approaches link the C exports from soils to the river network to the terrestrial C cycle, and thus allow for directly assessing the role of these C exports in the terrestrial C budget, its perturbation through changes in land use, climate and atmospheric chemistry and its impact on the terrestrial sink for anthropogenic CO₂ emissions.

In this study, we use the LSM ORCHILEAK (Lauerwald et al., 2017), a branch of the IPSL-LSCE LSM ORCHIDEE (Krinner et al., 2005), to quantify the DOC leaching from soils and its effects on the terrestrial C budget in Europe. ORCHILEAK not only simulates the vertical C cycling between vegetation, soils and the atmosphere in response to the climate, atmospheric CO₂ concentrations and land-use change, but also represents the lateral exports of DOC from soils to the river network, as well as the reactive transport of this DOC through the river network. To our knowledge, only one study (Kindler et al., 2011) has estimated soil DOC leaching flux based on runoff and direct observations of DOC concentrations in soil water for various locations across Europe. Thus, this empirical assessment will be used to evaluate the simulated DOC leaching fluxes in this study. Further, we evaluate simulated against observed riverine DOC fluxes, which are obtained from different water quality surveys and scientific publications. Assuming a realistic representation of DOC reactivity in the river network, which is to be evaluated against observations as well, this model–data comparison of riverine DOC fluxes represents a valuable and additional possibility to assess the validity of simulated soil DOC leaching.

So far, ORCHILEAK has been successfully tested and applied on large, near-natural river systems such as the Amazon (Lauerwald et al., 2017, 2020; Hastie et al., 2019), the Congo (Hastie et al. 2021) and the Lena rivers, with a version also including some specific permafrost-related mechanisms (Bowring et al., 2019, 2020). In this study, for the first time, ORCHILEAK is applied to, and evaluated for, the European river network, which is subject to direct impacts of agricultural land use in contrast to more natural river basins. For this reason, we devote special attention to manure application as an anthropogenic non-point source of DOC for the river network, while we assume that for the period of simulation (1979–2012), due to the quality of sewage water treatment, anthropogenic point sources of DOC are now negligible for most parts of Europe. Moreover, as shown for instance by Meybeck (1986), DOC from sewage is highly labile and only affects concentration within short distances downstream of water processing plants. Avoiding observational data from sites that are known to be directly impacted by sewage inputs, we are able to evaluate model performance with regard to fluvial transfers of soil-derived DOC, which is the focus of our study.

Making full use of the capabilities of the ORCHILEAK model, we study in detail the spatiotemporal patterns in DOC leaching and its quantitative contribution to the terrestrial C budget across Europe. We investigate how specific climate zones in Europe differ with regard to seasonality in DOC leaching fluxes, which are hypothesized to be controlled by hydrology, litter fall and temperature effects on litter and SOC decomposition. We will further try to quantify the effect of these controls in the different climate zones of Europe. Finally, we strive to find out in which climate zone DOC leaching affects the terrestrial C budget the most.

2 Methodology

2.1 ORCHILEAK

2.1.1 Model overview

ORCHILEAK (Lauerwald et al., 2017) is a branch of the model ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosystems; Krinner et al., 2005), the land surface component of the Institut Pierre-Simon Laplace (IPSL) Earth system model (ESM). ORCHIDEE simulates energy, water and C fluxes between the atmosphere and the land at a global scale. This LSM is based on four sub-modules. The first one, SECHIBA, simulates the energy budget (energy, carbon and water) between the atmosphere and the biosphere, as well as hydrology, which in the default setup used here are both represented using a 30 min time step. The second sub-module, adapted from the LPJ model (Sitch et al., 2003), represents the dynamics of vegetation distribution on long timescales (1 year), while the third one (STOMATE) simulates the C dynamics in vegetation and soils at a daily to sub-daily step (Krinner et al., 2005). Finally, the fourth sub-module handles the routing of water that is lost via surface runoff and drainage from soils to the ocean through the global river network (Polcher, 2003; Guimberteau et al., 2012), for which a daily time step is used. All processes are simulated on a horizontal model grid, the resolution of which can be adapted to that of the meteorological forcing files. In this study, simulations are run at a spatial resolution of 0.5^∘. Moreover, in the default setup, up to 13 plant function types (PFTs; bare soil, eight types of forest, two types of grassland and two types of cropland) can be distinguished for each cell, for which C budgets are simulated individually, while energy and water budgets are simulated at the grid-cell level.

ORCHIDEE represents the soil C dynamics distinguishing different pools of plant litter and soil organic C over a 2 m profile. A branch of ORCHIDEE, called ORCHIDEE-SOM (Camino-Serrano et al., 2018), added a vertical discretization of these carbon pools over 11 layers and included the representation of DOC production and cycling within the soil column (see Sect. 2.1.2 for more details). ORCHILEAK was built on this branch and accounts for the coupled reactive transport processes impacting the dissolved C inputs from soils to the river network, including both DOC leaching from soils and CO₂ produced by soil respiration, into the hydrologic routing scheme. Besides advective transport of carbon with the water flow, ORCHILEAK simulates the decomposition of DOC during riverine transport, the gas exchange of CO₂ at the interface between the inland water and the atmosphere and the exchange of C between water column and soil column in inundated floodplains. For the representation of in-river DOC decomposition, two pools with different decay rates are distinguished, a slow (refractory DOC) and a fast (labile DOC) pool. All of these fluxes are closely coupled to the model representation of hydrology that comprises interception of precipitation, throughfall, infiltration, percolation, surface runoff, drainage and the routing of discharge along the river–floodplain network.

2.1.2 Soil carbon module

The soil carbon module of ORCHILEAK (Fig. 1) is derived from the CENTURY soil carbon model of Parton et al. (1988). In the standard scheme (Krinner et al., 2005), C in the soil of each model grid cell, and of each PFT, is represented by four different litter pools and three different soil organic carbon (SOC) pools with different turnover rates. The four litter pools correspond to metabolic aboveground and belowground litter, and structural aboveground and belowground litter (Fig. 1). The SOC is subdivided into three pools, active, slow and passive, which have different default decomposition rates that are further modified at each time step according to the evolving soil moisture and soil temperature. In the CENTURY scheme, C from the decomposed structural litter enters the active and the slow pools with the fraction allocated to each pool depending on lignin content of the litter, while the entire metabolic litter pool and the remaining part of structural litter are allocated to the active SOC pool. The SOC pools then feed into each other, with the main C flux going from active to slow and to passive in order to represent microbial decomposition of detrital organic matter, and from a small return flux of slow and passive C back to the active pool in order to implicitly represent the C supply in the form of dead microbial biomass.

Figure 1 The new version of the soil module of ORCHIDEE-SOM. The left box represents the discretization of the soil column and the transport processes between layers. The right box shows an overview of all the biogeochemical transformation processes that occur in each layer.

Camino-Serrano et al. (2018) updated this scheme with a vertical discretization of distinct SOC and litter pools over a 2 m soil profile represented by 11 layers, with geometrically increasing thickness from top to bottom (Fig. 1). Camino-Serrano et al. (2018) further developed the soil C module by including an explicit representation of the fate of DOC along this vertically discretized soil profile. Processes accounted for are DOC production from the decomposition of SOC and litter, decomposition of DOC within the soil, sorption/desorption of DOC onto/from mineral surfaces, vertical advection and diffusion of DOC through the soil column and lateral, advective leaching of DOC out of the soil profile, along with surface runoff (water flux from the topsoil surface) and drainage (water flux from the last layer soil at 2 m depth). In each soil layer, ORCHILEAK explicitly simulates the fresh litter input (depending on the simulated vertical root distribution), decomposition of each organic matter pool (including litter and SOC), C transformation between different organic matter pools (shown by the blue and green arrows between different pools in Fig. 1), C transport and diffusion between neighboring soil layers and the loss of DOC due to leaching. For a specific organic C pool at each time step, only a fraction of the decayed C is respired as CO₂ to the atmosphere (orange arrows in Fig. 1), the remaining being transferred to other organic pools (to mimic microbial growth and mortality).

The DOC dynamics in the soil are simulated according to Eq. (1), which accounts for the dynamic interplay between production, decomposition, transport and sorption/desorption processes along the discretized 2 m soil column. All processes are simulated using a 30 min time step in following order: firstly, production and decomposition of DOC are calculated, and DOC stocks for each layer and pool are updated accordingly. Secondly, vertical exchange of DOC between soil layers is simulated in two steps: first for the process of DOC advection with the flow of water through the soil column, then for the diffusion of DOC. Lastly, the export of DOC through leaching from topsoil and bottom soil with runoff and drainage, respectively, is calculated.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{{\mathrm{dDOC}}_i}{\mathrm{d}t}}=\sum \left({\mathrm{Production}}_i-{\mathrm{Decomposition}}_i\right)+F_{\mathrm{A},i}+F_{\mathrm{D},i}\end{array}$$

In Eq. (1), i stands for the index of each layer. Each layer is connected to the adjacent layers by advective (F_A) and diffusive (F_D) fluxes. The total DOC transport flux is made of an advective component (Eq. 2) computed as the product of the water flux Q and the concentration of free DOC in the water solution of the ith layer and a diffusive component that follows Fick's first law (Eq. 3):

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{A}}=Q\cdot {\mathrm{DOC}}_i,\text{(3)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{D}}=-D\cdot {\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{DOC}}{\partial z^{\mathrm{2}}}},\end{array}$$

where i stands for the ith layer, z is the depth along the discretized soil profile and D stands for the molecular diffusion coefficient of DOC, which is assigned a value of $\mathrm{1.06}\times {\mathrm{10}}^{-\mathrm{5}}$ m² d⁻¹ (Ota et al., 2013).

The advective export of DOC to the river network is proportional to the top (first five layers, 4.5 cm) and bottom (11th layer) DOC concentrations, corresponding to water loss fluxes associated with runoff (for the near surface) and drainage (for deep soil layer). Diffusion of DOC between adjacent soil layers is proportional to the gradient in DOC concentrations in the soil solution (Eq. 3), moving towards equilibrium. In addition, we apply Fickian-type transport to represent the effect of bioturbation on SOC profiles. In this case, the transport is represented similar to Eq. (3), but follows the gradient in SOC concentration relative to the volume of the soil layer. Representing bioturbation as a diffusion-like process is a common approach in LSMs with vertically discretized SOC scheme (Camino-Serrano et al., 2018). However, bioturbation is much slower than diffusion of DOC in the soil solution, with a diffusion coefficient $D=\mathrm{2.74}\times {\mathrm{10}}^{-\mathrm{7}}$ m² d⁻¹ (Koven et al., 2013) compared to $D=\mathrm{1.06}\times {\mathrm{10}}^{-\mathrm{5}}$ m² d⁻¹ (Ota et al., 2013). Therefore, bioturbation impacts the vertical SOC profile while it has only a marginal influence on the DOC dynamics. The right-hand side of Fig. 1 summarizes the set of production/decomposition processes that occur in each layer. During litter decomposition, a fraction of the C is directly emitted back to the atmosphere as CO₂, while the remainder feeds the active and slow SOC pools:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}}\mathrm{respiration}\left(\mathrm{Litter}\right)=(\mathrm{1}-\mathrm{CUE})\cdot k_{\mathrm{L}}\cdot \left(\mathrm{1}-{\mathit{\omega }}_{\mathrm{L}}\right)\cdot \mathrm{Litter},\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{Litter}\mathrm{decomposition}=\mathrm{CUE}\cdot k_{\mathrm{L}}\cdot \left(\mathrm{1}-{\mathit{\omega }}_{\mathrm{L}}\right)\cdot \mathrm{Litter},\end{array}$$

where k_L is the kinetic rate for the litter decomposition (dependent on soil moisture and temperature; Camino-Serrano et al., 2018) and ω_L is the fraction of litter that is channeled into DOC production (as opposed to particulate SOC). This approach of relating DOC production directly to the decomposing litter is (following the ECOSSE model; Smith et al., 2010) inspired by Nakhavali et al. (2018), and is a major modification compared to the previous version of soil DOC and particulate organic carbon (POC) cycling of Camino-Serrano et al. (2018). In Eqs. (4) and (5), the partitioning between SOC production and respiration is defined by the carbon use efficiency (CUE).

In turn, active SOC is degraded into both slow and passive SOC and the respiration fluxes associated with these processes are also controlled by the CUE (Eqs. 6 and 7), with k_SOC as the kinetic rate for SOC decomposition, which depends on soil moisture and soil temperature, and ω_SOC as the fraction of decomposed SOC that is transformed into DOC):

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}}\mathrm{respiration}\left(\mathrm{SOC}\right)=(\mathrm{1}-\mathrm{CUE})\cdot k_{\mathrm{SOC}}\cdot \left(\mathrm{1}-{\mathit{\omega }}_{\mathrm{SOC}}\right)\cdot \mathrm{SOC},\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{SOC}\mathrm{decomposition}=\mathrm{CUE}\cdot k_{\mathrm{SOC}}\cdot \left(\mathrm{1}-{\mathit{\omega }}_{\mathrm{SOC}}\right)\cdot \mathrm{SOC}.\end{array}$$

The decomposition of the litter and SOC pools produces a small amount of DOC according to Eq. (8). The DOC pool is thus fed by seven contributing sources, one for each of the four decomposing litter pools and three from the decomposing SOC pools:

$$\begin{array}{}\text{(8)}& \mathrm{DOC}\mathrm{production}=k_{\mathrm{L}}\cdot {\mathit{\omega }}_{\mathrm{L}}\cdot \mathrm{Litter}+k_{\mathrm{SOC}}\cdot {\mathit{\omega }}_{\mathrm{SOC}}\cdot \mathrm{SOC}.\end{array}$$

In ORCHIDEE-SOM (Camino-Serrano et al., 2018), all decomposed litter and SOC which are not respired as CO₂ were first fed into the DOC pools, and only upon the decomposition of that DOC could the non-respired fraction of the decomposed DOC feed the other SOC pools. Such a formulation is in contrast to the adaption of the RothC SOC model in ECOSSE (Smith et al., 2010) and JULES (Nakhavali et al., 2018) that we followed here, where the major exchange of C is between the different litter and SOC pools, and the production of DOC is related to these SOC and litter pools by empirical rate constants, which were fitted to reproduce observed DOC turnover times (Kalbitz et al., 2003; Turgeon, 2008) and DOC concentrations in the soil. The much higher DOC production rates simulated by ORCHIDEE-SOM in its original configuration during preliminary tests over Europe led us to implement the new approach (Eqs. 4–7). While preserving the basic structure of ORCHIDEE-SOM, we thus adapted the model in a way that organic C exchange occurs mainly among the particulate litter and SOC pools, similar to the original CENTURY model, while preserving the basic structure of ORCHIDEE-SOM. The production of DOC is represented as a side product of this C exchange between the pools of litter and SOC, with production rates as used in ECOSSE. In the modified soil carbon module, we used the parameter ω (Eqs. 4–7) as a scaling factor that determines how much DOC is produced by the decomposition of litter and SOC. This parameter was calculated after Smith et al. (2010) as the ratio of the production of DOC from litter ($p_{{\mathrm{DOC}}_{\mathrm{L}}}$) and SOC pools ($p_{{\mathrm{DOC}}_{\mathrm{SOC}}}$) to the decomposition rates of litter (k_L) and SOC (k_SOC). The initial values for ω were 0.5 % and 3 % for the litter and SOC pools, respectively. Further optimization with regard to reproducing observed soil DOC concentrations led to ω values set at 0.2 % for the litter and 1.2 % for the SOC pools.

Once produced, the free DOC can then be adsorbed to soil mineral particles, and the adsorbed DOC can again be desorbed and return to the free DOC pool following a linear adsorption isotherm as described in Neff and Asner (2001) and Wu et al. (2014). We assume that equilibrium between the dissolved and absorbed phases is instantaneous. Moreover, the work by Kothawala et al. (2008) showed that this approach performed fairly well compared to the more complex approach of using Langmuir equations. The partitioning is controlled by K_D, the so-called equilibrium partition coefficient (Eq. 9), considered to be constant at $\mathrm{8.05}\times {\mathrm{10}}^{-\mathrm{5}}$ m³ water per kilogram of soil (Moore et al., 1992). All constants used here are listed in Table S1 in the Supplement.

$$\begin{array}{}\text{(9)}& {\mathrm{DOC}}_{\mathrm{adsorption}\text{-}\mathrm{desorption}\mathrm{eq}}=K_{\mathrm{D}}\cdot \mathrm{DOC}\end{array}$$

Finally, the DOC pool is subject to decomposition according to Eq. (10) and then partly feeds into the SOC pools (Eq. 11), where k_DOC is the DOC decay rate, which also depends on soil moisture and soil temperature.

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}}\mathrm{respiration}\left(\mathrm{DOC}\right)=(\mathrm{1}-\mathrm{CUE})\cdot k_{\mathrm{DOC}}\cdot \mathrm{DOC}\text{(11)}& {\displaystyle }& {\displaystyle }\mathrm{SOC}\mathrm{production}\mathrm{by}\mathrm{DOC}=\mathrm{CUE}\cdot k_{\mathrm{DOC}}\cdot \mathrm{DOC}\end{array}$$

2.1.3 Manure as an additional C source

In Europe, a large fraction of the landscape is dominated by agricultural and grazing activities, and manure application represents a significant additional C – in particular DOC – source for the soil in regions dominated by grasslands and croplands. Studies have shown an increase in riverine DOC concentration related to manure application, with frequency and intensity of storm events in spring directly after manure application exerting an important control on the amounts of additional DOC leached into the river network (e.g., Royer et al., 2006; Delpla et al., 2011; Singh et al., 2014; Humbert et al., 2020). The type of manure input can be very different from one region to another, and the physicochemical properties (e.g., C:N ratio and the ratio of dissolved and particulate organic matter) depend strongly on the specific type of manure input. However, current forcing data only provide the amount of total manure inputs, without information regarding specific composition and/or physicochemical properties of the manure. To constrain the C flux from manure infiltrating the soil, we used the gridded estimates of manure nitrogen (manure-N) applications produced by Zhang et al. (2017) as forcing file. Following the use of that data as forcing in the model branch ORCHIDEE-CNP developed by Sun et al. (2021), we assumed that 90 % of the total manure-N is in mineral form (i.e., ${\mathrm{NH}}_{\mathrm{4}}^{+}$ or ${\mathrm{NO}}_{\mathrm{3}}^{-}$) and the remaining 10 % is in organic form. To convert the organic manure-N into a manure-C flux, a C:N stoichiometric ratio of 13.7 was then applied (Vuichard et al., 2019). Finally, the particulate and dissolved organic manure-C were assumed to feed the litter and DOC pools, respectively (Fig. S1 in the Supplement). Consistent with ORCHIDEE-CNP (Goll et al., 2017), the fractions of particulate and dissolved manure-C were set to 0.9 and 0.1, respectively.

2.1.4 Hydrological processes

The representation of hydrological processes is handled in two distinct sub-modules. The first one, the hydrology sub-module, simulates the vertical exchange of water in the atmosphere–vegetation–soil system in each model grid cell, while the second one, the river routing module, simulates the horizontal transfers between grid cells. The hydrology is forced by several meteorological fields such as precipitation and air temperature. In the hydrology module, precipitation is divided into interception and throughfall, the latter being further subdivided into surface runoff and infiltration into the soil. The infiltration rate is controlled by the throughfall rate, the slope of the soil surface and the hydraulic conductivity of the soil which is a limiting factor for infiltration. The distribution of water within the soil is represented by the distribution of soil moisture over the discretized soil profile (de Rosnay et al., 2002; d'Orgeval et al., 2008). The water budget within the soil is thus determined by the infiltration rate and runoff from the topsoil, evaporation and transpiration from the soil and drainage at the bottom of the soil column. The infiltration rate and percolation through the soil profile are used to compute the advective flux of DOC (Eq. 2).

The second module deals with river routing and represents the horizontal transfers of water from the soil column to the aquatic system though surface runoff and drainage, and further through the river network and adjacent floodplains (Vörösmarty et al., 2000). Meybeck (1993a, b) found that soils were the major source of DOC for rivers, while autochthonous DOC was negligible at the global scale. Furthermore, autochthonous DOC has a short turnover time and is quickly recycled within the river (Farjalla et al., 2009; Fonte et al., 2013) and, therefore, does not contribute significantly to the net C budget of an entire river system. Since in this study the focus is on the role of fluvial DOC fluxes in the terrestrial C budget, autochthonous DOC is not accounted for. The representation of leaching processes in ORCHILEAK is simplified; the model represents flows of water from land to the stream network only through surface runoff and drainage from bottom soil. Leaching thus occurs either from the topsoil, which in our configuration represents the top 4.5 cm of the soil column, or from the bottom soil, i.e., the lowest 50 cm of the 2 m soil profile. DOC leaching from the topsoil is further limited by two reduction factors, a “general reduction factor” and a “reduction factor for export with runoff depending on the extents of saturated soils around headwaters”. The general reduction factor accounts for the fact that some of the runoff represents excess throughfall that never entered the soil and further corrects for the overestimated DOC concentration in the topsoil. The second reduction factor represents the connectivity between streams and their catchment through the extent of the water saturated riparian zone at sub-grid level, which varies with the amount of runoff and water stored in the “slow reservoir” of water, mimicking groundwater storage. Note that for larger river stretches (stream order 4 and higher), ORCHILEAK simulates the occasional inundation of the river's floodplains, where DOC production from the decomposition of litter and SOC in the topsoil, although reduced under inundated conditions, feeds directly into the DOC pool of the water column. For a detailed description of these features, please see Lauerwald et al. (2017).

2.2 Simulations

2.2.1 Model setup

Model domain, land cover and forcing data

The simulated model domain extends over the area (4.1×10⁶ km²) between 35 and 70^∘ N and 10^∘ W and 30^∘ E (Fig. 2). This domain includes 5600 model grid cells at 0.5×0.5^∘ resolution and encompasses six broad climate zones according to the Köppen–Geiger classification from Peel et al. (2017; Fig. 2a). The dominant PFTs within Europe include croplands (20 %, mainly C3), grasslands (31 %, of which 24 % are C3), and forests (39 %, of which 16 % and 9 % are needle-leaved evergreen and broadleaved summer-green boreal forests, respectively, while temperate broadleaved summer-green, needle-leaved and broadleaved evergreen forests comprise 8 %, 3 %, and 3 %, respectively; Fig. 2b). The spatial distribution of manure application on grasslands and croplands is shown in Fig. 2c. Finally, Fig. 2d illustrates the actual river network derived from the HydroSHEDS DEM data (Lehner et al., 2008) and the one corresponding to our river routing scheme at 0.5^∘ resolution, highlighting that the representation of the river network is not optimal due to the coarse spatial resolution of our model. This coarse resolution limits the possibility of model validation to the downstream parts of larger river networks. It is worth noting that the mouth of the Rhine is more than 100 km too far east, which further limits model validation for that river.

Figure 2 Spatial distribution for each 0.5^∘ grid cell of the continental European domain of (a) climate zones (according to the Köppen–Geiger classification); (b) dominant plant functional types (PFT); (c) manure application (in g C m⁻² yr⁻¹); (d) the routing network of ORCHILEAK (in blue). The real river network extracted from the European Environment Agency (https://www.eea.europa.eu/legal/copyright, last access: 3 March 2021) is also shown.

The forcing data applied in our study are listed in Table 1. They are the same as those used in Lauerwald et al. (2017) except for the meteorological forcing data and the land cover. The WFDEI meteorological forcing dataset used in this study was derived by applying the methodology originally used to create the WATCH Forcing Data (WFD) from the ERA-Interim reanalysis data (Weedon et al., 2014). The dataset has a 0.5^∘ spatial resolution and a 3-hourly time step covering from 1978 to 2014. The land cover forcing data set, which gives the areal proportion of the 13 PFTs within each 0.5^∘ grid cell, was taken from Peng et al. (2017). Note that the soil hydrology model in ORCHIDEE, which we adopted for ORCHILEAK, was developed and calibrated to work with the soil classes used in Reynolds et al. (1999). We thus kept this data source, while additional soil properties such as pH, and specific soil classes that we defined as “poor soils” (Histosols, Podzols) possessing lower C turnover times and DOC filtering, were taken from HWSD v1.1. A topographic index, which in ORCHIDEE controls the flow velocity in the river network of each cell, is taken from Vörösmarty et al. (2000). “Floodplains”, defined as the maximum areal proportion of a grid cell that can be flooded when the river exceeds its bankfull flow, and “Swamps”, representing groundwater-fed wetlands in the floodplain, were adopted from the Global Lake and Wetland database (Lehner and Döll, 2004). Depending on the areal extent of these swamps, a proportion of stream flow is simulated to feed into the soil moisture storage of the grid cell considered. Both parameters have an effect on the simulated river discharge and soil hydrology in the floodplains. For details, see Lauerwald et al. (2017).

Table 1 List of the forcing files used for our simulations, along with their spatiotemporal resolution.

Parametrization of hydrological processes

ORCHILEAK was previously parametrized and validated for the Amazon (Lauerwald et al., 2017, 2020; Hastie et al., 2019), Congo (Hastie et al., 2021) and Lena (Bowring et al., 2020) basins. In our study of the European river network, we updated ORCHILEAK with the more recent hydrology scheme of the recent standard version of ORCHIDEE (svn 5091). This hydrology scheme has been calibrated against observed runoff at a global scale (Ringeval et al., 2012; Yang et al., 2015). Furthermore, MacBean et al. (2020) has evaluated the model performance for simulating soil moisture in temperate ecosystems. This new hydrology scheme features a dynamic surface roughness of the vegetation, which decreases the aerodynamic resistance near the surface when vegetation cover is low, leading to lower ground temperatures and thus lower evaporation rates. This adjustment was deemed necessary in order to better capture the observed mean and seasonal variability of the discharge along the European river network. The two reduction factors controlling DOC leaching from the topsoil to the headwater streams were also adjusted (see Sect. 2.1.4).

Spin-up

Before the model can be used to simulate C dynamics over the past decades, a spin up is needed to reach an assumed steady state for the C fluxes during the pre-industrial period. This steady state is achieved by spinning up ORCHILEAK for 15 000 years. The spin up was realized by recursively looping over 4 years of climate forcing using the WFDEI forcing dataset over the 1979–1982 period (because the first year of the forcing, 1978, is incomplete) and constant land cover and atmospheric CO₂ concentration of 286 ppm (Guimberteau et al., 2018) corresponding to year 1861. After the end of the spin-up, the soil C stock across the entire European continent changed by less than 1 % over 1 century of simulation, which we considered close enough to steady state.

Transient runs

Using the steady-state outputs as initial condition, the first part of the transient simulation (1861–1978) was carried out with increasing atmospheric CO₂ concentration, changing land use and land cover and with river routing activated while still looping over 27 years (1979–2006) of the WFDEI forcing dataset. From 1979 onwards, the WFDEI atmospheric forcing data were applied over the entire period covered by this product with the changing land cover map and atmospheric CO₂ values applied for each year of simulation.

2.2.2 Model evaluation

Firstly, the simulated discharges were compared to times series of daily stream flow recorded at 11 gauging stations from “The Global Runoff Data Center (GRDC), 56068 Koblenz, Germany” dataset. For comparison, both the observed and simulated discharges were aggregated at the monthly temporal resolution over the years 1980 to 2006. Note that the river network in ORCHILEAK does not always match the real river network. The selected gauging stations were assigned to the cell best representing the parts of the river network to which the sampling location corresponds. However, important correction had to be performed for the most upstream stations in the Rhine and the Elbe river network. The period 1980 to 2006 was chosen based on the GRDC data coverage.

Model performance was further evaluated with respect to several variables of the terrestrial C cycle. Firstly, simulated net primary production (NPP) was compared to two different data products. The first one, the CARbon DAta MOdel fraMework (CARDAMOM; Bloom and Williams, 2015), was built from model data fusion analysis at 1^∘ resolution. The second one is the Global Inventory Modeling and Mapping Studies (GIMMS) at 0.5^∘ resolution based on AVHRR and MODIS sensors. GIMMS uses several atmospheric forcing data sets to derive NPP. These are CRUNCEP version 4 P1 and P2 (rainfall, cloudiness, relative humidity and temperature taken from the Climate Research Unit, CRU; the other fields, such as air pressure, longwave radiation and wind speed are directly derived from National Center for Atmospheric Research, NCEP), the European Centre for Medium-Range Weather Forecasts (ECMWF), the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2) and NCEPR2 (Reanalysis 2; https://www.esrl.noaa.gov/psd/, last access: 22 August 2019). For our comparison, we calculated the average of the NPP obtained with these five atmospheric forcing files. The NPP values from ORCHILEAK and GIMMS were averaged over the period 1982–2006 while CARDAMOM only covers a shorter time period, between 2001 and 2010. Modeled NPP was then compared to the NPP data products at the European scale and at the scale of five large European basins for which we also evaluated the simulated river discharge and DOC fluxes, which, taken together, represent 19 % of the model domain (Fig. 3): the Danube, Rhine, Elbe, Rhône and Seine. All five basins are located in an oceanic or humid continental climate (Fig. 2a), although the Rhône Basin extends further into the Mediterranean climate zone. The basin characteristics according to land cover types are as follows: the Danube and the Elbe basins have both a high proportion of croplands (around 40 %), the remainder being mostly covered by grasslands and boreal forests. Of the Rhône Basin, 50 % is covered by grasslands, while in the Seine Basin, croplands reach 50 %. The Rhine has a more diverse land cover, with a substantial proportion, about 30 %, comprised of boreal (10 %) and temperate (20 %) forests, 35 % of grasslands and 25 % of croplands. See Table S2 for further details.

Figure 3 Map of continental Europe delineating the (group of) catchments of focus in this study and the location of observed discharge and DOC concentrations. Catchments, from west to east, are All UK (light brown), Seine (orange), Rhône (yellow), Rhine (dark green), Elbe (violet), All Baltic (pink) and Danube (light green). Observations include the GRDC stations (red diamonds) in the Seine (S1 – Poses), Rhône (Ro1 – Beaucaire), Rhine (Ri1 – Lobith, Ri2 – Main in Frankfurt, Ri3 – Basel), Elbe (E1 – Neu Darchau, E2 – Dresden) and Danube (D1 – Ceatal Izmail, D2 – Svistov, D3 – Tisza in Senta, D4 – Bratislava) catchments and the river stations where DOC concentrations were measured (purple triangles): A1 – Douro, A2 – Sado, A3 – Gironde, A4 – Loire, A5 – Scheldt, A6 – Ems, A7 – Wales, A8 – Thames, M1 – Tech, M2 – Wales, M3 – Denmark and M4 – Finland (Abril et al., 2002; Mattsson et al., 2009).

The soil temperature is compared to the soil temperature generated using data from the ECMWF Reanalysis v5 (ERA5) dataset (Muñoz-Sabater et al., 2021). The soil heterotrophic respiration (SHR) is compared against the data-driven global SHR dataset published by Yao et al. (2021). The global SHR data set was produced using a random forest algorithm, upscaling from 455 data points from the global Soil Respiration Database (SRDB 4.0) based on gridded fields of climatic, edaphic and productivity-related variables as predictors (Yao et al., 2021). We compared the results of ORCHILEAK with the average, minimum and maximum values of SHR estimated by Yao et al. (2021). SOC from the Harmonized World Soil Database (HWSD) was used to evaluate the simulated SOC stocks. HWSD is a global soil database that contains up-to-date information on a large range of soil properties. For instance, this dataset reports the organic carbon content in the soil as well as the soil bulk density. The bulk density in HWSD was calculated in two different ways. The first one follows the method described in Saxton et al. (1986) where the bulk density is related to the soil texture – an approach tending to overestimate density in high porosity soils or in OC-rich soils. The second method uses the SOTWIS database in which the bulk density is calculated as a function of soil type and depth. In this database, all variables are reported for the topsoil (0–30 cm) and the sub-soil (30–100 cm) horizons. For comparison purposes, our simulated SOC stocks were thus integrated over the same depth intervals. We further assessed the extent to which our model can reproduce the main features in the observed soil DOC profiles. To that end, we compared our simulated DOC profile averaged over the entire European forest biome against the one established by Camino-Serrano et al. (2014) on the basis of a synthesis of local measurements. Although there are many studies on DOC concentrations in the soil, we selected that of Camino-Serrano et al. (2014) because it provides a synthesis at the pan-European scale, and is thus ideal for extracting “representative” concentration profiles over a sufficiently large domain, compatible with the regional scope of our study. Unfortunately, similar synthetic profiles based on observations have not been constructed for croplands and grasslands.

The key variables of interest in our study are the DOC leaching flux from the soil and the DOC export flux to the coast. These fluxes require accurate simulation of the water discharge fed by runoff and drainage, as well as of DOC concentrations in the leaching flux and in the riverine flux. For the leaching flux, our simulation results were compared to measured fluxes reported by Kindler et al. (2011) across different locations in Europe. Because the observed DOC leaching fluxes from both topsoil and bottom soil reported by Kindler et al. (2011) are based on local measurements that are not easily comparable to simulated fluxes at the coarse spatial resolution of our model (0.5^∘ or about 2×10³ km² at the corresponding latitude), we consider the comparison against measured river DOC fluxes to be more relevant for our purposes, as rivers are good integrators of mean, larger-scale catchment properties. For the riverine export fluxes, we assessed the modeled discharges and DOC concentrations separately. For evaluation of stream DOC concentrations, DOC data were extracted from the GLObal RIver CHemistry database (GLORICH, Hartmann et al., 2014) for the Rhine and Elbe basins and from the “Eau de France” database for the Seine and Rhône basins. These data were complemented by river DOC concentrations reported by Abril et al. (2002) for 9 river mouths (Sado, Thames, Ems, Scheldt, Gironde, Douro, Loire, Elbe and Rhine), and by Mattsson et al. (2009) for several river basins located in Finland, Denmark, Wales and France.

3 Results and discussion

3.1 Model evaluation at pan-European and catchment scales

3.1.1 Discharge

Figure 4 compares the simulated discharge against observations for selected stream gauging stations (Sect. 2). Those stations are located near the mouth of large rivers (Danube, Rhine, Rhône, Elbe and Seine) but also include a few locations further upstream the same rivers or at major tributaries (Fig. 3). The comparison is performed for the period 1990–2000, except for the Rhône at Beaucaire and the Danube at Svistov, for which the observed stream gauge data cover only a shorter period. Overall, the model reproduces the observations well, both in terms of amplitude and seasonality, except for the Elbe at Neu Darchau, for which the temporal variability is well captured but the absolute discharge is overestimated.

Figure 4 Modeled (black) and observed (red) time series of discharge at the GRDC gauging stations on the Danube (a–c) and its tributaries (d), the Elbe (e, f), the Rhine (h, i) and its tributaries (g) and the Rhône (j). Note the different time periods of measurements. See Fig. 3 for exact locations.

Note that the simulated catchment area often diverges (by −25 % to +30 %) from the observed value due to the coarse resolution ($\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$) of ORCHILEAK (Table S3). As a result of the model resolution, smaller tributaries are not represented individually and each grid cell was fully assigned to one larger river basin. The effect of the resolution is also shown in Fig. 2d, which compares the observed and modeled river network. Discrepancies between model and real world catchment area will translate into proportional biases in discharge simulation. Furthermore the 0.5^∘ resolution is too coarse to be able to perfectly represent the pathways of the river. Our model tends more often to underestimate the catchment area, while its yearly mean discharge is overestimated, except at the Beaucaire station along the Rhône River. The bias can be significant and cannot be explained by the model resolution alone.

To evaluate model performance for discharge, we used the Pearson's coefficient of determination (R²) and the Nash–Sutcliffe modeling efficiency (NSE, Nash and Sutcliffe, 1970). The R² only accounts for the correlation with regard to the temporal variability. With R² 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 discharges. The statistics confirm our previous observation that the model generally overestimates discharges (low NSEs) except for stations Elbe in Dresden, Rhône 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 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.

3.1.2 NPP, SHR, soil temperature, biomass and soil organic C stocks

We briefly compare simulated NPP with the gridded observation-based products GIMMS and CARDAMON (Sect. 2.2.2) 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 comprising 35 to 70^∘ N and 10^∘ W to 30^∘ E – 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 g C m⁻² yr⁻¹, a value in remarkable agreement with both GIMMS and CARDAMOM estimates of 430 and 460 g C m⁻² yr⁻¹, respectively. These two datasets entail an uncertainty that we assume to be 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 PgĊ or 2.3 kg C m⁻². This value is in good agreement with a recent estimate by Avitabile and Camia (2018), who report a biomass stock of around 16 Pg C. We estimate that the total soil carbon stock amounts to 58 Pg C. Averaged over the first meter of the soil horizon, this corresponds to a value of 9.5 kg C m⁻², which is comparable to that of HWSD (6 kg C m⁻²) when using the SOTWIS method to compute the bulk density, but significantly lower when applying Saxton's method (22 kg C m⁻²), plausibly because the latter overestimates the bulk density in OC-rich soils (Köchy et al., 2015). The results for the catchment scale are described in the Supplement.

Figure 5a shows the comparison of simulated and data-driven estimates of soil temperature. Soil temperature is overall well represented with a simulated mean temperature of 8.4 ^∘C against 9.3 ^∘C after ERA5. The overall slight underestimation is due to a substantial underestimation of soil temperatures in the northern regions. Furthermore, this underestimation is more important in the winter (January, February and March) with a difference relative to ERA5 values reaching 3.5 ^∘C, while in summer (July, August and September) this difference amounts to only 0.5 ^∘C. This underestimation is due to the poor representation of the isolating effect of snow cover in ORCHIDEE (Wang et al., 2013).

Figure 5 Difference in (a) simulated soil temperature (in ^∘C) against values reported by ERA5; (b) simulated soil heterotrophic respiration (in g C m⁻²,yr⁻¹) against values reported by Yao et al. (2021).

A comparison of simulated SHR against the Yao et al. (2021) estimates is shown in Fig. 5b. Over Europe, ORCHILEAK underestimates the SHR compared to the Yao et al. (2021) estimates by about 14 %, with a simulated average of 312 g C m⁻² yr⁻¹ against an average of 363 g C m⁻² yr⁻¹ (range from 317 to 417 g C m⁻² yr⁻¹) estimated by the latter for the period 1985–2013. Looking at specific climatic regions, some regions are well represented in ORCHILEAK, as the Mediterranean and humid continental regions show a mean SHR of 371 and 363 g C m⁻² yr⁻¹, against 385 and 354 g C m⁻² yr⁻¹ by Yao et al. (2021), respectively. SHR in semi-arid and tundra regions are conversely around 50 % lower than the Yao et al. (2021) estimate. For the tundra region, the underestimation in SHR is consistent with an underestimation of the NPP (see Table S4).

3.1.3 Soil DOC stocks

Comparison between observed and modeled DOC stocks and fluxes is more difficult than for biomass and SOC because the former have not been assessed at large spatial scales. Nevertheless, representative soil DOC concentration profiles for coniferous and broadleaved forests of Europe have been compiled by Camino-Serrano et al. (2014). These profiles were used to evaluate our model. Overall, we found that ORCHILEAK slightly overestimates DOC concentrations, especially in the very top soil layers with modeled values of around 100 mg L⁻¹ against 40–60 mg L⁻¹ in the observations (Fig. 6). We also simulated higher concentrations in broadleaved forests than in coniferous forests, while Camino-Serrano et al. (2014) obtained the opposite. When integrated over the first meter of the soil horizon of forested ecosystems (28 % of the European land area), the modeled and observed DOC stocks amount to 22.2 and 11.3 g C m⁻², respectively. Above, we have shown that over Europe, SOC stocks were underestimated, while the average DOC concentrations in the soil over all European forests are overestimated here. One explanation for the underestimation of SOC stocks and the likely overestimation of DOC stocks is thus that SOC decomposition rates in the new soil carbon module may be slightly too high. It is however difficult to generalize this conclusion because of the lack of synthesis data for other land cover types, especially croplands and grasslands, which together represent about 50 % of the total European land area. Modeled DOC stocks averaged over broad climate regions reveal the highest values for the oceanic climate, with 32 g C m⁻², and the Mediterranean climate, with 26 g C m⁻². Semi-arid and humid continental climates have similar stock densities of, respectively, 17.5 and 20 g C m⁻², and it is in the coldest climates (subarctic and tundra) that we find the lowest DOC stock densities, around 8 g C m⁻².

Figure 6 Modeled (blue) versus observed (red dashed) DOC concentration profiles averaged over the soils of the European coniferous and broadleaved forest biome. Data are from Camino-Serrano et al. (2014). The shaded area represents the 95 % bootstrap confidence interval for model and observations.

3.1.4 DOC leaching fluxes

The model simulates a yearly mean DOC leaching flux over Europe of 14.3 (±10) Tg C yr⁻¹ (Fig. 7), the standard deviation being here coarsely approximated by spatial variability. The average area specific flux rate amounts to 2.6 (±2.5) g C m⁻² yr⁻¹. 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 to March 2008. Comparing model results at 0.5^∘ resolution to point measurement is complicated, and thus in this section, we compare only our model-averaged result against the 17-site average from Kindler et al. (2011). Our modeled average of 2.6 (±2.5) g C m⁻² yr⁻¹ is of the same order of magnitude as that observed (4.2 g C m⁻² yr⁻¹). 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 model grid cells encapsulates the observational mean, highlighting a significant heterogeneity that is difficult to capture with local measurements alone. This comparison must be taken cautiously because of the limited number of observations and the resolution of our model. Furthermore, the DOC leaching flux at the coast is generally not well represented. The reason is that DOC leaching fluxes are normalized by the area of the whole cell. The area normalized flux at the coast is thus often lower.

Figure 7 Modeled yearly mean terrestrial DOC leaching flux (period 1979–2006) to the river European river network (in g C m⁻² yr⁻¹). The local observations from Kindler et al. (2011) are also reported, using the same scale. Note that the local observations cover a much shorter time period and may not be representative of the whole year.

The seasonal distribution of the DOC leaching flux is shown in Fig. 8. On average, the leaching flux per season averaged over Europe amounts to 1.6, 1.3, 0.5 and 1.4 Tg C per month in winter, spring, summer and autumn, respectively. If we exclude the high-latitude and high-altitude regions (Scandinavia, the Alps), a clear seasonality is observed, with the lowest fluxes in summer and spring and the highest fluxes in winter and autumn. In the high-latitude and high-altitude regions, the pattern is different, with the highest fluxes in spring, which extends to the summer in the Alps and corresponds to the snowmelt period. The highest fluxes per unit area are simulated in Scandinavia during the spring season, even though peatlands are not represented in the model. Some regions are leaching hotspots such as the Alps throughout the year, the western Balkans during autumn and the western flank of the UK in autumn and winter.

Figure 8 Seasonal distribution of the terrestrial DOC leaching flux (g C m⁻² per month) for (a) winter, (b) spring, (c) summer and (d) autumn, averaged over the period 1979–2006. A logarithmic scale is used to better highlight the spatiotemporal gradients.

3.1.5 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 and 0.01 d⁻¹ 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⁻¹, which corresponds to a half-life for riverine DOC of about 14 d (Table 2). This rate constant varies across Europe but always remains within the same order of magnitude, with half-lives ranging from 6 to 20 d (0.035–0.122 d⁻¹). These decomposition rates are in good agreement with the average rate reported by Berggren and Al-Kharusi (2020) of 0.037 d⁻¹, based on field experiments carried out at multiple river sampling locations across Europe. We thus conclude that DOC decomposition rates used in ORCHILEAK are reasonable, and fluvial DOC fluxes are a valid proxy to evaluate simulated DOC leaching fluxes.

Table 2 Estimated river DOC decay rates applied in ORCHILEAK. Values are reported for four large river basin and for the six dominant climate zones.

Figure 9 compares modeled versus observed multi-annual mean riverine DOC concentrations at specific locations or within a group of small river catchments. Local DOC measurements include data near the mouth of the Rhine, Elbe, Rhône and Seine rivers (discharge, DOC concentration and fluxes for the Rhine and Seine in Fig. 10). In addition, Abril et al. (2002) report DOC 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) and 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. The measured values are equal to 4.3 and 2.9 mg C L⁻¹ for the Rhine and 4.6 and 6.1 mg C L⁻¹ for the Elbe, highlighting the inherent variability in observational data. To complement these local samplings, we also compared our simulated DOC concentrations with those of Mattsson et al. (2009) 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 to 10 mg C L⁻¹, except in two regions in the north (Finland and basins flowing into the Baltic sea), where concentrations exceeded 10 mg C L⁻¹. For most of the data, the model slightly overestimated the river DOC concentrations. The model results also suggest that DOC 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 Mattsson et al. (2009). Such a pattern possibly results from the decreasing mean annual air temperature and runoff in northern Europe that favor incomplete decomposition of litter and soil DOC, thus favoring DOC production in the soil, while at the same time, DOC turnover rates in the soils decreased. The increased abundance of forests, in particular coniferous forests, is also a valid explanation for higher DOC leaching (Lauerwald et al., 2012). However, it is important to keep in mind that we are not representing peatlands, suggesting that we could lack part of the DOC leaching in subarctic and tundra regions, leading to even higher DOC fluxes further north. Finally, the comparison reveals that model performance tends to improve with catchment size, likely reflecting the difficulty of capturing the DOC dynamics at a small scale with the current resolution of ORCHILEAK. However, overall, our model is capable of reproducing observed yearly mean DOC concentrations for a wide range of river basins spread between Finland and Portugal.

Figure 9 Modeled river DOC concentration against observed values. The color code indicates the dominant climate zone for each catchment while the size of the diamond is proportional to the catchment area according to the following classes: <10 000, <50 000, <100 000 and >100 000 km². See Table S4 for further details.

Figure 10 Time series of discharge (a, d), DOC concentration (b, e) and DOC fluxes (c, f) in the Rhine at Lobith (a–c, period 1992–1996) and in the Seine at Poses (d–f, period 2002–2006. See Fig. 3 for location of stations. Black lines are for ORCHILEAK and red dots for observations.

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 (Rhine and Seine illustrated in Fig. 10). In terms of interannual variability (IAV), riverine DOC fluxes present the highest variability with a coefficient of variance (COV) of 0.62 for the Seine and 0.57 for the Rhine. For comparison, IAV of discharge and riverine DOC concentration shows COVs of 0.60 and 0.51 for the Seine and 0.40 and 0.45 for the Rhine, respectively. The higher IAV for the flux is due to the tendency for higher concentrations to coincide with higher discharge, which is due to the flushing effect, where higher discharges follow higher amounts of runoff from topsoils rich in DOC. The multi-year mean modeled DOC fluxes are estimated for the Rhine, Elbe, Rhône and Seine at 11.9, 7.2, 8.8 and 3.2 kg s⁻¹, respectively. The observations amount, respectively, to 7.9, 3.6, 4.6, 1.6 kg s⁻¹. For all stations, the model thus slightly overestimates 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 Rhône compared to that observed yields an R² of 0.6 and a mean error of 92 % (results for the Seine, Elbe and Rhine are reported in Table S6). In Fig. 9, we clearly see a large overestimation of the temporal variability in DOC concentrations and this could lead to an overestimation in DOC fluxes since there is a positive relationship between concentrations and discharge. The overestimation of DOC concentrations and consequently of DOC fluxes could be due to high DOC leaching.

The overestimation of DOC fluxes can also be due to the fact we choose to not recalibrate the hydrology scheme but instead we optimized the model for the discharge by adjusting the surface roughness of the vegetation (Sect. 2.2.1). Since those four stations are all located in the same region with the same type of land cover (western Europe), two other locations have been selected: the UK and the Baltic Sea. For those two locations, there are no time series data for DOC flux but some studies have measured DOC concentrations/fluxes. Worrall et al. (2012) estimated DOC concentration across the UK and Fransner et al. (2016) reported modeled DOC concentrations for all the catchments flowing into the Baltic Sea.

Finally, although the model–data comparison points to a slight overestimation of the river DOC export flux, our pan-European estimate amounts to 12.3 Tg C yr⁻¹. This estimate is in fact about 35 % lower than that reported in another model study (Li et al., 2019) based on TRIPLEX-HYDRA, a process-based model for which the DOC export flux reaches 19.3 Tg C yr⁻¹. Li et al. (2019) applied the model at the global scale, and simulation results were primarily evaluated against observations for the world's largest rivers and here 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, the export fluxes predicted by both models fall within the same order of magnitude.

Figure 11 Comparison of modeled yearly mean terrestrial DOC leaching flux (period 1979–2006) into the European river network (in g C m⁻² yr⁻¹), with (a) and without (b) the representation of manure application. (c) Increase of DOC leaching in percentage compared to DOC leaching without manure implementation. (d) Comparison of total DOC leaching (solid line) and DOC leaching through runoff only (dashed line) throughout Europe with and without manure application (period 2004–2006).

3.1.6 Manure implementation

The implementation of manure significantly affects DOC leaching from grasslands and croplands (Fig. 11), which cover more than half of the studied region. The average annual input rate of manure into the soil is around 2.5 g C m⁻² yr⁻¹ (Fig. 2c). With manure implementation, the DOC leaching rate increases drastically (an average of +72 % compared to 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 g C m⁻² yr⁻¹ and from 1.7 to 2.5 g C m⁻² yr⁻¹, respectively. For the whole of Europe, manure implementation leads to an increase of total DOC leaching into the river network from 9.8 to 14.3 Tg C yr⁻¹ (Fig. 11a and b). Figure 11d shows that the application of manure increases DOC leaching in particular in winter (January, February and March) while in summer (July, August and September) the increase is relatively low. In ORCHILEAK, the manure-derived DOC first enters the topsoil. There, a part of it is decomposed, and the rest is transported to deeper soil layers with percolating water. Finally, a variable part of the DOC derived from manure is flushed out of the soil column with the surface runoff and belowground drainage. As the manure first enters the topsoil, one could expect that it would mainly increase the DOC leaching from the topsoil. However, our results show that the application of manure does not influence the ratio of DOC leaching through surface runoff to belowground drainage. Over Europe, the average increase in topsoil DOC leaching due to manure is equal to half of the total increase, the rest of the increase being contributed by the drainage.

3.2 European-scale DOC leaching dynamics

3.2.1 Drivers of DOC leaching

Here, we analyze what controls the spatial distribution and temporal variability in DOC leaching. While the grid cell and the basin scales were the most relevant for the model evaluation, when searching for potential drivers of soil DOC leaching, such as temperature, runoff and drainage (driven by precipitation), a climatologic segmentation of the European domain can help to better explain the impact of these drivers. Figure 12 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.

Figure 12 Simulated DOC leaching flux (g C m⁻² per month; blue) and total runoff (mm per month; red) for the six largest climate zones in Europe (period 2004–2006).

To further explore the environmental controls of DOC leaching, we calculate the partial correlation factor for surface runoff, drainage, temperature and NPP, and they are, respectively, 0.43, 0.54, −0.17 and 0.18, highlighting that surface runoff and drainage can explain most of the spatiotemporal variability in DOC leaching fluxes, with temperature and NPP only playing a subordinate role. We decided to express DOC leaching as fraction of the annual terrestrial NPP (Fig. 13). In doing so, we assume that NPP, which is undoubtedly the ultimate C source for DOC production (since litter and SOC stocks, the sources of DOC, are fed by NPP), also acts as an important control of the DOC leaching flux. Moreover, in 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 13 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 drainage (Fig. 13b). This can be explained by the general decrease in soil DOC concentrations with depth (Fig. 6), 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 disproportionately high amount to the long-term DOC leaching (Idir et al., 1999; Raymond and Saiers, 2010). Finally, we found higher ratios of leaching to NPP at lower temperatures (Fig. 13a), suggesting that lower temperatures lead to longer turnover times of DOC in the soil, and thus higher concentrations in the leaching flux (Sect. 2.1.2).

Figure 13 Fraction (%) of terrestrial NPP that is leached as DOC into the river network as a function of total runoff. Each point represents the grid-cell average of both metrics for the entire simulation period (1979–2006). In panel (a) the color scale represents the grid-cell average temperature (^∘C) while in panel (b), the color scale represents the ratio of surface runoff to total runoff in percentage. Panel (c) shows the normalized predicted ratio of DOC leaching flux to NPP (Eq. 13) against the normalized simulated values.

Table 3 Key physical and biogeochemical characteristics of the six dominant climate zones of the European domain.

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). The idea behind this rationale is to highlight that once normalized to the terrestrial NPP, we can directly analyze which physical drivers impact the fraction of NPP that is lost to leaching from terrestrial ecosystems. To compare the importance of each predictor for the spatiotemporal patterns of DOC leaching, we normalized all variables, V_s, of Eq. (13) according to Eq. (14), where i is the cell index.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{DOC}\mathrm{leaching}}{\mathrm{NPP}}}& =K_{\mathrm{0}}+K_{\mathrm{R}}\cdot \mathrm{Runoff}+K_{\mathrm{D}}\cdot \mathrm{Drainage}\\ & +K_{\mathrm{T}}\cdot e^{T\left({}^{\circ }\mathrm{C}\right)}\end{array}\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{K}_{\mathrm{0}}=\mathrm{0.01}\pm \mathrm{7}\ast \cdot {\mathrm{10}}^{-\mathrm{4}}\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{R}}=\mathrm{0.342}\pm \mathrm{0.009}\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{K}_{\mathrm{D}}=\mathrm{0.276}\pm \mathrm{0.014}\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{T}}=-\mathrm{0.055}\pm \mathrm{0.013}\end{array}$$

(p value < $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{16}}$ except for temperature, where p value = $\mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{5}}$)

$$\begin{array}{}\text{(13)}& V_{i,\mathrm{N}}={\displaystyle \frac{V_i-V_{\min }}{V_{\max }-V_{\min }}}\end{array}$$

To rule out any significant multi-collinearity in the regression model, for each predictor we calculated 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, the VIFs of the runoff, drainage and temperature are, respectively, 1.13, 1.13 and 1.01, confirming that our prediction is robust and not biased by high multicollinearity.

In Fig. 13c, the DOC leaching simulated by ORCHILEAK is compared with the one predicted by Eq. (13). Our simple regression model is able to reproduce the simulations with a residual standard error of 0.68 % and an R² of 0.45. The coefficients of our regression model reveal that spatiotemporal variability in DOC leaching is mainly driven by surface runoff (K_R) and drainage (K_D). Air temperature, as a third control of DOC leaching, is of subordinate importance, as reflected by its low predictor coefficient (K_T).

Table 3 summarizes, for each climate zone in Europe, the DOC leaching fluxes (in absolute value and normalized by the 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 yr⁻¹) and temperatures are high.

The highest fraction of NPP exported to rivers as DOC is found in the tundra climate, reaching 1.22 %. This can be explained by high runoff and drainage (reaching 920 mm yr⁻¹) 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 ratio of DOC leaching to NPP, with a value of 0.84 %. The Mediterranean, oceanic and humid continental climate zones present ratios of intermediate DOC leaching to NPP 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 %.

3.2.2 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₂ 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 g C m⁻² yr⁻¹. Our results, 2.6±2.5 g C m⁻² yr⁻¹, support this hypothesis although we 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. Uncertainties in the processes included or omitted in the model could also explain some of the discrepancy. In terms of temporal variability, we found the highest DOC leaching in winter averaged over the continent (8.9 Tg C in total for the 6 months of winter October to March) and the lowest in summer (5.4 Tg C over the period April to September), consistent with the findings of Kindler et al. (2011). In terms of 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. (2019, 2021) found that, respectively, 12 % and 4 % of the NPP are exported each year to inland waters in the form of DOC and CO₂ – much higher than we report for Europe (0.6 %), but this value only accounts for DOC. Note that for these tropical lowland river basins, extensive riparian wetlands are an important source of DOC, which are of minor importance in Europe. For the Lena River basin located in the boreal region, Bowring et al. (2020) found a ratio of DOC leaching to NPP 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 USA can be calculated when dividing the average DOC leaching flux of 2.7 g C m⁻² yr⁻¹ simulated by Tian et al. (2015) by the average NPP of 780 g C m⁻² yr⁻¹ estimated by Zhao et al. (2005). Further, our value is quite similar to the one extracted from the global study by Nakhavali et al. (2021), which 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 wetlands, which have been reduced by major regulation of the European river network.

3.2.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 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 was 860 Tg C yr⁻¹ (123 g C m⁻² yr⁻¹), equalling about 28 % of the total NPP (Fig. 14b). The ratio of DOC leaching to NEE shows drastic spatial variation, varying from an average value of 0.4 % in the semi-arid region to a value of 5.7 % in the tundra region. For the whole Europe, the DOC leaching is about 3 % of the NEE.

Figure 14 Grid-cell average of (a) net primary production (NPP), (b) net ecosystem exchange (NEE) and (c) fraction of NPP leached into the river network as DOC (%) for the period 1979–2006.

3.3 Model limitations

ORCHILEAK is an LSM that simulates the impact of riverine DOC transfers on the terrestrial C budget, but it still suffers from several limitations. In fact, ORCHILEAK cannot represent all biogeochemical transformation processes affecting DOC in the soil column and the river network (Lauerwald et al., 2017). For instance, environmental controls such as soil pH and ionic strength have been demonstrated to have an impact on DOC solubility in soils (Monteith et al., 2007) and thus affect DOC leaching into streams. Unfortunately, these processes are not represented in our model, as there are still no reliable methods and forcing data to simulate the dynamics of soil pH and ionic strength in the soil solution at a large scale.

As mentioned before, peatlands are not included in the model, yet they cover a large part of northern Europe. Peatlands are known to play an important role in the C cycle, and are an important source of DOC for the river network. One of the major next steps would thus be to merge ORCHILEAK with ORCHIDEE-PEAT, a new branch of the LSM ORCHIDEE that simulates the development and C balance of peatlands (Qiu et al., 2019).

Another source of DOC originates from wastewater treatments plants, which are not included in the model due to the lack of forcing data related to sewage water treatment. It has been shown that DOC concentrations in sewage are important (Griffith et al., 2009). However, Meybeck (1986) showed that DOC from sewage is very labile and only affects concentrations within short distances downstream of water processing plants. Having avoided observational data from sites known to be impacted by sewage effluents directly, we assume that our model–data evaluation was not impacted by this potential DOC source. To assess the role of soil DOC leaching in the terrestrial C budget, sewage is not a contribution of direct interest.

While riparian zones are a major source of DOC for the river network (Inamdar and Mitchell, 2006; Grabs et al., 2012), the impact of riparian zones on DOC leaching through runoff into the river network is only implicitly represented in the model (as described in Lauerwald et al., 2017). Due to the coarse resolution of the model, riparian zones around small streams (order 1 to 3) cannot be explicitly included in the model. It is assumed that the extent of the riparian zones, from which most of the DOC stems, scales linearly to the surface area of these small streams, both in time as well as in space (i.e., between different grid cells of our model grid). While the surface area of these small streams is not directly represented, Lauerwald et al. (2017) assumed that the spatial and temporal variations in this stream surface area scale to the square root of the discharge that is flowing through these streams, roughly in line with empirical scaling laws (e.g., Raymond et al., 2012). For the larger rivers, for which the surface area is explicitly represented in the model, it is assumed that the inundated riparian zone can temporally comprise up to 10 % of the river water surface area, depending on the temporal variability of discharge. Here, DOC produced from decomposition of litter and SOC is directly injected into the river water body.

ORCHILEAK could further be improved through the implementation of lakes and reservoirs. It has been shown that dams have a direct impact on C retention efficiency in the inland water river network (Maavara et al., 2017). So far, ORCHILEAK does not represent lateral transport of POC at the continental scale, yet its non-negligible role in the terrestrial C budget has been demonstrated (Zhang et al., 2018, 2020; Naipal et al., 2018). Finally, the effect of nutrient limitation on the C cycle is not yet accounted for in ORCHILEAK. It has been demonstrated that the implementation of nitrogen (N) and phosphorus (P) could reduce the simulated land C sink (Goll et al., 2012; Sun et al., 2021). It can be assumed that nutrient limitation would similarly affect DOC leaching, and it could dampen its increase with the rising atmospheric CO₂ levels predicted by previous studies with ORCHILEAK (Lauerwald et al., 2020; Hastie et al., 2021).

4 Conclusion

We reconstructed the terrestrial and riverine C fluxes in Europe during the period 1979–2012 using the ORCHILEAK LSM. The total C leaching from soils into European rivers is 14.3 Tg C yr⁻¹ on average, about 0.6 % of the estimated NPP and 3 % of the terrestrial net uptake 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. For the whole of Europe, the DOC leaching rate is the highest in winter and lowest during the summer, mainly controlled by the seasonal variation of runoff. The implementation of manure leads 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 USA, an emergent view regarding the global role of DOC leaching in the terrestrial C balance, and its underlying drivers, is progressively emerging.