Simulating the effect of tillage practices with the global ecosystem model LPJmL (version 5.0-tillage)

The effects of tillage on soil properties, crop productivity, and global greenhouse gas emissions have been discussed in the last decades. Global ecosystem models have limited capacity to simulate the various effects of tillage. With respect to the decomposition of soil organic matter, they either assume a constant increase due to tillage or they ignore the effects of tillage. Hence, they do not allow for analysing the effects of tillage and cannot evaluate, for example, reduced tillage or no tillage (referred to here as “no-till”) practises as mitigation practices for climate change. In this paper, we describe the implementation of tillage-related practices in the global ecosystem model LPJmL. The extended model is evaluated against reported differences between tillage and no-till management on several soil properties. To this end, simulation results are compared with published meta-analyses on tillage effects. In general, the model is able to reproduce observed tillage effects on global, as well as regional, patterns of carbon and water fluxes. However, modelled N fluxes deviate from the literature values and need further study. The addition of the tillage module to LPJmL5 opens up opportunities to assess the impact of agricultural soil management practices under different scenarios with implications for agricultural productivity, carbon sequestration, greenhouse gas emissions, and other environmental indicators.

1 Introduction

Agricultural fields are tilled for various purposes, including seedbed preparation, incorporation of residues and fertilizers, water management, and weed control. Tillage effects a variety of biophysical processes that affect the environment, such as greenhouse gas emissions or soil carbon sequestration and can influence various forms of soil degradation (e.g. wind, water, and tillage erosion) (Armand et al., 2009; Govers et al., 1994; Holland, 2004). Reduced tillage or no tillage (hereafter referred to as “no-till”) is being promoted as a strategy to mitigate greenhouse gas (GHG) emissions in the agricultural sector (Six et al., 2004; Smith et al., 2008). However, there is an ongoing long-lasting debate about tillage and no-till effects on soil organic carbon (SOC) and GHG emissions (e.g. Lugato et al., 2018). In general, reduced tillage and no-till tend to increase SOC storage through a reduced decomposition and consequently reduces GHG emissions (Chen et al., 2009; Willekens et al., 2014). However, discrepancies exist on the effectiveness of reduced tillage or no-till on GHG emissions. For instance, Abdalla et al. (2016) found in a meta-analyses that on average no-till systems reduce CO₂ emissions by 21 % compared to conventional tillage, whereas Oorts et al. (2007) found that CO₂ emissions from no-till systems increased by 13 % compared to conventional tillage, and Aslam et al. (2000) found only minor differences in CO₂ emissions. These discrepancies are not surprising as tillage effects a complex set of biophysical factors, such as soil moisture and soil temperature (Snyder et al., 2009), which drive several soil processes, including the carbon and nitrogen dynamics and crop performance. Moreover, other factors such as management practices (e.g. fertilizer application and residue management) and climatic conditions have been shown to be important confounding factors (Abdalla et al., 2016; Oorts et al., 2007; van Kessel et al., 2013). For instance, Oorts et al. (2007) attributed the higher CO₂ emissions under no-till to higher soil moisture and decomposition of crop litter on top of the soil. Van Kessel et al. (2013) found that N₂O emissions were smaller under no-till in dry climates and that the depth of fertilizer application was important. Finally, Abdalla et al. (2016) found that no-till effects on CO₂ emissions are most effective in dryland soils.

In order to upscale this complexity and to study the role of tillage for global biogeochemical cycles, crop performance, and mitigation practices, the effects of tillage on soil properties need to be represented in global ecosystem models. Although tillage is already implemented in other ecosystem models at different levels of complexity (Lutz et al., 2019; Maharjan et al., 2018), tillage practices are currently underrepresented in global ecosystem models that are used for biogeochemical assessments. In these, the effects of tillage are either ignored or represented by a simple scaling factor of decomposition rates. Global ecosystem models that ignore the effects of tillage include, for example, JULES (Best et al., 2011; Clark et al., 2011), the Community Land Model (Levis et al., 2014; Oleson et al., 2010) PROMET (Mauser and Bach, 2009), and the Dynamic Land Ecosystem Model (DLEM; Tian et al., 2010). The models in which the effects of tillage are represented as an increase in decomposition include LPJ-GUESS (Olin et al., 2015; Pugh et al., 2015) and ORCHIDEE-STICS (Ciais et al., 2011).

The objective of this paper is to (1) extend the Lund Potsdam Jena managed Land (LPJmL5) model (von Bloh et al., 2018) so that the effects of tillage on biophysical processes and global biogeochemistry can be represented and studied and (2) evaluate the extended model against data reported in meta-analyses by using a set of stylized management scenarios. This extended model version allows for quantifying the effects of different tillage practices on biogeochemical cycles, crop performance, and for assessing questions related to agricultural mitigation practices. Despite uncertainties in the formalization and parameterization of processes, the process-based representation allows for enhancing our understanding of the complex response patterns as individual effects and feedbacks can be isolated or disabled to understand their importance. To our knowledge, some crop models that have been used at the global scale, e.g. EPIC (Williams et al., 1983) and DSSAT (White et al., 2010), have similarly detailed representations of tillage practices but models used to study the global biogeochemistry (Friend et al., 2014) have no or only very coarse representations of tillage effects.

2 Tillage effects on soil processes

Tillage affects different soil properties and soil processes, resulting in a complex system with various feedbacks on processes related to soil water, temperature, carbon (C), and nitrogen (N) (Fig. 1). The effect of tillage has to be implemented and analysed in conjunction with residue management as these management practices are often interrelated (Guérif et al., 2001; Strudley et al., 2008). The processes that were implemented into the model were chosen based on the importance of the process and its compatibility with the implementation of other processes within the model. Those processes are visualized in Fig. 1 with solid lines; processes that have been ignored in this implementation are visualized with dotted lines. To illustrate the complexity, we here describe selected processes in the model affected by tillage and residue management, using the numbered lines in Fig. 1.

Figure 1 Flow chart diagram of feedback processes caused by tillage, which are considered (solid lines) and not considered (dashed lines) in this implementation in LPJmL5.0-tillage. Blue lines highlight positive feedbacks, red negative, and black are ambiguous feedbacks. The numbers in the figure indicate the processes described in Sect. 2.

With tillage, surface litter is incorporated into the soil (1) and increases the soil organic matter (SOM) content of the tilled soil layer (2) (Guérif et al., 2001; White et al., 2010), while tillage also decreases the bulk density of this layer (3) (Green et al., 2003). An increase in SOM positively affects the porosity (4) and therefore the soil water holding capacity (whc) (5) (Minasny and McBratney, 2018). Tillage also affects the whc by increasing porosity (6) (Glab and Kulig, 2008). A change in whc affects several water-related processes through soil moisture (7). For instance, changes in soil moisture influence lateral runoff (8) and leaching (9) and affect infiltration. A wet (saturated) soil, for example, decreases infiltration (10), while infiltration can be enhanced if the soil is dry (Brady and Weil, 2008). Soil moisture affects primary production as it determines the amount of water which is available for the plants (11) and changes in plant productivity again determine the amount of residues left at the soil surface or to be incorporated into the soil (1) (feedback not shown).

The presence of crop residues on top of the soil (referred to as “surface litter” hereafter) enhances water infiltration into the soil (12) (Guérif et al., 2001; Jägermeyr et al., 2016; Ranaivoson et al., 2017), and thus increases soil moisture (13). That is because surface litter limits soil crusting, can constitute preferential pathways for water fluxes, and slows lateral water fluxes at the soil surface so that water has more time to infiltrate (Glab and Kulig, 2008). Consequently, surface litter reduces surface runoff (14) (Ranaivoson et al., 2017). Surface litter also intercepts part of the rainfall (15), reducing the amount of water reaching the soil surface, but also lowers soil evaporation (16) and thus reduces unproductive water losses to the atmosphere (Lal, 2008; Ranaivoson et al., 2017). Surface litter also reduces the amplitude of variations in soil temperature (17) (Enrique et al., 1999; Steinbach and Alvarez, 2006). The soil temperature is strongly related to soil moisture (18), through the heat capacity of the soil, i.e. a relatively wet soil heats up much slower than a relatively dry soil (Hillel, 2004). The rate of SOM mineralization is influenced by changes in soil moisture (19) and soil temperature (20) (Brady and Weil, 2008). The rate of mineralization affects the amount of CO₂ emitted from soils (21) and the inorganic N content of the soil. Inorganic N can then be taken up by plants (22), be lost as gaseous N (23), or transformed into other forms of N. The processes of nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) leaching, nitrification, denitrification, mineralization of SOM, and immobilization of mineral N forms are explicitly represented in the model (von Bloh et al., 2018). The degree to which soil properties and processes are affected by tillage mainly depends on the tillage intensity, which is a combination of tillage efficiency and mixing efficiency (explained in detail in Sect. 3.2 and 3.5.2). Tillage has a direct effect on the bulk density of the tilled soil layer. The type of tillage determines the mixing efficiency, which affects the amount of incorporating residues into the soil. Over time, soil properties reconsolidate after tillage, eventually returning to pre-tillage states. The speed of reconsolidation depends on soil texture and the kinetic energy of precipitation (Horton et al., 2016).

This implementation mainly focuses on two processes directly affected by tillage: (1) the incorporation of surface litter associated with tillage management and subsequent effects (Fig. 1, path 1 and following paths) and (2) the decrease in bulk density and the subsequent effects of changed soil water properties (Fig. 1, e.g. path 3 and following paths). In order to limit model complexity and associated uncertainty, tillage effects that are not directly compatible with the original model structure (such as subsoil compaction) or require very high spatial resolution are not taken into account in this initial tillage implementation, despite acknowledging that these processes can be important.

3 Implementation of tillage routines into LPJmL

3.1 LPJmL model description

The tillage implementation described in this paper was introduced into the dynamical global vegetation, hydrology, and crop-growth model LPJmL. This model was recently extended to also cover the terrestrial N cycle, accounting for N dynamics in soils and plants and N limitation of plant growth (LPJmL5; von Bloh et al., 2018). Previous comprehensive model descriptions and developments are described by Schaphoff et al. (2018a). The LPJmL model simulates the C, N, and water cycles by explicitly representing biophysical processes in plants (e.g. photosynthesis) and soils (e.g. mineralization of N and C). The water cycle is represented by the processes of rain water interception, soil and lake evaporation, plant transpiration, soil infiltration, lateral and surface runoff, percolation, seepage, routing of discharge through rivers, storage in dams and reservoirs, and water extraction for irrigation and other consumptive uses.

In LPJmL5, all organic matter pools (vegetation, litter, and soil) are represented as C pools and the corresponding N pools with variable C:N ratios. Carbon, water, and N pools in vegetation and soils are updated daily as the result of computed processes (e.g. photosynthesis, autotrophic respiration, growth, transpiration, evaporation, infiltration, percolation, mineralization, nitrification, and leaching; see von Bloh et al., 2018, for the full description). Litter pools are represented by the aboveground pool (e.g. crop residues, such as leaves and stubble) and the belowground pool (roots). The litter pools are subject to decomposition, after which the humified products are transferred to the two SOM pools that have different decomposition rates (Fig. S1a in the Supplement). The fraction of litter which is harvested from the field can range between almost fully harvested or not harvested, when all litter is left on the field (90 %; Bondeau et al., 2007). In the soil, pools of inorganic, reactive N forms (${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$) are also considered. Each organic soil pool consists of C and N pools and the resulting C:N ratios are flexible. Soil C:N ratios are considerably smaller than those of plants as immobilization by microorganisms concentrates N in SOM. In LPJmL, a soil C:N ratio of 15 is targeted by immobilization for all soil types (von Bloh et al., 2018). The SOM pools in the soil consist of a fast pool with a turnover time of 30 years, and a slow pool with a 1000-year turnover time (Schaphoff et al., 2018a). Soils in LPJmL5 are represented by five hydrologically active layers, each with a distinct layer thickness. The first soil layer, which is mostly affected by tillage, is 0.2 m thick. The following soil layers are 0.3, 0.5, 1.0, and 1.0 m thick, followed by a 10.0 m bedrock layer, which serves as a heat reservoir in the computation of soil temperatures (Schaphoff et al., 2013).

LPJmL5 has been evaluated extensively and demonstrated good skills in reproducing C, water, and N fluxes in both agricultural and natural vegetation on various scales (von Bloh et al., 2018; Schaphoff et al., 2018b).

3.2 Litter pools and decomposition

In order to address the residue management effects of tillage, the original aboveground litter pool is now separated into an incorporated litter pool (C_litter,inc) and a surface litter pool (C_litter,surf) for carbon, and the corresponding pools (N_litter,inc and N_litter,surf) for nitrogen (Fig. S1b in the Supplement). Crop residues not collected from the field are transferred to the surface litter pools. A fraction of residues from the surface litter pool are then partially or fully transferred to the incorporated litter pools, depending on the tillage practice:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{C}}_{\text{litter,inc},t+\mathrm{1}}={\mathrm{C}}_{\text{litter,inc},t}+{\mathrm{C}}_{\text{litter,surf},t}\cdot \text{TL}\text{ for carbon and}\\ \text{(1)}& {\displaystyle }& {\displaystyle }{\mathrm{N}}_{\text{litter,inc},t+\mathrm{1}}={\mathrm{N}}_{\text{litter,inc},t}+{\mathrm{N}}_{\text{litter,surf},t}\cdot \text{TL}\text{for nitrogen}.\end{array}$$

The C_litter,surf and N_litter,surf pools are reduced accordingly:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{C}}_{\text{litter,surf},t+\mathrm{1}}={\mathrm{C}}_{\text{litter,surf},t}\cdot (\mathrm{1}-\text{TL}),\\ \text{(2)}& {\displaystyle }& {\displaystyle }{\mathrm{N}}_{\text{litter,surf,}t+\mathrm{1}}={\mathrm{N}}_{\text{litter,surf},t}\cdot (\mathrm{1}-\text{TL}),\end{array}$$

where C_litter,inc and N_litter,inc are the amounts of incorporated surface litter C and N, respectively, in grammes per square metre (g m⁻²) at a time step t (days). The parameter TL is the tillage efficiency, which determines the fraction of residues that is incorporated by tillage (0–1). To account for the vertical displacement of litter through bioturbation under natural vegetation and under no-till conditions, we assume that 0.1897 % of the surface litter pool is transferred to the incorporated litter pool per day (equivalent to an annual bioturbation rate of 50 %).

The litter pools are subject to decomposition. The decomposition of litter depends on the temperature and moisture of its surroundings. The decomposition of the incorporated litter pools depends on soil moisture and temperature of the first soil layer (as described by von Bloh et al., 2018), whereas the decomposition of the surface litter pools depends on the litter's moisture and temperature, which are approximated by the model. The decomposition rate of litter (rdecom in gC m⁻² d⁻¹) is described by first-order kinetics, and is specific for each plant functional type (PFT) following Sitch et al. (2003);

$$\begin{array}{}\text{(3)}& {\displaystyle }{\text{rdecom}}_{\text{(PFT)}}=\mathrm{1}-\exp \left(-{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{10}\left(\text{PFT}\right)}}}\cdot g\left(T_{\mathrm{surf}}\right)\cdot F\left(\mathit{\theta }\right)\right),\end{array}$$

where τ₁₀ is the mean residence time for litter and F(θ) and g(T_surf) are response functions of the decay rate to litter moisture and litter temperature (T_surf), respectively. The response function to litter moisture F(θ) is defined as:

$$\begin{array}{}\text{(4)}& F\left(\mathit{\theta }\right)=\mathrm{0.0402}-\mathrm{5.005}\cdot {\mathit{\theta }}^{\mathrm{3}}+\mathrm{4.269}\cdot {\mathit{\theta }}^{\mathrm{2}}+\mathrm{0.7189}\cdot \mathit{\theta },\end{array}$$

where θ is the volume fraction of litter moisture which depends on the water holding capacity of the surface litter (whc_surf), the fraction of surface covered by litter (f_surf), the amount of water intercepted by the surface litter (I_surf) (Sect. 3.3.1), and lost through evaporation E_surf (Sect. 3.3.3).

The temperature function g(T_surf) describes the influence of temperature of surface litter on decomposition (von Bloh et al., 2018):

$$\begin{array}{}\text{(5)}& {\displaystyle }g\left(T_{\mathrm{surf}}\right)=\exp \left(\mathrm{308.56}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{66.02}}}-{\displaystyle \frac{\mathrm{1}}{\left(T_{\mathrm{surf}+\mathrm{56.02}}\right)}}\right),\end{array}$$

where T_surf is the temperature of surface litter (Sect. 3.4).

A fixed fraction (70 %) of the decomposed C_litter,surf is mineralized, i.e. emitted as CO₂, whereas the remaining humified C is transferred to the soil C pools, where it is then subject to the soil decomposition rules as described by von Bloh et al. (2018) and Schaphoff et al. (2018a). The mineralized N (also 70 % of the decomposed litter) is added to the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ pool of the first soil layer, where it is subjected to further transformations (von Bloh et al., 2018), whereas the humified organic N (30 % of the decomposed litter) is allocated to the different organic soil N pools in the same shares as the humified C. In order to maintain the desired C:N ratio of 15 within the soil (von Bloh et al., 2018), the mineralized N is subject to microbial immobilization, i.e. the transformation of mineral N to organic N directly reverting some of the N mineralization in the soil.

The presence of surface litter influences the soil water fluxes and soil temperature of the soil (see Sect. 3.3 and 3.4), and therefore affects the decomposition of the soil carbon and nitrogen pools, including the transformations of mineral N forms. Nitrogen fluxes such as N₂O from nitrification and denitrification, for instance, are partly driven by soil moisture (von Bloh et al., 2018):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{F}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{nitrification},l}=K_{\mathrm{2}}\cdot K_{max}\cdot F_{\mathrm{1}}\left(T_l\right)\cdot F_{\mathrm{1}}\left(W_{\mathrm{sat},l}\right)\\ {\displaystyle }& {\displaystyle }\cdot F\left(\mathrm{pH}\right)\cdot {\mathrm{NH}}_{\mathrm{4},l}^{+}\text{for nitrification and}\\ {\displaystyle }& {\displaystyle }{F}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{denitrification},l}=r_{\mathrm{mx}\mathrm{2}}\cdot F_{\mathrm{2}}\left(W_{\mathrm{sat},l}\right)\cdot F_{\mathrm{2}}\left(T_l,C_{\mathrm{org}}\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }\cdot {\mathrm{NO}}_{\mathrm{3},l}^{-}\text{for denitrification},\end{array}$$

where $F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{nitrification}}$ and $F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{denitrification}}$ are the N₂O flux related to nitrification and denitrification, respectively, in gN m⁻² d⁻¹ in layer l. K₂ is the fraction of nitrified N lost as N₂O (K₂=0.02), K_max is the maximum nitrification rate of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ ($K_{max}=\mathrm{0.1}$ d⁻¹). F₁(T_l) and F₁(W_sat,l) are response functions of soil temperature and water saturation, respectively, that limit the nitrification rate. F (pH) is the function describing the response of nitrification rates to soil pH, and ${\mathrm{NH}}_{\mathrm{4},l}^{+}$ and ${\mathrm{NO}}_{\mathrm{3},l}^{-}$ the soil ammonium and nitrate concentrations in gN m⁻², respectively. F₂(T_l,C_org) and F₂(W_sat,l) are reactions for soil temperature, soil carbon, and water saturation and r_mx2 is the fraction of denitrified N lost as N₂O (11 %, the remainder is lost as N₂). For a detailed description of the N-related processes implemented in LPJmL, we refer the reader to von Bloh et al. (2018).

3.3 Water fluxes

3.3.1 Litter interception

Precipitation and applied irrigation water in LPJmL5 is partitioned into interception, transpiration, soil evaporation, soil moisture, and runoff (Jägermeyr et al., 2015). To account for the interception and evaporation of water by surface litter, the water can now also be captured by surface litter through litter interception (I_surf) and be lost through litter evaporation, subsequently infiltrates into the soil and/or forms surface runoff. Litter moisture (θ) is calculated in the following way:

$$\begin{array}{}\text{(7)}& {\mathit{\theta }}_{t+\mathrm{1}}=min({\mathrm{whc}}_{\mathrm{surf}}-{\mathit{\theta }}_{\left(t\right)},I_{\mathrm{surf}}\cdot f_{\mathrm{surf}}).\end{array}$$

f_surf is calculated by adapting the equation from Gregory (1982) that relates the amount of surface litter (dry matter) per square metre (m²) to the fraction of soil covered:

$$\begin{array}{}\text{(8)}& f_{\mathrm{surf}}=\mathrm{1}-\exp (-A_{\mathrm{m}}\cdot {\mathrm{OM}}_{\text{litter,surf}}),\end{array}$$

where OM_litter,surf is the total mass of dry matter surface litter in grammes per square metre (g m⁻²) and A_m is the area covered per mass of crop specific residue (m² g⁻¹). The total mass of surface litter is calculated assuming a fixed C to organic matter (OM) ratio of 2.38 (CF_OM,litter), based on the assumption that 42 % of the organic matter is C, as suggested by Brady and Weil (2008):

$$\begin{array}{}\text{(9)}& {\mathrm{OM}}_{\text{litter,surf}}={\mathrm{C}}_{\text{litter,surf}}\cdot {\mathrm{CF}}_{\text{OM,litter}},\end{array}$$

where C_litter,surf is the amount of C stored in the surface litter pool in grammes of carbon per square metre (gC m⁻²). We apply the average value of 0.004 for A_m from Gregory (1982) to all materials, neglecting variations in surface litter for different materials. whc_surf (mm) is the water holding capacity of the surface litter and is calculated by multiplying the litter mass with a conversion factor of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ mm kg⁻¹ (OM_litter,surf) following Enrique et al. (1999).

3.3.2 Soil infiltration

The presence of surface litter enhances infiltration of precipitation or irrigation water into the soil, as soil crusting is reduced and preferential pathways are affected (Ranaivoson et al., 2017). In order to account for improved infiltration with the presence of surface litter, we follow the approach by Jägermeyr et al. (2016), which has been developed for implementing in situ water harvesting, e.g. by mulching in LPJmL. The infiltration rate (In in mm d⁻¹) depends on the soil water content of the first layer and the infiltration parameter p;

$$\begin{array}{}\text{(10)}& \text{In}=\text{prir}\cdot \sqrt[p]{\mathrm{1}-{\displaystyle \frac{W_{\mathrm{a}}}{W_{\mathrm{sat},l=\mathrm{1}}-W_{\mathrm{pwp},l=\mathrm{1}}}}},\end{array}$$

where prir is the daily precipitation and applied irrigation water in millimetres, W_a the available soil water content in the first soil layer, and $W_{\mathrm{sat},l=\mathrm{1}}$ and $W_{\mathrm{pwp},l=\mathrm{1}}$ the soil water contents at saturation and permanent wilting point of the first layer in millimetres. By default p=2, but four different levels are distinguished ($p=\mathrm{3},\mathrm{4},\mathrm{5},\mathrm{6}$) by Jägermeyr et al. (2016), in order to account for increased infiltration based on the management intervention. To account for the effects of surface litter, we here scale the infiltration parameter p between 2 and 6, based on the fraction of surface litter cover (f_surf);

$$\begin{array}{}\text{(11)}& p=\mathrm{2}\cdot (\mathrm{1}+f_{\mathrm{surf}}\cdot \mathrm{2}).\end{array}$$

Surplus water that cannot infiltrate forms surface runoff and enters the river system.

3.3.3 Litter and soil evaporation

Evaporation (E_surf in millimetres) from the surface litter cover (f_surf) is calculated in a similar manner as evaporation from the first soil layer (Schaphoff et al., 2018a). Evaporation depends on the vegetation cover (f_v), the radiation energy for the vaporization of water (PET), and the water stored in the surface litter that is available to evaporate (ω_surf) relative to whc_surf. Here, f_surf is also taken into account so that the fraction of soil uncovered is subject to soil evaporation as described in Schaphoff et al. (2018a);

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{surf}}=\text{PET}\cdot \mathit{\alpha }\cdot max(\mathrm{1}-f_{\mathrm{v}},\mathrm{0.05})\cdot {\mathit{\omega }}_{\mathrm{surf}}^{\mathrm{2}}\cdot f_{\mathrm{surf}},\text{(13)}& {\displaystyle }& {\displaystyle }{\mathit{\omega }}_{\mathrm{surf}}=\mathit{\theta }/{\mathrm{whc}}_{\mathrm{surf}},\end{array}$$

where PET is calculated based on the theory of equilibrium evapotranspiration (Jarvis and McNaughton, 1986) and α the empirically derived Priestley–Taylor coefficient (α=1.32) (Priestley and Taylor, 1972).

The presence of litter at the soil surface reduces the evaporation from the soil (E_soil). E_soil (millimetres) corresponds to the soil evaporation as described in Schaphoff et al. (2018a), and depends on the available energy for vaporization of water and the available water in the upper 0.3 m of the soil (ω_evap). However, with the implementation of tillage, the fraction of f_surf now also influences evaporation, i.e. greater soil cover (f_surf) results in a decrease in E_soil:

$$\begin{array}{}\text{(14)}& E_{\mathrm{soil}}=\text{PET}\cdot \mathit{\alpha }\cdot max(\mathrm{1}-f_{\mathrm{v}},\mathrm{0.05})\cdot {\mathit{\omega }}^{\mathrm{2}}\cdot (\mathrm{1}-f_{\mathrm{surf}}),\end{array}$$

where ω is calculated as the evaporation-available water (ω_evap) relative to the water holding capacity in that layer (whc_evap):

$$\begin{array}{}\text{(15)}& \mathit{\omega }=min\left(\mathrm{1},{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{evap}}}{{\mathrm{whc}}_{\mathrm{evap}}}}\right),\end{array}$$

where ω_evap is all the water above wilting point of the upper 0.3 m (Schaphoff et al., 2018a).

3.4 Heat flux

The temperature of the surface litter is calculated as the average of soil temperature of the previous day (t) of the first layer ($T_{\mathrm{soil},l=\mathrm{1}}$ in ^∘C) and actual air temperature ($T_{\mathrm{air},t+\mathrm{1}}$ in ^∘C), in the following way:

$$\begin{array}{}\text{(16)}& T_{\mathrm{litter},\mathrm{surf},t+\mathrm{1}}=\mathrm{0.5}(T_{\mathrm{air},t+\mathrm{1}}+T_{\mathrm{soil},l=\mathrm{1},t}).\end{array}$$

Equation (16) is an approximate solution for the heat exchange described by Schaphoff et al. (2013). The new upper boundary condition (T_upper in ^∘C) is now calculated by the average of T_air and T_surf weighted by f_surf. With the new boundary condition, the cover of the soil with surface litter diminishes the heat exchange between soil and atmosphere;

$$\begin{array}{}\text{(17)}& T_{\mathrm{upper}}=T_{\mathrm{air}}\cdot \left(\mathrm{1}-f_{\mathrm{surf}}\right)+T_{\mathrm{surf}}\cdot f_{\mathrm{surf}}.\end{array}$$

The remainder of the soil temperature computation remains unchanged from the description of Schaphoff et al. (2013).

3.5 Tillage effects on physical properties

3.5.1 Dynamic calculation of hydraulic properties

Previous versions of the LPJmL model used static soil hydraulic parameters as inputs, computed following the pedotransfer function (PTF) by Cosby et al. (1984). Different methods exist to calculate soil hydraulic properties from soil texture and SOM content for different points of the water retention curve (Balland et al., 2008; Saxton and Rawls, 2006; Wösten et al., 1999) or at continuous pressure levels (Van Genuchten, 1980; Vereecken et al., 2010). Extensive reviews of PTFs and their application in Earth system and soil modelling can be found in Van Looy et al. (2017) and Vereecken et al. (2016). We now introduce an approach following the PTF by Saxton and Rawls (2006), which was included in the model in order to dynamically simulate layer-specific hydraulic parameters that account for the amount of SOM in each layer, constituting an important mechanism of how hydraulic parameters are affected by tillage (Strudley et al., 2008).

As such, Saxton and Rawls (2006) define a PTF most suitable for our needs and capable of calculating all the necessary soil water properties for our approach: The PTF allows for a dynamic effect of SOM on soil hydraulic properties, and is also capable of representing changes in bulk density after tillage and was developed from a large number of data points. With this implementation, soil hydraulic properties are now all updated daily. Following Saxton and Rawls (2006), soil water properties are calculated as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{pwp},l}=-\mathrm{0.024}\cdot a+\mathrm{0.0487}\cdot \mathrm{Cl}+\mathrm{0.006}\cdot {\mathrm{SOM}}_l+\mathrm{0.005}\\ {\displaystyle }& {\displaystyle }\cdot \mathrm{Sa}\cdot {\mathrm{SOM}}_l-\mathrm{0.013}\cdot \mathrm{Cl}\cdot {\mathrm{SOM}}_l+\mathrm{0.068}\\ \text{(18)}& {\displaystyle }& {\displaystyle }\cdot \mathrm{Sa}\cdot \mathrm{Cl}+\mathrm{0.031},\text{(19)}& {\displaystyle }& {\displaystyle }{W}_{\mathrm{pwp},l}=\mathrm{1.14}\cdot {\mathit{\lambda }}_{\mathrm{pwp},l}-\mathrm{0.02},{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{fc},l}=-\mathrm{0.251}\cdot \mathrm{Sa}+\mathrm{0.195}\cdot \mathrm{Cl}+\mathrm{0.011}\cdot {\mathrm{SOM}}_l+\mathrm{0.006}\\ {\displaystyle }& {\displaystyle }\cdot \mathrm{Sa}\cdot {\mathrm{SOM}}_l-\mathrm{0.027}\cdot \mathrm{Cl}\cdot {\mathrm{SOM}}_l+\mathrm{0.452}\\ \text{(20)}& {\displaystyle }& {\displaystyle }\cdot \mathrm{Sa}\cdot \mathrm{Cl}+\mathrm{0.299},\text{(21)}& {\displaystyle }& {\displaystyle }{W}_{\mathrm{fc},l}=\mathrm{1.238}\cdot {\left({\mathit{\lambda }}_{\mathrm{fc},l}\right)}^{\mathrm{2}}+\mathrm{0.626}\cdot {\mathit{\lambda }}_{\mathrm{fc},l}-\mathrm{0.015},{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{sat},l}=\mathrm{0.278}\cdot \mathrm{Sa}+\mathrm{0.034}\cdot \mathrm{Cl}+\mathrm{0.022}\cdot {\mathrm{SOM}}_l-\mathrm{0.018}\\ {\displaystyle }& {\displaystyle }\cdot \mathrm{Sa}\cdot {\mathrm{SOM}}_l-\mathrm{0.027}\cdot \mathrm{Cl}\cdot {\mathrm{SOM}}_l-\mathrm{0.584}\\ \text{(22)}& {\displaystyle }& {\displaystyle }\cdot \mathrm{Sa}\cdot \mathrm{Cl}+\mathrm{0.078},\text{(23)}& {\displaystyle }& {\displaystyle }{W}_{\mathrm{sat},l}=W_{\mathrm{fc},l}+\mathrm{1.636}\cdot {\mathit{\lambda }}_{\mathrm{sat},l}-\mathrm{0.097}\cdot \mathrm{Sa}-\mathrm{0.064},\text{(24)}& {\displaystyle }& {\displaystyle }{\mathrm{BD}}_{\mathrm{soil},l}=(\mathrm{1}-W_{\mathrm{sat},l})\cdot \mathrm{MD}.\end{array}$$

SOM_l is the soil organic matter content in weight percent (wt %) of layer l; W_pwp,l is the moisture content at the permanent wilting point; W_fc,l is moisture contents at field capacity; W_sat,l is the moisture contents at saturation; λ_pwp,l, λ_fc,l, and λ_sat,l are the moisture contents for the first solution at permanent wilting point, field capacity, and saturation, respectively; Sa is the sand content in volume percent (vol %); Cl is the clay content in volume percent (vol %); BD_soil,l is the bulk density in kilogrammes per cubic metre (kg m⁻³); and MD is the mineral density of 2700 kg m⁻³. For SOM_l, total SOC content is translated into SOM of this layer:

$$\begin{array}{}\text{(25)}& {\mathrm{SOM}}_l={\displaystyle \frac{{\mathrm{CF}}_{\mathrm{OM},\mathrm{soil}}\cdot ({\mathrm{C}}_{\mathrm{fastSoil},l}+{\mathrm{C}}_{\mathrm{slowSoil},l})}{{\mathrm{BD}}_{\mathrm{soil},l}\cdot z_l}}\cdot \mathrm{100},\end{array}$$

where CF_OM,soil is the conversion factor of 2 as suggested by Pribyl (2010), assuming that SOM contains 50 % SOC; C_fastSoil,l is the fast decaying C pool in kilogrammes per square metre (kg m⁻²); C_slowSoil,l is the slow decaying C pool (kg m⁻²); BD_soil,l is the bulk density in kg m⁻³; and z is the thickness of layer l in metres. It was suggested by Saxton and Rawls (2006) that the PTF should not be used for SOM contents above 8 %, so we cap SOM_l at this maximum when computing soil hydraulic properties and thus treated soils with SOM_l content above this threshold as soils with 8 % SOM content. Saturated hydraulic conductivity is also calculated following Saxton and Rawls (2006) as

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }{\mathrm{Ks}}_l=\mathrm{1930}\cdot {\left(W_{{\mathrm{sat}}_{\left(l\right)}}-W_{{\mathrm{fc}}_{\left(l\right)}}\right)}^{\mathrm{3}-{\mathit{\varphi }}_l},\text{(27)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_l={\displaystyle \frac{\ln \left(W_{\mathrm{fc},l}\right)-\ln \left(W_{\mathrm{pwp},l}\right)}{\ln \left(\mathrm{1500}\right)-\ln \left(\mathrm{33}\right)}},\end{array}$$

where Ks_l is the saturated hydraulic conductivity in millimetres per hour (mm h⁻¹) and ϕ_l is the slope of the logarithmic tension–moisture curve of layer l.

3.5.2 Bulk density effect and reconsolidation

The effects of tillage on BD are adopted from the APEX model by Williams et al. (2015), which is a follow-up development of the EPIC model (Williams et al., 1983). Tillage causes changes in BD of the tillage layer (first topsoil layer of 0.2 m) after tillage. Soil moisture content for the tillage layer is updated using the fraction of change in BD. Ks_l is also updated based on the new moisture content after tillage. A mixing efficiency parameter (mE), depending on the intensity and type of tillage (0–1), determines the fraction of change in BD after tillage. A mE of 0.90, for example, represents a full inversion tillage practice, also known as conventional tillage (White et al., 2010). The parameter mE can be used in combination with residue management assumptions to simulate different tillage types. It should be noted that Williams et al. (1983) calculated direct effects of tillage on BD, while we changed the equation accordingly to account for the fraction at which BD is changed.

The fraction of BD change after tillage is calculated in the following way:

$$\begin{array}{}\text{(28)}& f_{\mathrm{BDtill},t+\mathrm{1}}=f_{\mathrm{BDtill},t}-\left(f_{\mathrm{BDtill},t}-\mathrm{0.667}\right)\cdot \mathrm{mE}.\end{array}$$

Tillage density effects on saturation and field capacity follow Saxton and Rawls (2006):

$$\begin{array}{}\text{(29)}& {\displaystyle }& {\displaystyle }{W}_{\mathrm{sat},\mathrm{till},l,t+\mathrm{1}}=\mathrm{1}-\left(\mathrm{1}-W_{\mathrm{sat},l,t}\right)\cdot f_{\mathrm{BDtill},t+\mathrm{1}},\text{(30)}& {\displaystyle }& {\displaystyle }{W}_{\mathrm{fc},\mathrm{till},l,t+\mathrm{1}}=W_{\mathrm{fc},l,t}-\mathrm{0.2}\cdot (W_{\mathrm{sat},l,t}-W_{\mathrm{sat},\mathrm{till},l,t+\mathrm{1}}),\end{array}$$

where $f_{\mathrm{BDtill},t+\mathrm{1}}$ is the fraction of density change of the topsoil layer after tillage, f_BDtill,t is the density effect before tillage, $W_{\mathrm{sat},\mathrm{till},l,t+\mathrm{1}}$ and $W_{\mathrm{fc},\mathrm{till},l,t+\mathrm{1}}$ are adjusted moisture contents at saturation and field capacity after tillage, and $W_{\mathrm{sat},l,t}$ and $W_{\mathrm{fc},l,t}$ are the moisture content at saturation and field capacity before tillage.

Reconsolidation of the tilled soil layer is accounted for following the same approach by Williams et al. (2015). The rate of reconsolidation depends on the rate of infiltration and the sand content of the soil. This ensures that the porosity and BD changes caused by tillage gradually return to their initial value before tillage. Reconsolidation is calculated the following way:

$$\begin{array}{}\text{(31)}& {\displaystyle }& {\displaystyle }\mathrm{sz}=\mathrm{0.2}\cdot \mathrm{In}\cdot {\displaystyle \frac{\mathrm{1}+\mathrm{2}\cdot \mathrm{Sa}/(\mathrm{Sa}+e^{\mathrm{8.597}-\mathrm{0.075}\cdot \mathrm{Sa}})}{z_{\mathrm{till}}^{\mathrm{0.6}}}},\text{(32)}& {\displaystyle }& {\displaystyle }f={\displaystyle \frac{\mathrm{sz}}{\mathrm{sz}+e^{\mathrm{3.92}-\mathrm{0.0226}\cdot \mathrm{sz}}}},\text{(33)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{BDtill},t+\mathrm{1}}=f_{\mathrm{BDtill},t}+f\cdot (\mathrm{1}-f_{\mathrm{BDtill},t}),\end{array}$$

where sz is the scaling factor for the tillage layer and z_till is the depth of the tilled layer in metres. This allows for a faster settling of recently tilled soils with high precipitation and for soils with a high sand content. In dry areas with low precipitation, and for soils with a low-sand content, the soil settles slower and might not consolidate back to its initial state. This is accounted for by taking the previous bulk density before tillage into account. The effect of tillage on BD can vary from year to year, but f_BDtill,t cannot be below 0.667 or above 1 so that unwanted amplification is not possible. We do not yet account for fluffy soil syndrome processes (i.e. when the soil does not settle over time) and negative implications from this, which results in an unfavourable soil particle distribution that can cause a decline in productivity (Daigh and DeJong-Hughes, 2017).

4 Model set-up

4.1 Model input, initialization, and spin-up

In order to bring vegetation patterns and SOM pools into a dynamic equilibrium stage, we make use of a 5000-year spin-up simulation of only natural vegetation, which recycles the first 30 years of climate input following the procedures of von Bloh et al. (2018). For simulations with land-use inputs and to account for agricultural management, a second spin-up of 390 years is conducted to account for historical land-use change, which is introduced in the year 1700. The spatial resolution of all input data and model simulations is 0.5^∘. Land-use data are based on crop-specific shares of MIRCA2000 (Portmann et al., 2010) and cropland and grassland time series since 1700 from HYDE3 (Klein Goldewijk et al., 2010) as described by Fader et al. (2010). As per default setting, intercrops are grown on all set-aside stands in all simulations (Bondeau et al., 2007). As we are here interested in the effects of tillage on cropland, we ignore all natural vegetation in grid cells with cropland by scaling existing cropland shares to 100 %. We drive the model with daily mean temperature from the Climate Research Unit (CRU TS version 3.23; University of East Anglia Climate Research Unit, 2015; Harris et al., 2014), monthly precipitation data from the Global Precipitation Climatology Centre (GPCC Full Data Reanalysis version 7.0; Becker et al., 2013), and shortwave downward and net longwave downward radiation data from the ERA-Interim data set (Dee et al., 2011). Static soil texture classes are taken from the Harmonized World Soil Database (HWSD) version 1.1 (Nachtergaele et al., 2009) and aggregated to 0.5^∘ resolution by using the dominant soil type. Twelve different soil textural classes are distinguished according to the USDA soil texture classification and one unproductive soil type, which is referred to as “rock and ice”. Soil pH data are taken from the WISE data set (Batjes, 2005). The NOAA/ESRL Mauna Loa station (Tans and Keeling, 2015) provides atmospheric CO₂ concentrations. Deposition of N was taken from the ACCMIP database (Lamarque et al., 2013).

4.2 Simulation options and evaluation set-up

The new tillage management implementation allows for specifying different tillage and residue systems. We conducted four contrasting simulations on current cropland area with or without the application of tillage and with or without removal of residues (Table 1). The default setting for conventional tillage is mE =0.9 and TL =0.95. In the tillage scenario, tillage is conducted twice a year, at sowing and after harvest. Soil water properties are updated on a daily basis, enabling the tillage effect to be effective from the subsequent day onwards until it wears off due to soil settling processes. The four different management settings (MSs) for global simulations are as the following: (1) full tillage and residues left on the field (T_R), (2) full tillage and residues are removed (T_NR), (3) no-till practice and residues are retained on the field (NT_R), and (4) no-till practice and residues are removed from the field (NT_NR). The specific parameters for these four settings are listed in Table 1. The default MS is T_R and was introduced in the second spin-up from the year 1700 onwards, as soon as human land use is introduced in the individual grid cells (Fader et al. 2010). All of the four MS simulations were run for 109 years, starting from year 1900. Unless specified differently, the outputs of the four different MS simulations were analysed using the relative differences between each output variable using T_R as the baseline MS;

$$\begin{array}{}\text{(34)}& {\mathrm{RD}}_X={\displaystyle \frac{X_{\mathrm{MS}}}{X_{\mathrm{T}\mathrm{\_}\mathrm{R}}}}-\mathrm{1},\end{array}$$

where RD_X is the relative difference between the management scenarios for variable X and X_MS and X_{T_R} are the values of variable X of the MS of interest and the baseline management systems: conventional tillage with residues left on the field (T_R). Spin-up simulations and relative differences for Eq. (34) were adjusted if a different MS was used as reference system, e.g. if reference data are available for comparisons of different MSs. The effects were analysed for different time scales: the 3-year average of years 1 to 3 for short-term effects, the average after years 9 to 11 for mid-term effects, and the average of years 19 to 21 for long-term effects. Depending on available reference data in the literature, the specific duration and default MS of the experiment were chosen. The results of the simulations are compared to literature values from selected meta-analyses. Meta-analyses allow for the comparison of globally modelled results to a set of combined results of individual studies from all around the world, assuming that the data basis presented in meta-analyses is representative. A comparison to individual site-specific studies would require detailed site-specific simulations making use of climatic records for that site and details on the specific land-use history. Results of individual site-specific experiments can differ substantially between sites, which hampers the interpretation at larger scales. We calculated the median and the 5th and 95th percentiles (values within brackets) between MS in order to compare the model results to the meta-analyses, where averages and 95 % confidence intervals (CIs) are mostly reported. We chose medians rather than arithmetic averages to reduce outlier effects, which is especially important for relative changes that strongly depend on the baseline value. If region-specific values were reported in the meta-analyses, e.g. climate zones, we compared model results of these individual regions, following the same approach for each study, to the reported regional value ranges.

Table 1 LPJmL simulation settings and tillage parameters used in the stylized simulations for model evaluation.

^a Litter cover is calculated following Gregory (1982). ^b Litter amounts and litter cover are modelled internally.

To analyse the effectiveness of selected individual processes (see Fig. 1) without confounding feedback processes, we conducted additional simulations of the four different MSs on bare soil with uniform dry matter litter input (simulation NT_NR_bs and NT_R_bs1 to NT_R_bs5) of uniform composition (C:N ratio of 20), no atmospheric N deposition and static fertilizer input (Elliott et al., 2015). This helps with isolating soil processes, as any feedbacks via vegetation performance are eliminated in this setting.

5 Evaluation and discussion

5.1 Tillage effects on hydraulic properties

Table 2 presents the calculated soil hydraulic properties of tillage for each of the soil classes prior to and after tillage (mE of 0.9), combined with a SOM content in the tilled soil layer of 0 % and 8 %. In general, both tillage and a higher SOM content tend to increase whc, W_sat,l, W_fc,l, and Ks_l. Clay soils are an exception, since higher SOM content decreases whc, W_sat,l, and W_fc,l, and increases Ks_l. The effect of increasing SOM content on whc, W_sat,l, and W_fc,l is greatest in the soil classes sand and loamy sand. The increasing effects of tillage on the hydraulic properties are generally weaker compared to an increase in SOM by 8 % (maximum SOM content for computing soil hydraulic properties in the model). While tillage (mE of 0.9, 0 % SOM) in sandy soils increase whc by 83 %, 8 % of SOM can increase whc in an untilled soil by 105 % and in a tilled soil by 84 %. As comparison, in silty loam soils with 0 % SOM, tillage (mE of 0.9) increases whc by 16 %, while 8 % SOM can increase whc by 31 % and by 26 % for untilled and tilled soil, respectively.

Table 2 Percentage values for each soil textural class of silt, sand and clay content used in LPJmL and correspondent hydraulic parameters before and after tillage with 0 % and 8 % SOM using the Saxton and Rawls (2006) pedotransfer function.

^a Soil class rock is not affected by SOM changes and tillage practices. ^b For SOM we only consider the C part in SOM in grammes of carbon per square metre (gC m⁻²). ^c Tillage with a mE of 0.9 for conventional tillage. ^d whc
is calculated as whc $=W_{\mathrm{fc}}-W_{\mathrm{pwp}}$ in all cases.

The PTF by Saxton and Rawls (2006) uses an empirical relationship between SOM, soil texture, and hydraulic properties derived from the USDA soil database, implying that the PTF is likely to be more accurate within the US than outside. A PTF developed for global-scale application is, to our knowledge, not yet developed. Nevertheless PTFs are used in a variety of global applications, despite the limitations to validate at this scale (Van Looy et al., 2017).

5.2 Productivity

In our simulations, adopting NT_R slightly increases productivity for all rain-fed crops simulated (wheat, maize, pulses, rapeseed) on average, but ranges from increases to decreases across all cropland globally. This increase can be observed for the first 3 years (Fig. S2 in the Supplement), and for the first 10 years (Fig. 2a and b). All the results shown here and in the subsequent sections are calculated as RD following Eq. (34), unless otherwise stated. The numbers discussed in this section refer to the productivity after 10 years (average of year 9–11). The largest positive impact can be found for rapeseed, where NT_R results in a median increase of +3.5 % (5th, 95th percentiles: −24.5 %, +57.8 %). The positive impact is lowest for maize, with median increases by +1.8 % (5th, 95th percentiles: −24.6 %, +56.2 %). The median productivity of wheat increases slightly by +2.5 % (5th, 95th percentiles: −15.2 %, +53.5 %) under NT_R. The slight increases in median productivity under NT_R are contrasting to the values reported by Pittelkow et al. (2015b), who report slight decreases in productivity for wheat and maize and small increases for rapeseed (Table 3). They report both positive and negative effects for wheat and rapeseed, but only negative effects for maize. Pittelkow et al. (2015b) identify aridity and crop type as the most important factors influencing the responses of productivity to the introduction of no-till systems with residues left on the field. The aridity index was determined by dividing the mean annual precipitation by potential evaporation. No-till performed best under rain-fed conditions in dry climates (aridity index <0.65), by which the overall response was equal or positive compared to T_R.

Figure 2 Relative yield changes for rain-fed wheat (a) and rain-fed maize (b) compared to aridity indexes after 10 years NT_R vs. T_R. Low aridity index values indicate arid conditions as the index is defined as mean annual precipitation divided by potential evapotranspiration, following Pittelkow et al. (2015a). Substantial increases in crop yields only occur in arid regions, with aridity indices <0.75.

Table 3 Comparison of simulated model output and literature values from meta-analyses. Values for modelled results are calculated according to Eq. (34) with adjusted default management.

^a Estimated from graph. ^b Time horizon of the study is unclear in the meta-analysis. The average over the first 3 years of model results is taken. ^c Includes conservation tillage. ^d Residue management
for conventional tillage unsure ^e Time horizon not explicitly mentioned by author.

The positive effects on productivity under NT_R in dry regions can also be found in our simulations. For instance, wheat productivity increases substantially under NT_R, whereas this effect diminishes with increases in aridity indexes (Fig. 2a). Similar results are found for maize productivity (Fig. 2b). This positive effect can be attributed to the presence of surface litter, which leads to higher soil moisture conservation through increased water infiltration into the soil and decreases in evaporation. Areas where crop productivity is limited by soil water could therefore potentially benefit from NT_R (Pittelkow et al., 2015a). The influence of climatic condition of no-till effects on productivity was already found by several other studies (e.g. Ogle et al., 2012; Pittelkow et al., 2015a; van Kessel et al., 2013). Ogle et al. (2012) found declines in productivity but these declines were larger in the cooler and wetter climates. Pittelkow et al. (2015a) found only small declines in productivity in dry areas but emphasized that increases in yield can be found when no-till is combined with residues and crop rotation. This was not the case for humid areas (aridity index >0.65); there declines in productivity were larger under no-till regardless of whether residues and crop rotations were applied. Finally, van Kessel et al. (2013) found declines in productivity after adapting to no-till in dry areas (−11 %) and humid areas (−3 %). However, in their analysis it is not clear how crop residues are treated in no-till and tillage (i.e. removed or retained).

Negative effects of NT_R on productivity can be observed in mainly tropical areas. As soil moisture increases in tropical areas under NT_R as well (Fig. 5c), the decline is resulting from a decrease in N availability in the soil (Fig. 5d). Soil moisture drives many N-related processes that can cause a decline in N. For instance, the increase in soil moisture can lead to an increase in denitrification, which decreases the amount of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (which will be discussed further in Sect. 5.5). On the other hand, mineralization can also be reduced if soil moisture is too high. However, the soil-moisture–N availability and yield feedback is complex as many processes are involved.

5.3 Soil C stocks and fluxes

We evaluate the effects of tillage and residue management on simulated soil C dynamics and fluxes for CO₂ emissions from cropland soils, relative change in C input, SOC turnover time, and relative changes in soil and litter C stocks of the topsoil (0.3 m). In our simulation CO₂ emissions initially decrease for the average of the first 3 years by a median value of −11.9 % (5th, 95th percentile: −24.1 %, +2.0 %) after introducing no-till (NT_R vs. T_R) (Fig. S3a in the Supplement) and soil and litter C stocks increase. After 10 years duration (average of year 9–11), however, both CO₂ emissions and soil and litter C stocks are higher under NT_R than under T_R (Fig. 3a, d). Median CO₂ emissions from NT_R compared to T_R increase by +1.7 % (5th, 95th percentile: −17.4 %, +32.4 %) (Fig. 3a), while at the same time median topsoil and litter C also increase by +5.3 % (5th, 95th percentile: +1.4 %, +12.8 %) (Fig. 3d), i.e. the soil and litter C stock has already increased enough to sustain higher CO₂ emissions. There are two explanations for CO₂ increase in the long term: (1) more C input from increased net primary production (NPP) for NT_R or (2) a higher decomposition rate over time under NT_R, due to changes in, for example, soil moisture or temperature. Initially CO₂ emissions decrease almost globally due to increased turnover times under T_R (Fig. S3c in the Supplement), but after 10 years, CO₂ emissions start to increase in drier regions, while they still decrease in most humid regions (Fig. 3a). The median of the relative differences in mean residence time of soil carbon for NT_R compared to T_R is small, but variable (+0.0 % after 10 years, 5th, 95th percentile: −22.9 %, +23.7 %) (Fig. 3c), and mean residence time shows similar spatial patterns, i.e. it decreases in drier areas but increases in more humid areas. The drier regions are also the areas where we observe a positive effect of reduced evaporation and increased infiltration on plant growth, i.e. in these regions the C input into soils is substantially increased under NT_R compared to T_R (Fig. 3b) (see also Sect. 5.2 for productivity). As such, both mechanisms that affect CO₂ emissions are reinforcing each other in many regions. This is in agreement with the meta-analyses conducted by Pittelkow et al. (2015b), who report a positive effect on yields (and thus general productivity and thus C input) of no-till compared to conventional tillage in dry climates. Their results show that, in general, no-till performs best relative to conventional tillage under water-limited conditions, due to enhanced water-use efficiencies when residues are retained.

Figure 3 Relative C dynamics for NT_R vs. T_R comparison after 10 years of simulation experiment (average of year 9–11) for relative CO₂ change (a), relative C input change (b), relative change of soil C turnover time (c), and relative topsoil and litter C change (d).

Abdalla et al. (2016) reviewed the effect of tillage, no-till, and residue management and found that if residues are returned, no-till compared to conventional tillage increases soil and litter C content by 5.0 % (95th CI: −1.0 %, +9.2 %) and decreases CO₂ emissions from soils by −23.0 % (95th CI: −35.0 %, −13.8 %) (Table 3). These findings of Abdalla et al. (2016) are in line with our findings for CO₂ emissions if we consider the first 3 years of duration for CO₂ emissions and 10 years duration for topsoil and litter C. Abdalla et al. (2016) do not explicitly specify a time of duration for these results. If we only analyse the tillage effect without taking residues into account (T_NR vs. NT_NR), we find in our simulation that topsoil and litter C decreases by −18.0 % (5th, 95th percentile: −42.5 %, −0.5 %) after 20 years, while CO₂ emissions increase by +21.3 % (5th, 95th percentile: −1.1 %, +125.2 %) mostly in humid regions (Table 3). Abdalla et al. (2016) also reported soil and litter C changes from a T_NR vs. NT_NR comparison and reported a decrease in soil and litter C under T_NR of −12.0 % (95th CI: −15.3 %, −5.1 %) and a CO₂ increase of +18.0 % (95th CI: +9.4 %, +27.3 %), which is well in line with our model results.

Ogle et al. (2005) conducted a meta-analysis and reported SOC changes from NT_R compared to T_R system with medium C input, grouped for different climatic zones. They found a +23 %, +17 %, +16 %, and +10 % mean increase in SOC after converting from a conventional tillage to a no-till system for more than 20 years for tropical moist, tropical dry, temperate moist, and temperate dry climates, respectively. We only find a +4.8 %, +8.3 %, +3.5 %, and +5.8 % mean increase in topsoil and litter C for these regions, respectively. However, Ogle et al. (2005) analysed the data by comparing a no-till system with high C inputs from rotation and residues to a conventional tillage system with medium C input from rotation and residues. We compare two similarly productive systems with each other, where residues are either left on the field or incorporated through tillage (NT_R vs. T_R), which may explain why we see smaller relative effects in the simulations. Comparing a high input system with a medium or a low input system will essentially lead to an amplification of soil and litter C changes over time; nevertheless, we are still able to generally reproduce a SOC increase over longer periods.

Unfortunately there are high discrepancies in the literature with regard to no-till effects on soil and litter C, since the high increases found by Ogle et al. (2005) are not supported by the findings of Abdalla et al. (2016). Ranaivoson et al. (2017) found that crop residues left on the field increases soil and litter C content, which is in agreement with our simulation results.

5.4 Water fluxes

We evaluate the effects of tillage and residue management on water fluxes by analysing soil evaporation and surface runoff. Our results show that evaporation and surface runoff under NT_R compared to T_R are generally reduced by −44.3 % (5th, 95th percentiles: −64.5, −17.4 %) and by −57.8 % (5th, 95th percentiles: −74.6 %, −26.1 %), respectively (Fig. S4a and b in the Supplement). We also analysed soil evaporation and surface runoff for different amounts of surface litter loads and cover on bare soil without vegetation in order to compare our results to literature estimates from field experiments. We find that both the reduction in evaporation and surface runoff are dependent on the residue load, which translates into different rates of surface litter cover.

On the process side, water fluxes highly influence plant productivity and are affected by tillage and residue management (Fig. 1). Surface litter, which is left on the surface of the soil, creates a barrier that reduces evaporation and also increases the rate of infiltration into the soil. Litter that is incorporated into the soil through tillage loses this function to cover the soil. Both the reduction of soil evaporation and the increase in rainfall infiltration contribute to increased soil moisture and hence plant water availability. The model accounts for both processes. Scopel et al. (2004) modelled the effect of maize residues on soil evaporation calibrated from two tropical sites and found that a presence of 100 g m⁻² of surface litter decreases soil evaporation by −10 % to −15 % in the data, whereas our model shows a median decrease in evaporation of −6.6 % (5th, 95th percentiles: −26.1 %, +20.3 %) globally (Fig. S5a in the Supplement). The effect of a higher amount of surface litter is much more dominate as Scopel et al. (2004) found that 600 g m⁻² surface litter reduced evaporation by approx. −50 %. For the same litter load, our model shows a median decrease in evaporation of −72.6 % (5th, 95th percentiles: −81.5 %, −49.1 %) (Fig. S5b in the Supplement), which is higher than the results found by Scopel et al. (2004). We further analyse and compare our model results to the meta-analysis from Ranaivoson et al. (2017), who reviewed the effect of surface litter on evaporation and surface runoff and other agroecological functions. Ranaivoson et al. (2017) and the studies compiled by them not explicitly distinguish between the different compartments of runoff (e.g. lateral-, surface-runoff). We assume that they measured surface runoff, since lateral runoff is difficult to measure and has to be considered in relation to plot size. In Fig. 4, modelled global results for relative evaporation and surface runoff change for 10 %, 30 %, 50 %, 70 %, and 90 % soil cover on bare soil are compared to literature values from Ranaivoson et al. (2017). Concerning the effect of soil cover on evaporation (Fig. 4a), we find that we are well in line with literature estimates from Ranaivoson et al. (2017) for up to 70 % soil cover, especially when analysing humid climates. For higher soil cover ≥70 %, the model seems to be more in line with literature values for arid regions. Overall for high soil cover of 90 %, the model seems to overestimate the reduction of evaporation. It should be noted that the estimates from Ranaivoson et al. (2017) are only taken from two field studies, which are only representative for the local climatic and soil conditions, since global data on the effect of surface litter on evaporation are not available. The general effect of surface litter on the reduction in soil evaporation is thus captured by the model, but the model seems to overestimate the response at high litter loads. It is not entirely clear from the literature if these experiments have been carried out on bare soil without vegetation. If crops are also grown in the experiments, water can be used for transpiration which is otherwise available for evaporation, which could explain why the model overestimates the effect of surface litter on evaporation on bare soil without any vegetation.

Figure 4 Relative change in evaporation (a) and surface runoff (b) relative to soil cover from surface residues for different soil cover values of 10 %, 30 %, 50 %, 70 %, and 90 % (simulation NT_R_bs1 to NT_R_bs5 vs. NT_NR_bs, respectively). For better visibility, the red and blue boxplots are plotted next to the overall boxplots, but correspond to the soil cover value of the overall simulation (empty boxes).

Ranaivoson et al. (2017) also investigated the runoff reduction under soil cover, but the results do not show a clear picture. In theory, surface litter reduces surface runoff and the literature generally supports this assumption (Kurothe et al., 2014; Wilson et al., 2008), but the magnitude of the effect varies. Figure 4b compares our modelled results under different soil cover to the literature values from Ranaivoson et al. (2017). This shows that modelled results across all global cropland are on the upper end of the effect of surface runoff reduction from soil cover, but they are still well within the range reported by Ranaivoson et al. (2017). The amount of water which is infiltrated (and thus not going into surface runoff) is affected by the parameter p in Eq. (11), which is dependent on the amount of surface litter cover (f_surf). The parameterization of p is chosen to be at the upper end of the approach by Jägermeyr et al. (2016) at full surface litter cover, as this should substantially reduce surface runoff (Tapia-Vargas et al., 2001) and thus increase infiltration rates (Strudley et al., 2008). The parametrization of p can be adjusted if better site-specific information on slope, soils crusting, and rainfall intensity is available.

5.5 N₂O fluxes

Switching from tillage to no-till management with leaving residues on the fields (NT_R vs. T_R) increases N₂O emissions by a median of +20.8 % (5th, 95th percentile: −3.6 %, +325.5 %) (Fig. S6a in the Supplement). The strongest increase is found in the cool temperate zone where the average increase is +23.5 % (5th, 95th percentile: −0.1 %, +664.4 %) (Fig. S6e in the Supplement). The lowest increase is found in the tropical zone +15.8 % (5th, 95th percentile: −7.3 %, +72.1 %) (Fig. S6c in the Supplement).

The increase in N₂O emissions after switching to no-till is in agreement with several literature studies (Linn and Doran, 1984; Mei et al., 2018; van Kessel et al., 2013; Zhao et al., 2016) (Table 3). Mei et al. (2018) reports an overall increase of +17.3 % (95th CI: +4.6 %, +31.1 %), which is in agreement with our median estimate. However, the regional patterns over the different climatic regimes are in less agreement. LPJmL simulations strongly underestimate the increase in N₂O emissions in the tropical zone, whereas simulations overestimate the response in cool temperate and humid zones and to some extent in the warm temperate zone (Table 3).

In general, N₂O emissions are formed in two separate processes: nitrification and denitrification. The increase in N₂O emissions after adapting to NT_R is mainly resulting from denitrification in our simulations (+55.9 %, Fig. 5a). This increase is visible in most of the regions. The N₂O emissions resulting from nitrification decrease mostly (median of −6.0 %, Fig. 5b) but tends to increase in dry areas. The increase in denitrification and decrease in nitrification, results in a decrease in ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (median of −26.4 %), which appears to be stronger in the tropical areas as well (Fig. 5d). The transformation of mineral N to N₂O is not only affected by the nitrification and denitrification rates, but also by substrate availability (${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$). These in turn are affected by nitrification and denitrification rates, but also by other processes, such as plant uptake and leaching. In the Sahel zone, for example, denitrification decreases and nitrification increases, but ${\mathrm{NO}}_{\mathrm{3}}^{-}$ stocks decline, because leaching increases more strongly (Fig. S7 in the Supplement).

Figure 5 Relative changes for the average of the first 3 years of NT_R vs. T_R for denitrification (a), nitrification (b), soil water content (c), and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (d).

In LPJmL, denitrification and nitrification rates are mainly driven by soil moisture and to a lesser extent by soil temperature, soil C (denitrification), and soil pH (nitrification). A strong increase in annually averaged soil moisture can be observed after adapting NT_R (median of +18.9 %, Fig. 5c). Denitrification, as an anoxic process, increases non-linearly beyond a soil moisture threshold (von Bloh et al. 2018), whereas there is an optimum soil moisture for nitrification, which is reduced at low and high soil moisture contents. In wet regions, as in the tropical and humid areas, nitrification is thus reduced by no-till practices, whereas it increases in dryer regions. The increase in soil moisture under NT_R is caused by higher water infiltration rates and reduced soil evaporation (see Sect. 5.4). Also, no-till practices tend to increase bulk density and thus higher relative soil moisture contents (Fig. 1) also affecting nitrification and denitrification rates and therefore N₂O emissions (van Kessel et al., 2013; Linn and Doran, 1984).

Empirical evidence shows that the introduction of no-till practices on N₂O emissions can cause both increases and decreases in N₂O emissions (van Kessel et al., 2013). This variation in response is not surprising, as tillage affects several biophysical factors that influence N₂O emissions (Fig. 1) in possibly contrasting manners (van Kessel et al., 2013; Snyder et al., 2009). For instance, no-till can lower soil temperature exchange between soil and atmosphere through the presence of litter residues, which can reduce N₂O emissions (Enrique et al., 1999). Reduced N₂O emissions under no-till compared to tillage MS can also be observed in the model results, for instance in northern Europe and areas in Brazil (Fig. S6a in the Supplement).

As several biophysical factors are affected, N₂O emissions are characterized by significant spatial and temporal variability. As a result, the estimation of N₂O emissions are accompanied with high uncertainties (Butterbach-Bahl et al., 2013), which hamper the evaluation of the model results (Chatskikh et al., 2008; Mangalassery et al., 2015).

The deviations from the model results compared to the meta-analyses especially for specific climatic regimes (i.e. tropical and cool temperate) require further investigations and verification, including model simulations for specific sites at which experiments have been conducted. The sensitivity of N₂O emissions highlights the importance of correctly simulating soil moisture. However, simulating soil moisture is subject to strong feedback with vegetation performance and comes with uncertainties, as addressed by, for example, Seneviratne et al. (2010). The effects of different management settings (as conducted here) on N₂O emissions and soil moisture therefore requires further analyses, ideally in different climate regimes, soil types, and in combination with other management settings (e.g. N fertilizers). We expect that further studies using this tillage implementation in LPJmL will increase the understanding of management effects on soil nitrogen dynamics. The great diversity in observed responses in N₂O emissions to management options (Mei et al., 2018) renders modelling these effects as challenging, but we trust that the ability of LPJmL5.0-tillage to represent the different components can also help to better understand their interaction under different environmental conditions.

5.6 General discussion

The implementation of tillage into the global ecosystem model LPJmL opens up opportunities to assess the effects of different tillage practices on agricultural productivity and its environmental impacts, such as nutrient cycles, water consumption, GHG emissions, and C sequestration and is a general model improvement to the previous version of LPJmL (von Bloh et al., 2018). The implementation involved (1) the introduction of a surface litter pool that is incorporated into the soil column at tillage events and the subsequent effects on soil evaporation and infiltration, (2) dynamically accounting for SOM content in computing soil hydraulic properties, and (3) simulating tillage effects on bulk density and the subsequent effects of changed soil water properties and all water-dependent processes (Fig. 1).

In general, a global model implementation on tillage practices is difficult to evaluate as effects are often reported to be quite variable, depending on local soil and climatic conditions. The model results were evaluated with data compiled from meta-analyses, which implies several limitations. Due to the limited amount of available meta-analyses, not all fluxes and stocks could be evaluated within the different management scenarios. For the evaluation we focused on productivity, soil and litter C stocks and fluxes, water fluxes, and N₂O dynamics. The sample size in some of these meta-analyses was sometimes low, which may result in biases if an unrepresentative set of climate and soil combinations was tested. Clearly a comparison of a small sample size to simulations of the global cropland is challenging. Nevertheless, the meta-analyses gave the best overview of the overall effects of tillage practices that have been reported for various individual experiments.

We find that the model results for NT_R compared to T_R are generally in agreement with literature with regard to magnitude and direction of the effects on C stocks and fluxes. Despite some disagreement between reported ranges in effects and model simulations, we find that the diversity in modelled responses across environmental gradients is an asset of the model. The underlying model mechanisms, such as the initial decrease in CO₂ emissions after introduction of no-till practices, can be maintained for longer time periods in moist regions, but are inverted in dry regions due to the feedback of higher water availability on plant productivity and reduced turnover times, and generally increasing soil carbon stocks (Fig. 3) are plausible and in line with general process understanding. Certainly, the interaction of the different processes may not be captured correctly and further research on this is needed. We trust that this model implementation representing this complexity allows for further research in this direction. For water fluxes, the model seems to overestimate the effect of surface residue cover on evaporation for high surface cover, but the evaluation is also constrained by the small number of suitable field studies. Effects can also change over time so that a comparison needs to consider the timing, history, and duration of management changes and specific local climatic and soil conditions. The overall effect of NT_R compared to T_R on N₂O emissions is in agreement with literature as well. However, the regional patterns over the different climatic regimes are in less agreement. N₂O emissions are highly variable in space and time and are very sensitive to soil water dynamics (Butterbach-Bahl et al., 2013). The simulation of soil water dynamics differs per soil type as the calculation of the hydraulic parameters is texture specific. Moreover, these parameters are now changed after a tillage event. The effects of tillage on N₂O emissions, as well as other processes that are driven by soil water (e.g. CO₂, water dynamics), can therefore be different per soil type. The soil-specific effects of tillage on N₂O and CO₂ emissions was already studied by Abdalla et al. (2016) and Mei et al. (2018). Abdalla et al. (2016) found that differences in CO₂ emissions between tilled and untilled soils are largest in sandy soils (+29 %), whereas the differences in clayey soils are much smaller (+12 %). Mei et al. (2018) found that clay content <20 % significantly increases N₂O emissions (+42.9 %) after adapting to conservation tillage, whereas this effect for clay content >20 % is smaller (+2.9 %). These studies show that soil-type-specific tillage effects on several processes can be of importance and should be investigated in more detail in future studies. The interaction of all relevant processes is complex, as seen in Fig. 1, which can also lead to high uncertainties in the model. Again, we think that this model implementation captures substantial aspects of this complexity and thus lays the foundation for further research.

It is important to note that not all processes related to tillage and no-till are taken into account in the current model implementation. For instance, NT_R can improve soil structure (e.g. aggregates) due to increased faunal activity (Martins et al., 2009), which can result in a decrease in BD. Although tillage can have several advantages for the farmer, e.g. residue incorporation and topsoil loosening, it can also have several disadvantages. For instance, tillage can cause compaction of the subsoil (Bertolino et al., 2010), which results in an increase in BD (Podder et al., 2012) and creates a barrier for percolating water, leading to ponding and an oversaturated topsoil. Strudley et al. (2008), however, observed diverging effects of tillage and no-till on hydraulic properties, such as BD, Ks, and whc for different locations. They argue that affected processes of agricultural management have complex coupled effects on soil hydraulic properties; also, that variations in space and time often lead to higher differences than the measured differences between the management treatments. They further argue that characteristics of soil type and climate are unique for each location, which cannot simply be transferred from one field location to another. A process-based representation of tillage effects as in this extension of LPJmL allows for further studying of management effects across diverse environmental conditions, but also to refine model parameters and implementations where experimental evidence suggests disagreement.

One of the primary reasons for tillage, weed control, is also not accounted for in LPJmL5.0-tillage or in other ecosystem models. As such, different tillage and residue management strategies can only be assessed with respect to their biogeochemical effects, but only partly with respect to their effects on productivity and not with respect to some environmental effects (e.g. pesticide use). Our model simulations show that crop yields increase under no-till practices in dry areas but decrease in wetter regions (Fig. 2). However, the median response is positive, which may be in part because the water saving effects from increased soil cover with residues are overestimated or because detrimental effects, such as competition with weeds, are not accounted for.

The included processes now allow us to analyse long-term feedbacks of productivity on soil and litter C stocks and N dynamics. Nevertheless the results need to be interpreted carefully, due to the capacity of the model and implemented processes. We also find that the modelled impacts of tillage are very diverse in space as a result of different framing conditions (soil, climate, management) and feedback mechanisms, such as improved productivity in dry areas if residue cover increases plant-available water. In LPJmL5.0-tillage, the process-based representation of tillage and residue management, and the effects on water fluxes such as evaporation and infiltration at the global scale, is unique in the context of global biophysical models (e.g. Friend et al., 2014; LeQuéré et al., 2018). Future research on improved parameterization and the implementation of a more detailed representation of tillage processes, the effects on soil water processes, changes in porosity and subsoil compaction, and effects on biodiversity and soil N dynamics are needed in order to better assess the impacts of tillage and residue management at the global scale. The spatial resolution needed to resolve processes, such as erosion, data availability, and model structure, need to be considered in further model development (Lutz et al., 2019). As such, some processes, such as a detailed representation of soil crusting processes, may remain out of reach for global-scale modelling.

6 Conclusions

We described the implementation of tillage-related processes into the global ecosystem model LPJmL5.0-tillage. The extended model was tested under different management scenarios and evaluated by comparing to reported impact ranges from meta-analyses on C, water, and N dynamics, as well as on crop yields.

We find that mostly arid regions benefit from a no-till management with leaving residues on the field, due to the water saving effects of surface litter. We are able to broadly reproduce reported tillage effects on global stocks and fluxes, as well as regional patterns of these changes with LPJmL5.0-tillage, but deviations in N fluxes need to be further examined. Not all effects of tillage – including one of its primary reasons, weed control – could be accounted for in this implementation. Uncertainties mainly arise because of the multiple feedback mechanisms affecting the overall response to tillage, especially as most processes are affected by soil moisture. The processes and feedbacks presented in this implementation are complex and evaluation of effects is often limited to the availability of reference data. Nonetheless, the implementation of more detailed tillage-related mechanics into the LPJmL global ecosystem model improves our ability to represent different agricultural systems and to understand management options for climate change adaptation, agricultural mitigation of GHG emissions, and sustainable intensification. We trust that this model implementation and the publication of the underlying source code will promote research on the role of tillage for agricultural production, its environmental impact, and global biogeochemical cycles.