Soil nitrous oxide emissions from global land ecosystems and their drivers within the LPJ-GUESS model (v4.1)

Nitrogen (N) transformation processes by soil microbes account for significant nitrous oxide (N₂O) emissions from natural ecosystems and cropland. However, understanding and quantifying global soil N₂O emissions and their responses to changing environmental conditions remain challenging. Here, we implemented a soil nitrification–denitrification module into the dynamic vegetation model LPJ-GUESS to estimate N₂O emissions from global lands. The performance of this new development is examined using observed N₂O fluxes from natural-soil and cropland field trials and independent global-scale estimates. LPJ-GUESS broadly reproduces the cumulative N₂O emissions under different climate conditions and N fertilizer applications that are observed in the field experiments, with some deviations in emission seasonality. Globally, simulated soil N₂O emissions from terrestrial ecosystems increase from 5.6±0.2 Tg N yr⁻¹ in the 1960s to 9.9±0.3 Tg N yr⁻¹ in the 2010s, with croplands contributing about two-thirds of the total increase. East Asia and South Asia show the fastest growth rates in N₂O emissions over the study period due to the expansion of fertilized croplands. On a global scale, N fertilization (including synthetic fertilizer and manure use), atmospheric N deposition, and climate change contribute 58 %, 46 %, and 24 %, respectively, to the simulated soil N₂O emissions in the 2010s. Rising CO₂ levels in the atmosphere reduce the simulated emissions by 32 % through increased plant N uptake, whereas land use changes have varied spatial effects on emissions depending on N management intensity after land cover conversion. Our estimates only account for the direct soil N₂O emissions, excluding those from fertilized pastures. This study highlights the importance of environmental factors in influencing global soil N₂O emissions, particularly for assessing greenhouse gas mitigation potential in agricultural ecosystems.

1 Introduction

Understanding how anthropogenic activities that influence greenhouse gas emissions and exchanges affect the climate system is crucial to address the threats posed by global change. While much of the attention has been on the carbon (C) cycle and carbon dioxide (CO₂), interest in non-CO₂ greenhouse gases has grown over the last 2 decades, particularly for nitrous oxide (N₂O) due to its rapidly increasing concentration in the atmosphere (from ∼290 ppb in 1940 to ∼336 ppb in 2022, with a marked rapid growth of 1.3 ppb yr⁻¹ from 1980 to 2022; Tian et al., 2024). This human-induced N₂O increase contributes significantly to global warming as N₂O is about 273 times more potent than CO₂ with regard to warming the atmosphere from a 100-year perspective (Canadell et al., 2021). Additionally, N₂O is an important stratospheric-ozone-depleting substance, potentially increasing surface levels of harmful ultraviolet radiation (World Meteorological Organization, 2022).

The main sources of terrestrial N₂O emissions to the atmosphere are nitrogen (N) transformation processes by soil microbes. Other major non-soil contributors include fossil fuel combustion, inland waters, and biomass burning (Tian et al., 2020). From 2000 to 2016, the global soil release of N₂O through nitrification and denitrification from land ecosystems has been estimated to be 9–13 Tg N yr⁻¹ (Bouwman et al., 2013; Tian et al., 2019; Zaehle, 2013), with tropical rainforest soils being the most important sources due to their high organic matter turnover rates in their warm and moist environment (Stehfest and Bouwman, 2006). Agricultural soils (within and outside the tropics) are also critical owing to their high reactive N levels. A significant portion of the increase in atmospheric N₂O from agriculture can be attributed to N management practices and N applications to croplands, especially mineral fertilizer from the Haber–Bosch process. Estimates indicate that global cropland N₂O emissions have increased by 1–3 times in the last several decades, ranging from 0.4–1.4 Tg N yr⁻¹ in the 1960s to 1.3–3.3 Tg N yr⁻¹ in the 2010s (Tian et al., 2019; Wang et al., 2020; Xu et al., 2020). These soil N₂O emissions are expected to continue to rise due to the growing food demand associated with the increasing human population and changes in per capita consumption (Davidson and Kanter, 2014; Fowler et al., 2013). Significant reductions in agricultural N₂O emissions are required to achieve ambitious climate targets, particularly through changes in land management practices and, ultimately, enhanced N use efficiency – that is, increasing the fraction of N input that is harvested as products (Gu et al., 2023; Springmann et al., 2018; Zhang et al., 2021).

N input, especially the use of industrial N fertilizer, is one of the most reliable predictors to quantify past and present N₂O emissions in agricultural ecosystems, with an estimated emission factor (EF) varying from 0.2 % to 1.8 % of the N applied to the soil (default value in IPCC Tier 1; Hergoualc'h et al., 2019). While large-scale estimates of N₂O emissions have been made based on regional and/or global mean EF values (primarily estimated from various reactive N sources; see Crutzen et al., 2008; Smith et al., 2012), all EF methods assume a linear increase in N₂O in response to N input. However, field experimental evidence shows that emission trends are usually exponential when N fertilizer rates exceed plant needs (Shcherbak et al., 2014; Song et al., 2018). Process-based ecological models that capture soil–vegetation interactions and the effects of environmental drivers on the physiological and biogeochemical processes of plants or crops, soils, and microbial communities provide an approach for quantifying historical and future N₂O emissions on large spatial scales (Butterbach-Bahl et al., 2013). This is due to their mechanistic parameterization of C–N dynamics between vegetation and soils under changing environmental conditions and land use management (Pongratz et al., 2018).

Compared with some well-developed site-specific models from the 1990s (e.g., DayCent, Parton et al., 1996; DNDC, Li et al., 1992), the implementation of N₂O-related processes in global biosphere models, such as dynamic global vegetation models (DGVMs; Cramer et al., 2001), only began in the early years of this century. For instance, Xu-Ri and Prentice (2008) adopted the descriptions of nitrification and denitrification processes from the DNDC model (Li et al., 2000) and introduced them in a simplified way into the DGVM LPJ-DyN. They assumed that N₂O fluxes are mainly regulated by soil moisture and temperature, carbon supply, soil characteristics (such as aeration and texture), and reactive N availability. Likewise, Zaehle et al. (2011) largely followed the approach of Xu-Ri and Prentice (2008) and incorporated inorganic soil N dynamics into the O–CN vegetation model, enhancing it with the inclusion of soil pH and chemo-nitrification processes. To better represent the influence of soil microbes on N transformation, some DGVMs further considered the activity of nitrifiers and denitrifiers by simulating the growth and mortality of the responsible bacteria (e.g., CLM3.5, Saikawa et al., 2013; IBIS, Ma et al., 2022b; TRIPLEX, Zhang et al., 2017b). However, most of these DGVM developments focused predominantly on N₂O emissions in natural soils given the extensive coverage of natural vegetation on Earth's land surface (e.g., forest and grassland; see Huang and Gerber, 2015; Saikawa et al., 2013; Xu-Ri et al., 2012). Considering the growing contribution of the agricultural sector to the global N cycle, DGVMs are also being equipped to account for cropland management options – such as N fertilizer and manure use (Tian et al., 2012; Von Bloh et al., 2018), lime and basalt application (Val Martin et al., 2023), tillage (Ciais et al., 2011; Lutz et al., 2020), cover crops (Ma et al., 2023; Olin et al., 2015b; Porwollik et al., 2022), and residue retention (Ren et al., 2020). These strategies have been shown through field experiments to play crucial roles in regulating N₂O emissions (Abalos et al., 2022; Li et al., 2023; Quemada et al., 2020; Yangjin et al., 2021).

In this study, we implemented soil nitrification–denitrification processes in the DGVM LPJ-GUESS, into which detailed agricultural management practices have been incorporated previously (Olin et al., 2015a; Smith et al., 2014). Building on the concepts introduced by Xu-Ri and Prentice (2008), these processes are parameterized within our model to estimate soil N₂O emissions across major natural-vegetation and cropland land cover classes, allowing us to investigate how N₂O responds to environmental changes and N management. The performance of the updated model is evaluated using field experimental data, published results from other modeling studies, and global inversion estimates. Our objective is to quantify the temporal and spatial pattern of global N₂O fluxes across natural vegetation, pastures, and croplands while exploring the environmental factors that drive changes in N₂O emissions over the historical period.

2 Methods and data

2.1 Model description

The LPJ-GUESS model is a comprehensive, process-based global vegetation simulator designed to study plant–soil interactions and their influence on ecosystem biogeochemical cycling, including C–N dynamics of natural and managed ecosystems under environmental changes (e.g., global warming and rising CO₂ levels; see Smith et al., 2014) and land management (capturing, e.g., conservation agriculture, forest thinning, and clear-cutting; see Lindeskog et al., 2021; Olin et al., 2015b). This is achieved by simulating the physiological and biogeochemical processes of plants on a daily basis (Smith et al., 2014). In the default global configuration, the model's natural-vegetation component includes 12 plant functional types (PFTs), specifically 10 forest and 2 grass types, each characterized by distinct phenological and morphological traits; bioclimatic constraints; and specific strategies for establishment, growth, and mortality. Pasture ecosystems are modeled as a competition between C₃ and C₄ grasses, with 50 % of the aboveground biomass being removed annually to represent grazing effects. To account for internal manure deposition from livestock in pastures, the model assumes that 75 % of the N from the harvested biomass is returned to the soil (Lindeskog et al., 2013). Cropland in the version of LPJ-GUESS used in this study is characterized by six crop functional types (CFTs): two temperate C₃ cereals sown in spring and autumn–winter, a C₄ crop representing maize, a tropical C₃ crop representing rice, and two N-fixing C₃ grain legumes representing soybean and pulses (Ma et al., 2022a; Olin et al., 2015a). These CFTs are simulated as either rain-fed or irrigated, determined by the prescribed fractions provided as external input. Crops are harvested once every year when the required heat units are reached. Agricultural practices – such as tillage intensity, N mineral fertilizer and manure application, crop residue removal, and leguminous and non-leguminous cover crops – are also included (Ma et al., 2023; Olin et al., 2015b; Pugh et al., 2015). For large-scale applications, to reflect the widespread adoption of conventional practices in current global agriculture (Porwollik et al., 2019), the model assumes that all croplands are under tillage management without cover-cropping systems and that 25 % of aboveground crop residue is retained in the fields after harvest. Industrial N fertilizer is added to soils at three different stages of crop growth, with application rates varying by CFT. In contrast, all manure is applied as a single input at crop sowing to reflect the time required for manure N to become available for plant uptake in real-world practices (Olin et al., 2015b).

In the model, the C transfers induced by decomposition between 11 soil organic matter (SOM) pools drive N mineralization–immobilization processes in order to maintain prescribed C:N ratios and mass balances in both the C-donor and C-receiver pools (Smith et al., 2014). Soil mineral N after the processes of mineralization and immobilization (i.e., NH${}_{\mathrm{4}}^{+}$) and biological N fixation (NH${}_{\mathrm{4}}^{+}$), chemical N fertilizer input (NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$), and atmospheric N deposition (NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$) jointly determine the total size of the soil reactive N pool. This pool is initially depleted by vegetation uptake, followed by hydrological N losses and gaseous N emissions from the soils (Wårlind et al., 2014). Mineral N leaching in the model increases linearly with soil NO${}_{\mathrm{3}}^{-}$ concentration and is dynamically adjusted by soil water-holding capacity and percolation rate. Gaseous N losses through soil N transformation are implemented in this study and are described in detail below.

2.2 Representation of gaseous N emissions from the soil

2.2.1 Ammonia (NH₃) volatilization

N losses through NH₃ volatilization significantly affect the concentration of NH${}_{\mathrm{4}}^{+}$ in soils, thereby influencing N₂O emissions through nitrification. Following Xu-Ri and Prentice (2008), production of NH₃ from NH${}_{\mathrm{4}}^{+}$ is simulated as a response function to (1) the exchangeable NH${}_{\mathrm{4}}^{+}$ in the soil solution, (2) soil temperature, and (3) pH level:

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}{\mathrm{NH}}_{\mathrm{3}\mathrm{soil}}={\mathrm{NH}}_{\mathrm{4}\mathrm{soil}\mathrm{\_}\mathrm{solution}}^{+}\cdot f_{\mathrm{t}\mathit{\_}{\mathrm{NH}}_{\mathrm{3}}}\cdot f_{\mathrm{pH}\mathit{\_}{\mathrm{NH}}_{\mathrm{3}}}\\ {\mathrm{NH}}_{\mathrm{4}\mathrm{soil}\mathrm{\_}\mathrm{solution}}^{+}={\mathrm{NH}}_{\mathrm{4}\mathrm{soil}}^{+}\cdot W_{\mathrm{soil}\mathrm{\_}\mathrm{relative}}\end{array},\right.\end{array}$$

where NH_3soil is the NH₃ concentration in the soil (kg N m⁻² d⁻¹); NH${}_{\mathrm{4}\mathrm{soil}\mathrm{\_}\mathrm{solution}}^{+}$ is the exchangeable NH${}_{\mathrm{4}}^{+}$ concentration determined by the pool size of NH${}_{\mathrm{4}\mathrm{soil}}^{+}$ and the relative soil water content (W_{soil_relative}); and $f_{\mathrm{t}\mathit{\_}{\mathrm{NH}}_{\mathrm{3}}}$ and $f_{\mathrm{pH}\mathit{\_}{\mathrm{NH}}_{\mathrm{3}}}$ are the limitations (ranging from 0 to 1) placed on volatilization by soil temperature and pH, respectively, and are taken from the LPJ-DyN (Xu-Ri and Prentice, 2008) and CLM5.0 (Val Martin et al., 2023) models.

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{t}\mathit{\_}{\mathrm{NH}}_{\mathrm{3}}}=\min \left[\mathrm{1},e^{\mathrm{308.56}\cdot \left(\frac{\mathrm{1}}{\mathrm{71.02}}-\frac{\mathrm{1}}{T_{\mathrm{soil}}+\mathrm{46.02}}\right)}\right]\text{(3)}& {\displaystyle }& {\displaystyle }{f}_{{\mathrm{pH}}_{{\mathrm{NH}}_{\mathrm{3}}}}=\left\{\begin{array}{ll}\mathrm{0.6},& \forall \mathrm{pH}\in [\mathrm{0},\mathrm{5})\\ \mathrm{0.6}+\frac{\mathrm{0.4}}{\mathrm{3}}\cdot (\mathrm{pH}-\mathrm{5})& \forall \mathrm{pH}\in [\mathrm{5},\mathrm{8})\\ \mathrm{1.0}& \forall \mathrm{pH}\in [\mathrm{8},\mathrm{\infty })\end{array}\right.\end{array}$$

In the above, T_soil is the soil temperature (°C) at a depth of 25 cm, representing the mean temperature of the 0–50 cm top layer in the two-layer implementation of LPJ-GUESS we used, and pH is the soil pH values provided together with soil mineral properties as external input.

Figure 1 Representation of key transformations of soil mineral N in LPJ-GUESS. N species in white boxes represent the concentration of the transformed gases accumulated in the soil. See the equations in Sect. 2.2.1–2.2.4 for details on the limitation factors during each process.

After NH₃ volatilization, the remaining soil NH${}_{\mathrm{4}}^{+}$ pool will continue to be depleted by nitrification processes in the model (see Fig. 1). Since nitrification and denitrification can occur simultaneously under aerobic and anaerobic conditions, respectively, the concept of the “anaerobic balloon” proposed by Li et al. (2000) is adopted to partition reactive N species (e.g., NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$) into these two soil states, with soil water-filled pore space (WFPS) being used as an indicator (Xu-Ri and Prentice, 2008). According to Li et al. (2000), the size of the anaerobic balloon increases exponentially with soil moisture. Therefore, in the model, the fraction of anaerobic substrates (f_anaero) is simply estimated as an exponential function of WFPS:

$$\begin{array}{}\text{(4)}& \left\{\begin{array}{l}{f}_{\mathrm{anaero}}=\mathrm{0.05}+\frac{\mathrm{0.95}-\mathrm{0.05}}{\mathrm{1}+e^{-\mathrm{7.5}\cdot (\mathrm{WFPS}-\mathrm{0.5})}}\\ \mathrm{WFPS}=\frac{W_{\mathrm{soil}\mathrm{\_}\mathrm{volume}}}{{\mathrm{POR}}_{\mathrm{soil}}}\end{array},\right.\end{array}$$

where W_{soil_volume} is the soil volumetric water content in the top layer (m³ m⁻³; 0–50 cm), and POR_soil is the soil porosity determined by soil physical properties (m³ m⁻³). The partitioning of soil substrates under aerobic (X_{soil_aero}) and anaerobic (X_{soil_anaero}) conditions is estimated using

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{l}{X}_{\mathrm{soil}\mathrm{\_}\mathrm{aero}}=\left(\mathrm{1}-f_{\mathrm{anaero}}\right)\cdot X_{\mathrm{soil}}\\ X_{\mathrm{soil}\mathrm{\_}\mathrm{anaero}}=f_{\mathrm{anaero}}\cdot X_{\mathrm{soil}}\end{array},\right.\end{array}$$

where X_soil is the soil substrate concentration (kg N (or C) m⁻² d⁻¹), representing any of the following in this study: NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$, NO${}_{\mathrm{2}}^{-}$, or labile carbon.

2.2.2 Nitrification

Autotrophic nitrification and heterotrophic nitrification are two distinct biological processes involved in the N transformation in soil ecosystems. We focused solely on representing autotrophic nitrification, which is the dominant process in most natural and agricultural soils (Chapin et al., 2011). The heterotrophic pathway is also more challenging to model as it requires estimation of dissolved organic nitrogen as the main substrate for the responsible nitrifying bacteria. Autotrophic nitrification is an aerobic process wherein NH${}_{\mathrm{4}}^{+}$ undergoes sequential oxidation to NO${}_{\mathrm{2}}^{-}$ and then to NO${}_{\mathrm{3}}^{-}$, producing NO_x and N₂O as intermediates and/or by-products. The initial oxidation step involves two distinct groups of nitrifiers: ammonia-oxidizing bacteria and archaea. The subsequent oxidation of NO${}_{\mathrm{2}}^{-}$ to NO${}_{\mathrm{3}}^{-}$ is facilitated by nitrite-oxidizing bacteria. Due to the current limitation in the model's ability to simulate the growth and mortality of soil microbes, we integrate these two oxidation steps into one single process – i.e., NH${}_{\mathrm{4}}^{+}$ is oxidized to NO${}_{\mathrm{3}}^{-}$ directly – to collectively represent nitrification in LPJ-GUESS (see Fig. 1). The production of NO${}_{\mathrm{3}}^{-}$ through nitrification ($R_{{\mathrm{NO}}_{\mathrm{3}\mathrm{\_}\mathrm{nit}}^{-}}$) is formulated as follows:

$$\begin{array}{}\text{(6)}& R_{{\mathrm{NO}}_{\mathrm{3}\mathrm{\_}\mathrm{nit}}^{-}}=k_{\max \mathrm{\_}\mathrm{nit}}\cdot {\mathrm{NH}}_{\mathrm{4}\mathrm{soil}\mathrm{\_}\mathrm{aero}}^{+}\cdot f_{\mathrm{t}\mathrm{\_}\mathrm{nit}}\cdot f_{\mathrm{WFPS}\mathrm{\_}\mathrm{nit}}\cdot f_{\mathrm{pH}\mathrm{\_}\mathrm{nit}},\end{array}$$

where k_{max_nit} is the maximum nitrification coefficient and is set as a constant of 0.1 based on the experimental data from Khalil et al. (2004); NH${}_{\mathrm{4}\mathrm{soil}\mathrm{\_}\mathrm{aero}}^{+}$ is the aerobic soil NH${}_{\mathrm{4}}^{+}$ concentration after NH₃ volatilization (kg N m⁻² d⁻¹; see Sect. 2.2.1); and f_{t_nit}, f_{WFPS_nit}, and f_{pH_nit} are the limitation factors of soil temperature, moisture, and pH, respectively. Soil temperature plays a crucial role in regulating microbial activities. For nitrite-oxidizing bacteria, 37–39 °C is found to be optimal for substrate oxidation (Taylor et al., 2019), and for ammonia-oxidizing bacteria and archaea, the optimal soil temperature can range from 31 to 42 °C (Ouyang et al., 2017). In the model, the maximum nitrification rate is thus assumed to occur at 38 °C as the average optimal temperature for these three groups of nitrifiers:

$$\begin{array}{}\text{(7)}& f_{\mathrm{t}\mathrm{\_}\mathrm{nit}}={\left({\displaystyle \frac{\mathrm{70}-T_{\mathrm{soil}}}{\mathrm{70}-\mathrm{38}}}\right)}^{\mathrm{12}}\cdot e^{\mathrm{12}\cdot \frac{T_{\mathrm{soil}}-\mathrm{38}}{\mathrm{70}-\mathrm{38}}}.\end{array}$$

Besides soil temperature, soil moisture and pH are also key factors affecting nitrification rates. Gleeson et al. (2010) demonstrated that the activity of nitrifying bacteria decreases rapidly when soil WFPS exceeds 0.6 and stops completely when it surpasses 0.8. Consequently, a three-threshold limitation function is incorporated into LPJ-GUESS to simulate soil moisture influence (Eq. 8). For soil pH constraints, the response function of nitrification rate is adopted from Parton et al. (2001), as implemented in the DayCent model (Eq. 9):

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}& f_{\mathrm{WFPS}\mathrm{\_}\mathrm{nit}}=\\ & \left\{\begin{array}{ll}\mathrm{0.239}\cdot e^{\mathrm{2.38}\cdot \mathrm{WFPS}}& \forall \mathrm{WFPS}\in [\mathrm{0.0},\mathrm{0.6}]\\ \mathrm{1}-\frac{\mathrm{1}}{\mathrm{0.2}}\cdot (\mathrm{WFPS}-\mathrm{0.6})& \forall \mathrm{WFPS}\in (\mathrm{0.6},\mathrm{0.8})\\ \mathrm{0}& \forall \mathrm{WFPS}\in [\mathrm{0.8},\mathrm{1.0}]\end{array},\right.\end{array}\text{(9)}& {\displaystyle }{f}_{\mathrm{pH}\mathrm{\_}\mathrm{nit}}=\mathrm{0.56}+{\displaystyle \frac{\arctan [\mathrm{0.45}\cdot \mathit{\pi }\cdot (\mathrm{pH}-\mathrm{5}\left)\right]}{\mathit{\pi }}}.\end{array}$$

Since the mechanisms of N gas emissions during nitrification are not yet fully understood (Butterbach-Bahl et al., 2013), we adopted the same assumption as other ecosystem and crop models, namely that the gaseous N concentration through nitrification (i.e., R_{gas_nit}) is proportional to the nitrification rate, and we estimate the emissions using

$$\begin{array}{}\text{(10)}& R_{\mathrm{gas}\mathrm{\_}\mathrm{nit}}=F_{\max \mathrm{\_}\mathrm{gas}\mathrm{\_}\mathrm{nit}}\cdot R_{{\mathrm{NO}}_{\mathrm{3}\mathrm{\_}\mathrm{nit}}^{-}},\end{array}$$

where $R_{{\mathrm{NO}}_{\mathrm{3}\mathrm{\_}\mathrm{nit}}^{-}}$ is the soil NO${}_{\mathrm{3}}^{-}$ concentration produced through nitrification in Eq. (6) (kg N m⁻² d⁻¹), and F_{max_gas_nit} is the maximum fraction of nitrified N lost as NO_x and N₂O, a parameter that varies widely between models due to differences in their structural configuration for simulating N gas species and specific nitrification processes (Gabbrielli et al., 2024). We assume this parameter to be a constant of 0.25 based on the MICN model (Ma et al., 2022b) to broadly account for the potential N gas fluxes driven by different groups of nitrifying bacteria. Previous studies have revealed that, during nitrification, NO_x dominates at a WFPS below 0.3, and the ratio of NO_x : N₂O is often close to 1.0 at a WFPS of 0.6 (Davidson et al., 2000; Pilegaard, 2013). Thus, in the model, nitrified gaseous N concentration (i.e., R_{gas_nit}) is partitioned into NO_x and N₂O species using an empirical function of soil WFPS:

$$\begin{array}{}\text{(11)}& \left\{\begin{array}{l}{\mathrm{NO}}_{\mathrm{x}\mathrm{\_}\mathrm{nit}}=P_{{\mathrm{NO}}_x:\left({\mathrm{NO}}_x+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{nit}}\cdot R_{\mathrm{gas}\mathrm{\_}\mathrm{nit}}\\ {{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{\_}\mathrm{nit}}=\left(\mathrm{1}-P_{{\mathrm{NO}}_x:\left({\mathrm{NO}}_x+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathrm{\_}\mathrm{nit}}\right)\cdot R_{\mathrm{gas}\mathrm{\_}\mathrm{nit}}\\ P_{{\mathrm{NO}}_x}:\left({\mathrm{NO}}_x+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathrm{\_}\mathrm{nit}=\mathrm{1}-\frac{\mathrm{0.5}}{\mathrm{1}+e^{-\mathrm{20}\cdot (\mathrm{WFPS}-\mathrm{0.375})}}\end{array},\right.\end{array}$$

where NO_{x_nit} and N₂O_{_nit} are nitrified NO_x and N₂O gases in the soil, respectively (kg N m⁻² d⁻¹), and $P_{{\mathrm{NO}}_x:({\mathrm{NO}}_x+{\mathrm{N}}_{\mathrm{2}}\mathrm{O})\mathit{\_}\mathrm{nit}}$ is the partitioning scheme between NO_x and N₂O taken from previous studies (Goldberg and Gebauer, 2009; Pilegaard, 2013).

2.2.3 Denitrification

Denitrification is a series of reduction reactions driven by different groups of microorganisms in anaerobic conditions. Heterotrophic denitrifying bacteria facilitate the full reduction chain from NO${}_{\mathrm{3}}^{-}$ to molecular nitrogen (N₂) (${\mathrm{NO}}_{\mathrm{3}}^{-}\to {\mathrm{NO}}_{\mathrm{2}}^{-}\to {\mathrm{NO}}_x\to {\mathrm{N}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{N}}_{\mathrm{2}}$), a process known as denitrifier denitrification. In contrast, autotrophic nitrifiers typically convert NH${}_{\mathrm{4}}^{+}$ to NO${}_{\mathrm{2}}^{-}$ under aerobic conditions, but when oxygen becomes scarce, they switch to reducing NO${}_{\mathrm{2}}^{-}$ to NO_x and N₂O and, finally, to N₂, a process known as nitrifier denitrification. Given the high reactivity of NO_x under the reducing conditions that facilitate denitrification (Parton et al., 2001; Schlüter et al., 2025), accurately simulating every single transformation step from NO${}_{\mathrm{2}}^{-}$ to N₂ is challenging due to the interdependent nature of these processes (Ma et al., 2022b). In this study, following the concept of “holes in a pipe” (e.g., Firestone and Davidson, 1989; Davidson et al., 2000; Val Martin et al., 2023), we combine the entire reduction chain from NO${}_{\mathrm{2}}^{-}$ to N₂ into a single step (see Fig. 1) to broadly represent the emissions of all N gases produced during denitrification:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathrm{DENIT}}_{\text{LPJ-GUESS}}& ={\mathrm{NO}}_{\mathrm{3}\mathrm{soil}\mathrm{\_}\mathrm{anaero}}^{-}\to R_{{\mathrm{NO}}_{\mathrm{2}\mathrm{\_}\mathrm{denit}}^{-}}\\ & \to R_{\mathrm{gas}\mathrm{\_}\mathrm{denit}},\end{array}\end{array}$$

where DENIT_LPJ-GUESS is the N reduction chain during denitrification represented in LPJ-GUESS; NO${}_{\mathrm{3}\mathrm{soil}\mathrm{\_}\mathrm{anaero}}^{-}$ is the soil NO${}_{\mathrm{3}}^{-}$ concentration under anaerobic conditions calculated in Eq. (5) (kg N m⁻² d⁻¹); and $R_{{\mathrm{NO}}_{\mathrm{2}\mathrm{\_}\mathrm{denit}}^{-}}$ and R_{gas_denit} are the concentrations of denitrified soil NO${}_{\mathrm{2}}^{-}$ and all N gases (kg N m⁻² d⁻¹), respectively. Since the availability of the particular N oxide (NO${}_{\mathrm{3}}^{-}$ or NO${}_{\mathrm{2}}^{-}$) and the availability of soil labile C are the two dominant drivers controlling the activity of denitrifiers (Weier et al., 1993; Chapin et al., 2011), we incorporate these two factors into the model to reflect their limitations on the denitrification rate, along with soil temperature as an additional constraint:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{R}_{{\mathrm{NO}}_{\mathrm{2}\mathrm{\_}\mathrm{denit}}^{-}}& =k_{\max \mathrm{\_}\mathrm{denit}}\cdot {\mathrm{NO}}_{\mathrm{3}\mathrm{soil}\mathrm{\_}\mathrm{anaero}}^{-}\cdot f_{\mathrm{t}\mathrm{\_}\mathrm{denit}}\cdot f_{\mathrm{LC}\mathrm{\_}\mathrm{denit}}\\ & \cdot f_{{\mathrm{NO}}_{\mathrm{3}}^{-}\mathit{\_}\mathrm{denit}},\end{array}\end{array}$$

where f_{t_denit}, f_{LC_denit}, and $f_{{\mathrm{NO}}_{\mathrm{3}}^{-}\mathit{\_}\mathrm{denit}}$ are limitations (ranging from 0 to 1) on denitrification by soil temperature, available C, and soil NO${}_{\mathrm{3}}^{-}$ concentration, respectively. k_{max_denit} is the maximum denitrification coefficient and can reach 1.0 when there are no limitations by environmental factors or the populations of denitrifying bacteria (Gabbrielli et al., 2024). In LPJ-GUESS, to reflect the absence of limitation due to the growth of denitrifiers, k_{max_denit} is assumed to be a constant value of 0.5, chosen as the middle point of the possible range from 0 to 1. For f_{LC_denit} and $f_{{\mathrm{NO}}_{\mathrm{3}}^{-}\mathit{\_}\mathrm{denit}}$, both response functions are adopted from the DNDC model (Li et al., 1992), following the Michaelis–Menten equation. For f_{t_denit}, an empirical sigmoid function, built on experimental observations (Benoit et al., 2015; Ma et al., 2022b), is used for the parameterization of temperature effects:

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}{f}_{\mathrm{LC}\mathrm{\_}\mathrm{denit}}=\frac{{\mathrm{Rh}}_{\mathrm{soil}\mathrm{\_}\mathrm{anaero}}}{K_{\mathrm{C}}\cdot W_{\mathrm{soil}\mathrm{\_}\mathrm{volume}}}+{\mathrm{Rh}}_{\mathrm{soil}\mathrm{\_}\mathrm{anaero}}\\ f_{{\mathrm{NO}}_{\mathrm{3}}^{-}\mathit{\_}\mathrm{denit}}=\frac{{\mathrm{NO}}_{\mathrm{3}\mathrm{soil}\mathrm{\_}\mathrm{anaero}}^{-}}{K_{\mathrm{N}}\cdot W_{\mathrm{soil}\mathrm{\_}\mathrm{volume}+{{\mathrm{NO}}_{\mathrm{3}\mathrm{soil}\mathrm{\_}\mathrm{anaero}}}^{-}}}\\ f_{\mathrm{t}\mathrm{\_}\mathrm{denit}}=e^{(-\mathrm{1})\cdot \frac{{\left(T_{\mathrm{soil}}-\mathrm{37}\right)}^{\mathrm{2}}}{{\mathrm{25}}^{\mathrm{2}}}}\end{array},\right.\end{array}$$

where K_C and K_N are Michaelis–Menten constants of 0.017 kg C m⁻³ and 0.083 kg N cm⁻³ for labile C and N oxides, respectively (Li et al., 1992). W_{soil_volume} is volumetric soil water content in the top layer (m³ m⁻³; 0–50 cm). Since labile C is not explicitly modeled in LPJ-GUESS, we use soil heterotrophic respiration under anaerobic conditions (Rh_{soil_anaero}; kg C m⁻² d⁻¹) as a surrogate for C availability, following Parton et al. (2001) and Xu-Ri and Prentice (2008).

In addition to soil labile C and N oxides, much experimental evidence has shown that soil pH and moisture are also critical in regulating the denitrification rate, particularly during the transformation process from NO${}_{\mathrm{2}}^{-}$ to N₂ (e.g., Bao et al., 2012; Bergaust et al., 2010; Kool et al., 2011; Lim et al., 2018). Therefore, in the model, after the first reduction step (${\mathrm{NO}}_{\mathrm{3}\mathrm{soil}\mathrm{\_}\mathrm{anaero}}^{-}\to R_{{\mathrm{NO}}_{\mathrm{2}\mathrm{\_}\mathrm{denit}}^{-}}$), soil NO${}_{\mathrm{2}}^{-}$ is further denitrified as N gases (R_{gas_denit}) using the reduction equation as the first step but with added limitations by soil pH and WFPS (Eqs. 15–17). Both of these two response functions ($f_{\mathrm{pH}\mathit{\_}{\mathrm{NO}}_{\mathrm{2}}^{-}\mathit{\_}\mathrm{denit}}$ and $f_{\mathrm{WFPS}\mathit{\_}{\mathrm{NO}}_{\mathrm{2}}^{-}\mathit{\_}\mathrm{denit}}$ below) are established based on experimental data and are taken from Blanc-Betes et al. (2021) and Ma et al. (2022b), respectively:

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}{R}_{\mathrm{gas}\mathrm{\_}\mathrm{denit}}& =k_{\max \mathrm{\_}\mathrm{denit}}\cdot R_{{\mathrm{NO}}_{\mathrm{2}\mathrm{\_}\mathrm{denit}}^{-}}\cdot f_{\mathrm{t}\mathrm{\_}\mathrm{denit}}\cdot f_{\mathrm{LC}\mathrm{\_}\mathrm{denit}}\\ & \cdot f_{{\mathrm{NO}}_{\mathrm{2}}^{-}\mathrm{\_}\mathrm{denit}}\cdot f_{\mathrm{pH}\mathit{\_}{\mathrm{NO}}_{\mathrm{2}}^{-}\mathit{\_}\mathrm{denit}}\cdot f_{\mathrm{WFPS}\mathit{\_}{\mathrm{NO}}_{\mathrm{2}}^{-}\mathit{\_}\mathrm{denit}}\end{array},\text{(16)}& {\displaystyle }{f}_{\mathrm{pH}\mathit{\_}{\mathrm{NO}}_{\mathrm{2}}^{-}\mathit{\_}\mathrm{denit}}=\left\{\begin{array}{ll}\mathrm{0.001}& \forall \mathrm{pH}\in [\mathrm{0},\mathrm{4}]\\ \frac{\mathrm{pH}-\mathrm{4}}{\mathrm{3}}& \forall \mathrm{pH}\in (\mathrm{4},\mathrm{7})\\ \mathrm{1.0}& \forall \mathrm{pH}\in [\mathrm{7},\mathrm{\infty })\end{array},\right.\text{(17)}& {\displaystyle }\begin{array}{rl}& f_{\mathrm{WFPS}\mathit{\_}{\mathrm{NO}}_{\mathrm{2}}^{-}\mathit{\_}\mathrm{denit}}=\mathrm{0.624}\\ & +{\displaystyle \frac{\mathrm{0.8}\cdot \arctan [\mathrm{0.45}\cdot \mathit{\pi }\cdot (\mathrm{10}\cdot \mathrm{WFPS}-\mathrm{8}\left)\right]}{\mathrm{2.85}}}.\end{array}\end{array}$$

Davidson et al. (2000) and Pilegaard (2013) pointed out that NO_x emissions exponentially decrease when soil WFPS exceeds 0.3 and cease entirely at 0.7. We therefore assume that the denitrified gaseous N (R_{gas_denit} in Eq. 18) produces NO_x and N₂O species only when WFPS is below 0.7. Above this threshold, the production shifts to N₂O and N₂ gases instead.

$$\begin{array}{}\text{(18)}& R_{\mathrm{gas}\mathrm{\_}\mathrm{denit}}=\left\{\begin{array}{ll}{\mathrm{NO}}_{\mathrm{x}\mathrm{\_}\mathrm{denit}}+{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{denit}}& \forall \mathrm{WFPS}\in [\mathrm{0.0},\mathrm{0.7})\\ {{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}_{\mathit{\_}\mathrm{denit}}+{\mathrm{N}}_{\mathrm{2}\mathrm{\_}\mathrm{denit}}& \forall \mathrm{WFPS}\in [\mathrm{0.7},\mathrm{1.0}]\end{array}\right.\end{array}$$

In the above, NO_{x_denit}, N₂O_{_denit}, and N_{2_denit} are denitrified NO_x, N₂O, and N₂ in the soil, respectively (kg N m⁻² d⁻¹). For soil WFPS between 0 and 0.7, the denitrified gaseous N is partitioned into NO_x and N₂O using

$$\begin{array}{}\text{(19)}& \left\{\begin{array}{l}{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}_{\mathit{\_}\mathrm{denit}}=P_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}:\left({\mathrm{NO}}_x+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}\cdot R_{\mathrm{gas}\mathrm{\_}\mathrm{denit}}\\ {\mathrm{NO}}_{\mathrm{x}\mathrm{\_}\mathrm{denit}}=\left(\mathrm{1}-P_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}:\left({\mathrm{NO}}_x+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}\right)\cdot R_{\mathrm{gas}\mathrm{\_}\mathrm{denit}}\\ P_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{NO}}_x+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}}=\max [\mathrm{0},(\mathrm{3.2}\cdot \mathrm{WFPS}-\mathrm{0.92})/\\ (\mathrm{3.2}\cdot \mathrm{WFPS}-\mathrm{0.08})]\end{array},\right.\end{array}$$

where $P_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:({\mathrm{NO}}_x+{\mathrm{N}}_{\mathrm{2}}\mathrm{O})\mathit{\_}\mathrm{denit}}$ is the partitioning ratio between NO_x and N₂O, built on the data provided in Davidson et al. (2000) and Pilegaard (2013). Previous studies indicated that low temperature, combined with low pH and low soil water content, reduces the activity of N₂O reductase, thereby increasing the ratio of N₂O to (N₂+N₂O) in the last step of denitrification (Weier et al., 1993; Siljanen et al., 2020). Accordingly, for soil WFPS above 0.7, the partitioning of R_{gas_denit} into N₂O and N₂ is determined by soil temperature, moisture, and pH levels, with emissions being estimated using

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{l}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\mathit{\_}\mathrm{denit}=P_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}}\cdot R_{\mathrm{gas}\mathrm{\_}\mathrm{denit}}\\ {\mathrm{N}}_{\mathrm{2}\mathrm{\_}\mathrm{denit}}=\left(\mathrm{1}-P_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}}\right)\cdot R_{\mathrm{gas}\mathrm{\_}\mathrm{denit}}\\ P_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}}=f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}\mathrm{\_}\mathrm{t}}\\ \cdot f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}\mathrm{\_}\mathrm{WFPS}}\cdot f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}\mathrm{\_}\mathrm{pH}}\end{array},\right.\text{(21)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{l}{f}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}\mathrm{\_}\mathrm{t}}=\frac{\mathrm{1}}{\mathrm{1}+e^{\left(T_{\mathrm{soil}}-\mathrm{5}\right)}/\mathrm{10}}\\ f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}\mathrm{\_}\mathrm{WFPS}}=\mathrm{0.85}-\frac{\mathrm{0.85}-\mathrm{0.10}}{\mathrm{1}+e^{-\mathrm{23}\cdot (\mathrm{WFPS}-\mathrm{0.75})}}\\ f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:\left({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right)\mathit{\_}\mathrm{denit}\mathrm{\_}\mathrm{pH}}=\min \left(\mathrm{1},\mathrm{7.23}\cdot \frac{\mathrm{1}}{e^{\mathrm{0.497}\cdot \mathrm{pH}}}\right)\end{array},\right.\end{array}$$

where $P_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O})\mathit{\_}\mathrm{denit}}$ is the partitioning ratio of N₂O to (N₂+N₂O), jointly constrained by soil temperature ($f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O})\mathit{\_}\mathrm{denit}\mathrm{\_}\mathrm{t}}$), WFPS ($f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O})\mathit{\_}\mathrm{denit}\mathrm{\_}\mathrm{WFPS}}$), and pH value ($f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:({\mathrm{N}}_{\mathrm{2}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O})\mathit{\_}\mathrm{denit}\mathrm{\_}\mathrm{pH}}$). The field-based observed data – soil temperature from Weier et al. (1993), soil WFPS from Maag and Vinther (1996), and soil pH from Liu et al. (2010) and Rochester (2003) – are used to establish the three limitation functions.

2.2.4 Gas diffusion

In the model, N gases transformed through ammonia volatilization, nitrification, and denitrification accumulate in the soil and are ultimately released into the atmosphere. Since soil aeration status and temperature are key factors influencing gas diffusion (Li et al., 1992; Zhang et al., 2017b), a straightforward equation based on these two variables, as suggested by Xu-Ri and Prentice (2008), is adopted to estimate gas transport from the top soil layer to the atmosphere:

$$\begin{array}{}\text{(22)}& f_{\mathrm{t}\mathrm{\_}\mathrm{dif}}=\min \left[\mathrm{1},e^{\mathrm{308.56}\cdot \left(\frac{\mathrm{1}}{\mathrm{68}}-\frac{\mathrm{1}}{T_{\mathrm{soil}}+\mathrm{46}}\right)}\right],\end{array}$$

$$\begin{array}{}\text{(23)}& \left\{\begin{array}{l}{\mathrm{NH}}_{\mathrm{3}\mathrm{\_}\mathrm{gas}}={\mathrm{NH}}_{\mathrm{3}\mathrm{soil}}\cdot f_{\mathrm{t}\mathrm{\_}\mathrm{dif}}\cdot (\mathrm{1}-\mathrm{WFPS})\\ {\mathrm{NO}}_{\mathrm{x}\mathrm{\_}\mathrm{gas}}=\left({\mathrm{NO}}_{\mathrm{x}\mathrm{\_}\mathrm{nit}}+{\mathrm{NO}}_{\mathrm{x}\mathrm{\_}\mathrm{denit}}\right)\cdot f_{\mathrm{t}\mathrm{\_}\mathrm{dif}}\cdot (\mathrm{1}-\mathrm{WFPS})\\ {\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{\_}\mathrm{gas}}=\left({\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{\_}\mathrm{nit}}+{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{\_}\mathrm{denit}}\right)\cdot f_{\mathrm{t}\mathrm{\_}\mathrm{dif}}\cdot (\mathrm{1}-\mathrm{WFPS})\\ {\mathrm{N}}_{\mathrm{2}\mathrm{\_}\mathrm{gas}}={\mathrm{N}}_{\mathrm{2}\mathrm{\_}\mathrm{denit}}\cdot f_{\mathrm{t}\mathrm{\_}\mathrm{dif}}\cdot (\mathrm{1}-\mathrm{WFPS})\end{array},\right.\end{array}$$

where f_{t_dif} is the soil temperature limitation function on gas diffusion, taken from Xu-Ri and Prentice (2008), and NH_{3_gas}, NO_{x_gas}, N₂O_{_gas}, and N_{2_gas} are NH₃, NO_x, N₂O, and N₂ gases released from the soil into the atmosphere, respectively. It should be noted that other N gas processes, such as gas diffusion fluxes between soil layers and NO_x losses due to its rapid oxidation near the plant canopy, are not implemented in the model at the moment. A flowchart of key transformations of soil mineral N in LPJ-GUESS is shown in Fig. 1.

2.3 Model experimental protocol

This study includes three model applications or experiments. First, we evaluate the model performance in reproducing observed N₂O fluxes from various natural and cropland field sites worldwide. Next, a global simulation is performed over the historical period to assess the model results against estimates from other modeling studies and inversion approaches. In the last part, we analyze and discuss the environmental factors driving temporal and spatial changes in historical N₂O emissions in global natural-vegetation, pasture, and cropland ecosystems.

To initialize the model state, we adopted the spin-up procedure following the protocol in our earlier study (Ma et al., 2022a). All simulations are initialized with a 500-year spin-up using atmospheric CO₂ concentrations and N deposition from 1901, along with repeating the de-trended 1901–1920 climate (data sources described below). The spin-up starts with simulations of potential natural vegetation, representing unmanaged forest and grassland ecosystems. In the last 30 years of the spin-up, the model gradually introduces cropland by linearly increasing the cropland fraction from zero to the 1901 historic value. Details of the model experiment protocol are provided below.

2.3.1 Model evaluation at site scale

To assess model performance, field-based N₂O observations from natural and cropland soils were compiled from the published literature. This dataset includes 49 natural-vegetation sites and 55 cropland sites situated between ∼43° S and ∼61° N (Fig. 2a and b). For natural lands, field studies typically reported annual N₂O emissions (kg N ha⁻¹ yr⁻¹) across five ecosystems: tropical forests (13 sites), temperate forests (21 sites), boreal forests (5 sites), grasslands (8 sites), and savannas (2 sites). These measurements were recorded during the period from 1981 to 2010. In contrast, observed data from croplands focused solely on cumulative N₂O emissions over the crop-growing season within four major cropping systems (kg N ha⁻¹ per season; wheat, maize, rice, and legumes). This dataset, covering the period from 1994 to 2020, examined how N₂O fluxes respond to varying N fertilizer inputs and management practices in agricultural soils (Fig. 2c). Detailed information for these sites – including geographic coordinates, experimental periods, and cropping management systems – is provided in Tables S1 and S2 in the Supplement.

Figure 2 Distribution of in situ N₂O observation studies used for model evaluation based on (a) natural vegetation and (b) cropland. For cropland field experiments, overall N₂O flux measurement information – cropping systems, measurement duration, and N fertilizer and conservation agriculture managements (i.e., residue, tillage, and cover crops) – is shown in (c).

Due to the lack of weather and N deposition data for most study sites, an observation-based gridded climate dataset, CRU JRA v2.4 (Harris et al., 2020; Kobayashi et al., 2015), and an atmospheric N deposition dataset simulated by CCMI (NH_x–N and NO_y–N; Tian et al., 2018) were used as inputs to drive LPJ-GUESS, selecting the representative grid cell (0.5°×0.5°) for each experimental site. To maintain equilibrium in soil C and N pools after model spin-up, natural-vegetation and cropland systems were simulated continuously from 1901 onwards at their respective experimental sites. Since N fertilizer management had been established for several years before the start of the N₂O emission trials at most cropland sites, we assumed a 5-year period of N fertilizer use in croplands prior to the field trials, with application rates being consistent with those used during the trials. Throughout the experimental period, simulations based on croplands were performed based on the management details reported in the literature (Table S2), whereas natural-vegetation sites remained unmanaged (i.e., growing under rain-fed and unfertilized conditions without timber or biomass harvests). Additionally, to estimate soil hydraulic properties and evaluate our developed N transformation processes, soil physical characteristics – such as soil pH and texture (i.e., content of sand, silt, and clay) – were collected from the literature and kept constant during the simulation period. The accuracy of the simulated N₂O fluxes was statistically assessed using the adjusted R² (the goodness of fit for the linear regression), mean error (ME), mean absolute error (MAE), and root mean square error (RMSE) across all sites.

2.3.2 Global soil N₂O emissions and their drivers

For global-scale applications, climate variables from the CRU JRA v2.4 dataset – daily temperature (minimum, mean, and maximum), precipitation, solar radiation, specific humidity, and wind speed – were used for driving model simulations, covering the period from 1901 to 2020 at a resolution of 0.5°×0.5° (Harris et al., 2020; Kobayashi et al., 2015). Historical annual atmospheric CO₂ concentration and monthly N deposition rates over the same period were derived from Meinshausen et al. (2020) and Tian et al. (2018), respectively (Fig. S1 in the Supplement). Land use and land cover data spanning from 1901 to 2020 were sourced from HILDA+ (Winkler et al., 2020), initially at a 0.01° resolution and later being aggregated to 0.5°. This dataset provides annually varying proportions of natural vegetation, pasture, and cropland for each grid cell. The crop distribution map, including rain-fed and irrigated fractions per grid cell around the year 2000, was extracted from the MIRCA dataset (Portmann et al., 2010) and aggregated to match the six CFTs simulated in the model (see Sect. 2.1 above). To parameterize soil water characteristics, global gridded soil profile data at 0.5° resolution from WISE3 (Batjes, 2009) were used. Synthetic N fertilizer and manure application rates to crops were obtained from AG Grid (Elliott et al., 2015) and Zhang et al. (2017a), respectively, covering the period of 1901–2015 (Fig. S1). As N fertilizer data were available only until 2015, we assumed that N application rates to croplands during 2016–2020 remained steady at 2015 levels. At present, pasture ecosystems represented in LPJ-GUESS do not receive any N fertilizer inputs.

Table 1 Simulation setups representing the contribution of environmental factors to global soil N₂O emissions (see Sect. 2.3.2).

^a Historical climate (1901–1920) with temperature de-trended, repeated throughout the period of 1901–2020. ^b N fertilization on croplands, including mineral N fertilizer and manure application. Fertilized pasture is not simulated in this study. ^c Historical N inputs between 1901–2015, with the 2015 data being extended to cover the period of 2016–2020.

All model experiments for this part of the study spanned from 1901 to 2020. However, the focus of our analyses was primarily on the period from 1960 to 2020, during which time N fertilizer use became prevalent. In the “Reference” simulation (referred to as S0 in Table 1), the model was driven by a constant recycled 1901–1920 climate, together with 1901 CO₂ concentration, N deposition, land use, and fertilizer inputs, to monitor model drift and internal variation. The result of this run was used to generate background N₂O emissions with minimal human influence. Conversely, the “All_Combined” run (referred to as S1 in Table 1) incorporated dynamic inputs for all these factors from 1901 to 2020, reflecting realistic N₂O emissions due to anthropogenic perturbation and environmental change. In each subsequent simulation (referred to as S2–S6 in Table 1), all but one of the factors were allowed to vary dynamically over time, with one factor being held constant at its 1901 level. This setup was designed to isolate and identify the individual impact of each factor on N₂O emissions by comparing S1 with any runs in S2–S6 (denoted as ΔN₂O_si in Eq. 24). The total change caused by these five factors was calculated as the difference between S1 and S0 simulations (ΔN₂O_all). The relative contribution (%) of every single factor to the total change in N₂O emissions was then determined by the ratio of ΔN₂O_si to ΔN₂O_all using Eq. (24):

$$\begin{array}{}\text{(24)}& \left\{\begin{array}{l}\mathrm{\Delta }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{si}}={\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{s}\mathrm{1}}-{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{si}}i\in [\mathrm{2},\mathrm{6}]\\ \mathrm{\Delta }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{all}}={\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{s}\mathrm{1}}-{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{s}\mathrm{0}}\\ \mathrm{\Delta }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{si}\mathrm{\%}}=\mathrm{\Delta }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{si}}/\mathrm{\Delta }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{all}}\times \mathrm{100}\mathit{\%}i\in [\mathrm{2},\mathrm{6}]\\ \mathrm{\Delta }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{interactive}\mathrm{\%}}=\left(\mathrm{\Delta }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{all}}-\sum \mathrm{\Delta }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{si}}\right)/\\ \mathrm{\Delta }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{all}}\times \mathrm{100}\mathit{\%}i\in [\mathrm{2},\mathrm{6}]\end{array},\right.\end{array}$$

where ΔN₂O_si% represents the relative contribution of each factor to the total change in N₂O emissions, and ΔN₂O_interactive% denotes the interactive effects between these factors, presented as a percentage. Si refers to the five environmental factor simulations of Const_Climate, Const_CO₂, Const_Ndep, Const_Nfert, and Const_LUC given in Table 1.

3 Results

3.1 Model–observation comparisons at site scale

3.1.1 Model performance across all sites

The simulated cumulative N₂O emissions generally showed a good agreement with measurements, with regression slopes ranging from 0.72 to 0.87. However, the model tended to overestimate measured emissions by 65 % globally for natural vegetation and by 11 % for cropland sites (Fig. 3a and b). According to in situ observations of natural lands (Fig. S2), tropical forests were identified as the primary N₂O-emitting sources, showing a mean cumulative flux of 1.23 kg N ha⁻¹ annually. Temperate forests followed with an average flux of 0.52 kg N ha⁻¹ yr⁻¹, and boreal forests had the lowest emissions at 0.12 kg N ha⁻¹ yr⁻¹. The model broadly reproduced a similar regional pattern, with the highest N₂O emissions in the warm tropics and the lowest emissions in the cold boreal region, although it underrepresented N₂O sources in the tropics by 15 % while overestimating emissions in temperate and boreal forests by 84 % and 50 %, respectively. Both field measurements and model experiments indicated grasslands to be a weak N₂O source, with the simulations being 24 % lower than the observed emissions of 0.19 kg N ha⁻¹ yr⁻¹ on average. This pattern of underestimation was also found at the grass-dominant savanna site, where LPJ-GUESS underestimated the reported N₂O fluxes of the field trials by 78 % (Fig. S2).

Figure 3 Comparison of modeled and observed cumulative N₂O emissions from (a) natural-vegetation, (b) cropland, and (c) crop-specific systems across all studied sites. The dashed line is the 1:1 line, and the bold black line is a fitted linear regression. ME (mean error) and MAE (mean absolute error) are shown as percentages, while RMSE (root mean square error) is in kg N ha⁻¹ yr⁻¹ for natural vegetation and in kg N ha⁻¹ per season for cropland. Error bars denote the standard error from different field trial replicates collected from the literature. In (c), marker sizes from large to small indicate descending N fertilizer rates applied to crops.

Compared with natural vegetation, cropping systems showed higher N₂O emissions throughout the growing period, primarily due to their high levels of N fertilizer use (Fig. 3c). While the linear regression slopes for the four simulated cropping systems were not far from 1.0, LPJ-GUESS overestimated the seasonal N₂O fluxes in most cropland measurements. This overestimation was particularly pronounced in field-grown rice trials, wherein the modeled emissions were approximately 50 % higher than the field observations. Conversely, grain legumes – crops that fix atmospheric N and typically receive little N fertilizer in the fields – were the only cropping system showing lower simulated N₂O emissions than measured, with an underestimation of 45 % (Fig. 3c). A positive exponential relationship between N fertilizer input and cumulative N₂O emissions was found across various field trials in observations. Although LPJ-GUESS simulations effectively captured this observed increase in N₂O fluxes with rising N fertilizer levels, some discrepancies were noted. Specifically, the model underestimated emissions in the unfertilized soils by 28 % and overestimated emissions in most highly fertilized trials (>200 kg N ha⁻¹) by 65 %. This overestimation tended to become more pronounced with higher N application rates (Fig. S2).

3.1.2 Seasonal pattern of N₂O emissions in natural vegetation

The seasonal pattern of N₂O emissions showed significant variation across different vegetation ecosystems and between individual years in both field experiments and model simulations. For a tropical rainforest in Brazil (Fig. 4a), observed fluxes of N₂O increased during the rainy season and rapidly decreased over the dry period. The model was able to reproduce this rainfall-induced mean response; however, the simulated peak flux was delayed by 2 months compared with the observations, likely due to a temporal mismatch between the modeled and reported soil moisture during the wet season (Fig. S3; Davidson et al., 2008). In contrast, a tropical montane forest in Indonesia (Fig. 4b), where rainfall is more evenly distributed throughout the year (Purbopuspito et al., 2006), exhibited no distinct seasonal pattern in N₂O fluxes in the simulations or the observations. Soil temperature and water availability jointly influenced the magnitude of N₂O emissions in temperate and boreal ecosystems, with the largest fluxes being observed during summer (June–August) and weaker sources (or occasional sinks, depending on the sites) being observed during the winter season (Fig. 4c–g). LPJ-GUESS did not capture these negative fluxes recorded in the field trials and instead produced near-zero N₂O emissions in cold and dry conditions (Fig. 4f and g), mainly as a result of negative temperatures deactivating nitrification and denitrification during the winter. Soil moisture was identified as the dominant factor controlling the seasonal dynamics of N₂O fluxes in a semi-arid grassland with a sandy loam soil texture, where observed WFPS ranged from 0.01 to 0.48 between June and August (Du et al., 2006). While the model effectively represented this N₂O rise due to increasing WFPS under aerobic conditions, it overestimated total emissions by ∼95 % over the summer season (Fig. 4h). This overestimation primarily resulted from the model simulating a higher WFPS value of 0.37 for this sand-dominant soil compared with the observed average of 0.20 in these 3 months (Fig. S3; Du et al., 2006).

Figure 4 Comparison of modeled and observed monthly N₂O emissions at eight natural-vegetation sites: (a) Tapajós National Forest, Brazil (2.9° S, 54.9° W; 2000–2004; Davidson et al., 2008); (b) tropical montane forest, central Sulawesi, Indonesia (1.4° S, 120.3° E; November 2002–October 2003 in Wuasa; Purbopuspito et al., 2006); (c) White Mountain National Forest, USA (43.9° N, 71.8° W; 1998–2000; Groffman et al., 2006); (d) temperate Cedrus forest, Fukushima, Japan (37.4° N, 140.3° E; 2003–2004; Morishita et al., 2007); (e) temperate Cedrus forest, Oita-1 site, Japan (33.2° N, 131.3° E; 2003–2004; Morishita et al., 2007); (f) boreal spruce forest, Finland (61.3° N, 24.9° E; June–August 2000, 2001, and 2003; Maljanen et al., 2006); (g) alpine grassland, Tibetan Plateau, China (35.1° N, 93.1° E; August 2000–July 2001; Pei et al., 2004); and (h) native semi-arid grassland, Inner Mongolia, China (43.5° N, 114.7° E; 1995, 1998, and 2001–2003; Du et al., 2006). The dashed dark-red lines denote the multi-year average of simulations over the observation period. Dashed lines in lighter colors represent the simulations for individual years. Open circles indicate the observed N₂O fluxes averaged over all measurement years, with error bars showing the maximum and minimum values.

Figure 5 Modeled and observed seasonal pattern of (a–c) N₂O emissions in response to three levels of N fertilizer inputs (0, 150, and 250 kg N ha⁻¹, referred to as N0, N150, and N250), (d) soil moisture (WFPS), and (e) soil temperature in °C at a rain-fed maize field site in China for the cropping season of 2005 (47.6° N, 133.5° E; Song and Zhang, 2009). In (a)–(c), totals of 40 % and 60 % of the mineral N fertilizer are applied at the time of the maize-sowing and maize-jointing stages, respectively, with the implemented managements of conventional tillage, zero residue retention, and no cover crops (Song and Zhang, 2009). The thick dashed blue lines in (d) and (e) denote the simulated mean of three N fertilizer inputs. The thinner dashed lines represent the simulations for individual N treatments.

3.1.3 Cropland N₂O emission response to N fertilization

The model's ability to simulate the observed N₂O flux response to N fertilizer application was assessed using seasonal data from a rain-fed maize field site in northeastern China (47.6° N, 133.5° E). Alongside three levels of N fertilizer inputs (0, 150, and 250 kg N ha⁻¹, denoted as N0, N150, and N250 below), this cropping system was managed with conventional tillage, zero residue retention, and no cover crops (see Song and Zhang, 2009). Over the maize-growing season, the cumulative N₂O fluxes were measured at 0.4, 2.0, and 4.8 kg N ha⁻¹ for the N0, N150, and N250 treatments, respectively, exceeding the simulated results of 0.2, 1.8, and 4.1 kg N ha⁻¹. The model's estimates of the fertilizer-induced N₂O emission factors were 1.1 % for the N150 treatment and 1.6 % for the N250 treatment, which closely matched the measured range of 1.0 %–1.8 %, suggesting a good overall agreement between the model simulations and field experiments regarding the N₂O response to N addition (Fig. 5a–c). However, LPJ-GUESS failed to capture the peak N₂O fluxes at the maize-flowering stage, particularly in the highly fertilized N250 treatment (Fig. 5c). The remaining difference between modeled and measured seasonal dynamics was found within the 3-week period after each application of N fertilizer, with the simulated N₂O rates being much higher than the observed values. This overestimation indicated that the soil N transformation processes in the model were overly sensitive to reactive N input. At this site with silt–clay-textured soil, the simulated seasonal trends of soil WFPS and temperature broadly aligned with the observed variations (Pearson correlation coefficients of 0.42–0.57, p<0.05 for both variables) despite the modeled WFPS values being consistently higher than the observed ones (Fig. 5d and e).

3.2 Global soil N₂O emissions

The modeled global N₂O emissions from the soil to the atmosphere increased steadily from 1960 to 2020, with estimates growing from 5.6±0.2 Tg N yr⁻¹ in the 1960s to 9.9±0.3 Tg N yr⁻¹ by the 2010s (Fig. 6). While natural soils remained as the major sources of N₂O, their contributions to global total emissions declined from 81% to 59 % over this period. In contrast, simulated N₂O emissions in croplands showed a clear upward trend from 1960, coinciding with the widespread use of synthetic N fertilizer. Croplands reached their highest average emission rate, 3.6±0.2 Tg N yr⁻¹, in the 2010s, representing 37 % of global land emissions. Pasture ecosystems were identified as weak sources of N₂O in our simulations, with historical estimates varying between 0.3±0.04 to 0.4±0.05 Tg N yr⁻¹ (Table S3); we note, however, that fertilized pastures were not simulated in the model. Overall, the model's estimates of global soil N₂O emissions since the 1980s – in terms of both magnitude and interannual variability – were broadly consistent with other studies using bottom-up approaches (Global Nitrous Oxide Budget; Tian et al., 2024), process-based modeling (global N₂O Model Intercomparison Project, NMIP; Tian et al., 2019), and atmospheric-inversion methods (Thompson et al., 2019).

Figure 6 Modeled global soil N₂O emissions from natural vegetation, pasture, and cropland by LPJ-GUESS for the period of 1960–2020 (see S1 run in Table 1) compared with global land estimates from the literature. Reported data from Tian et al. (2019) in open circles indicate 10-year average emissions simulated by seven process-based vegetation models, with error bars representing 1 standard deviation. The decadal mean emissions from Tian et al. (2024), shown in red triangles, are derived from bottom-up estimates and exclude the emissions from non-soil components (such as inland water, fossil fuels and industry, and biomass burning). Annual N₂O emissions between 1998–2016 from Thompson et al. (2019), shown with the dashed line, represent the average of three atmospheric-inversion frameworks, with the range indicated by the blue-shaded area.

The modeled map of soil N₂O emissions revealed large spatial variation in the 1960s (Fig. 7a). Simulated N₂O rates as high as 1.5–2.5 kg N ha⁻¹ yr⁻¹ were found in the tropics (such as the Congo Basin) and parts of Europe and the United States, where neither water nor temperature was a critical constraint for nitrification and denitrification processes. Conversely, regions with arid climates or at high latitudes experienced N₂O emissions as low as 0–0.5 kg N ha⁻¹ yr⁻¹ as soil water content or temperature limitations restricted the turnover rates of soil N pools in LPJ-GUESS. At a regional scale, Africa and South America, with their extensive areas of natural vegetation, together accounted for 50 % of simulated global land N₂O emissions in the 1960s. North America and Southeast Asia followed, with contributions of 13 % and 9 %, respectively (Table S3).

Figure 7 Maps of soil N₂O emissions modeled by LPJ-GUESS in the (a) 1960s and (b) 2010s and (c–l) time series of simulated total soil N₂O emissions from natural vegetation, pasture, and cropland at a continental level for the period of 1960–2020 (see S1 run in Table 1). The division of the 10 continents used in this study is given in Fig. S4.

Compared with the 1960s, soil N₂O emission rates increased in most parts of the world during 2011–2020, mainly as a result of the combined effects of environmental changes and N management practices (Fig. 7b). The regions with high N deposition and intensive fertilizer use – such as northern China, India, central Europe, and the eastern United States (Fig. S1) – were simulated to have the highest N₂O rates, ranging from 3.5 to 4.5 kg N ha⁻¹ yr⁻¹. Compared with other regions, East Asia and South Asia showed the fastest growth in emissions between the 1960s and 2010s (Fig. 7c–l), largely due to their expansion of fertilized croplands. From 2011 to 2020, these two regions jointly contributed 31 % to global total emissions, which was slightly higher than the combined 27 % contribution from North America and Europe (Table S3).

3.3 Drivers of increased N₂O emissions

Changes in the temporal–spatial patterns of soil N₂O emissions were influenced by a combination of land use factors, climate variation, and atmospheric-composition changes. Between 1960–2020, the global increase in soil N₂O emissions was primarily driven by the growing use of N fertilizer and manure, as well as N deposition and climate change, which elevated soil levels of reactive N available for N₂O production (Fig. 8). During the 2010s, N fertilization alone contributed 3.2±0.2 Tg N yr⁻¹, representing 58 % of the increased global terrestrial emissions. N deposition and climate change followed, with estimated contributions of 46 % and 24 %, respectively (Table S4). In contrast, rising CO₂ concentrations were found to lower soil N₂O emissions, with the negative effect increasing over time. This reduction was particularly significant in natural-vegetation and pasture ecosystems and was less pronounced in croplands. From 2011 to 2020, the CO₂ effect was simulated to reduce global soil N₂O emissions by 1.83±0.1 Tg N yr⁻¹, roughly offsetting half of the increased emissions due to N fertilizer use.

Figure 8 Modeled contributions (ΔN₂O) of environmental factors (climate change, rising CO₂ levels, N fertilization, N deposition, and land use change) to global soil N₂O emissions between 1960–2020 across various vegetation types: (a) natural vegetation, (b) pasture, (c) cropland, and (d) the aggregate of all three ecosystems. The white bar with slashes is the net emission from all factors' contributions. See Eq. (24) for ΔN₂O calculation.

In model simulations, the impact of land use change on soil N₂O emissions showed significant spatial variation depending on N management intensity after land cover conversion (Fig. 9). For instance, increased N₂O due to land use change was typically found in regions where soils received high reactive N input (such as northern China and central Europe). Conversely, regions with low levels of N fertilizer and manure use – like most countries in Africa and South America – were beneficial in reducing N₂O emissions after the conversion from natural vegetation to croplands. Additionally, a net expansion of natural lands in the northern temperate regions, such as the southeastern US and eastern Europe, was found to contribute to the mitigation of emissions. On a global scale, a land-use-induced reduction in N₂O emissions on natural lands was estimated to be $-\mathrm{0.60}\pm \mathrm{0.05}$ Tg N yr⁻¹ during 2011–2020, while cropland experienced an increase of 0.58±0.04 Tg N yr⁻¹ over the same period.

Figure 9 Modeled maps of the contributions of (a) N fertilization, (b) N deposition, (c) rising CO₂, (d) climate change, and (e) land use change to soil N₂O emissions between 2011–2020. See S1–S6 runs in Table 1 for simulation setups and Eq. (24) for ΔN₂O calculation.

4 Discussion

4.1 Model uncertainties at site scale

Incorporating specific nitrification and denitrification processes, together with agricultural management practices, into LPJ-GUESS led to a good agreement between modeled and observed N₂O emissions at cropland sites despite some overestimations relative to highly fertilized trials. One factor contributing to this overestimation is that some of the processes taking place during crop growth are not well represented in the model. In previous studies (Ma et al., 2022a, 2023), crop yields simulated by LPJ-GUESS under high N fertilizer inputs were lower than observations, indicating an underestimation of both plant N demand and uptake. Consequently, the excess N remaining in the soil would facilitate higher gaseous loss in the model. This can also explain the significant overestimations of cumulative N₂O emissions in rice-cropping systems (Fig. 3), where the simulated growing season was about 1 month shorter than field experiments since the growth phase between rice sowing and transplanting has not been implemented in LPJ-GUESS. Compared to observations, such a reduction in the simulated growing period was expected to produce lower N uptake and higher N₂O emissions.

We found that the model generally underrepresented N₂O sources in the tropics while simultaneously overestimating annual emissions in temperate and boreal forests across all evaluated natural sites (Fig. S2). This discrepancy can be partially attributed to the high levels of soil WFPS simulated for most humid tropical climates with fine- or medium-textured soils (>0.75, not shown), leading to large amounts of N₂ gas (instead of N₂O) being produced at low oxygen concentrations (Davidson et al., 2000; Pilegaard, 2013). A field-based synthetic analysis estimated that the global mean ratio of N₂O to (N₂O + N₂) in natural soils was 0.125 during the denitrification processes (Scheer et al., 2020). Although our simulated ratios of 0.08–0.10 in the tropics and of 0.12–0.15 in a typical temperate site were close to this global mean estimate (Fig. S3), they may show a potential 35 % underestimation of N₂O in tropical climates and a 15 % overestimation in temperate natural vegetation. This issue could be addressed by adjusting the partitioning scheme between N₂O and N₂ in the denitrification processes (see Eq. 21). However, this is currently prevented by the lack of measured N₂ data as accurately determining N₂ fluxes from the soil is challenging due to the high background concentration in the atmosphere (Butterbach-Bahl et al., 2013; dos Reis Martins et al., 2022).

LPJ-GUESS simulated the cumulative N₂O emissions in response to variations in soil moisture, temperature, and N fertilizer input reasonably well, but it did not capture the measured seasonal dynamics well, particularly for the representation of peak N₂O daily fluxes (Figs. 4 and 5). This issue was not unique to our study and has also been reported in other studies using different ecological models (e.g., Gaillard et al., 2018; Huang and Gerber, 2015; Ma et al., 2022b; Val Martin et al., 2023). The poor performance in simulating the variability in daily N₂O fluxes was likely to be due to the missing or incomplete N transformation processes in the model. For instance, N₂O uptake has often been reported in field studies during low-temperature seasons (Fig. 3), especially in natural ecosystems at high latitudes (Brummell et al., 2014). Recent studies indicated that an N₂O sink could result from efficient N₂O consumption by anaerobic microsites during the final reduction step of denitrification processes (Hiis et al., 2024; Sihi et al., 2020; Siljanen et al., 2020). However, this underlying mechanism is not well understood, which limits the model's ability to characterize N₂O uptake. In addition, the large differences between simulated and measured N₂O monthly fluxes were seen during the spring, when the soil underwent seasonal freezing and thawing (Fig. 3). These N₂O increases associated with freeze–thaw cycles are difficult to simulate because of challenges in parameterizing the transient pulses of nutrient availability in the micro-environment (Zhang et al., 2017b). The impacts of soil freeze events on nutrient release – stimulated by microbial mortality and the physical breakdown of soil aggregates – have not been represented in either LPJ-GUESS or other process-based models (Tian et al., 2019) despite many field measurements reporting their importance for annual N₂O emissions (Kazmi et al., 2025; Wagner-Riddle et al., 2024).

Similarly to other modeling studies, LPJ-GUESS adopted empirical modifiers, k_{max_nit} (Eq. 6) and k_{max_denit} (Eq. 13), to simulate maximum nitrification and denitrification rates in order to account for the limitations that cannot currently be clarified. These parameters vary to a large degree between models, and this variation likely represents each model's unique structure and the specific N transformation processes it includes (Gabbrielli et al., 2024). In this study, we used prescribed constants for k_{max_nit} and k_{max_denit} based on values from the literature (see Sect. 2.2). Assessing the impact of this implementation on seasonal dynamics of N₂O is challenging as these experiment-dependent parameters are typically unavailable for each test site. Using fixed values to represent all unclarified nitrification–denitrification situations, such as in our global-uniform parameterization, cannot reflect the high spatial–temporal variability across different sites, which may further affect the evaluated seasonal pattern of N₂O emission at the site level.

4.2 Global soil N₂O emissions and their drivers

Global simulations by LPJ-GUESS indicated that soil N₂O emissions from natural vegetation were 5.9±0.13 Tg N yr⁻¹ in the 2010s, which was comparable with the 6.4 Tg N yr⁻¹ reported by the Global Nitrous Oxide Budget (Tian et al., 2024). Saikawa et al. (2014) used the top-down atmospheric-inversion method to estimate N₂O emissions from natural soils over the period of 1995–2008, yielding an average emission rate of 7.1 Tg N yr⁻¹. Our estimate of 5.8±0.2 Tg N yr⁻¹, although lower, remains within the range of these other findings (4.7–8.4 Tg N yr⁻¹) for that same 14-year time frame.

Cropland N₂O emissions were simulated by LPJ-GUESS to be 3.5±0.25 Tg N yr⁻¹ globally between 2007–2016, which was consistent with the ensemble mean of 3.3 Tg N yr⁻¹ obtained by NMIP models over the same period (Tian et al., 2019). Our estimates were also close to those of Xu et al. (2020) and Val Martin et al. (2023), who reported a range of 3.1–3.2 Tg N yr⁻¹ in N₂O emissions during 2010–2014 using the DLEM and CLM5 models, respectively. However, these process-based models broadly produced higher N₂O emission rates compared with estimates using the global inventory approach. For instance, based on the IPCC Tier-1 emission factor (EF), FAOSTAT predicted a 10-year average of 1.7 Tg N yr⁻¹ in direct soil N₂O emissions from 2011 to 2020 through the combination of emissions from N fertilizer use, manure application, and crop residues (FAOSTAT, 2024). Similarly, the latest Emissions Database for Global Atmospheric Research (EDGAR v8.0; Crippa et al., 2024) reported direct soil N₂O emissions in agricultural lands of 3.0 Tg N yr⁻¹ for the same period, including large contributions from livestock excreta on grazing pastures. The discrepancies between our results and EF-based approaches were most likely due to the background “legacy effect” – resulting from soil mineralization and residual N accumulation from previous years – not being accounted for by IPCC guidelines. Considering this legacy effect (∼1.5 Tg N yr⁻¹; Kim et al., 2013), a recent study re-estimated global cropland N₂O emissions to be 2.6 Tg N yr⁻¹ over 2010–2014 based on an N-rate-dependent EF method (Wang et al., 2020). While our model results are higher, they can thus be seen to be comparable with this statistically derived estimate.

Emission factors from IPCC Tier 1 (default values of 0.4 % for flooded rice and 1 % for other crops; Hergoualc'h et al., 2019) are often used to quantify large-scale soil N₂O emissions caused by reactive N inputs. Davidson (2009) implemented an EF of 2.5 % and reported global synthetic-fertilizer-induced soil emissions to be 2.2 Tg N yr⁻¹ in 2005, including both direct and indirect sources. In this study, we only simulated direct N₂O emissions from fertilization, which amounted to 2.6 Tg N yr⁻¹ in the same year. The higher simulated fertilization effect can be partially attributed to the inclusion of emissions from manure application in the model (S5 in Table 1), whereas Davidson (2009) did not report this agricultural source independently. According to our simulation results, in 2005, the ratio of fertilization effect (2.6 Tg N yr⁻¹) to total fertilizer inputs (synthetic N and manure; 109 Tg N yr⁻¹) was estimated to be 2.4 % globally, consistently with the findings in Davidson (2009). Notably, both studies imply that the default EF of 1 % in IPCC Tier 1 might lower the estimation of direct N₂O emissions in agricultural soils. Similarly with N-applied croplands, recent global meta-analysis has shown that N enrichment (e.g., atmospheric N deposition) significantly increased N₂O emissions by 80 %–101 % in natural ecosystems, particularly in temperate and boreal forests (Deng et al., 2020; Shen and Zhu, 2022). We identified N deposition to be the most important environmental driver of N₂O emissions in our model (S4 in Table 1), contributing to a 126 % increase for natural lands during 2011–2020 (Fig. 8). This finding was in line with the results of the meta-data analyses referenced above.

Understanding of the CO₂ effect on soil N₂O emissions is still incomplete at the global scale. A quantitative assessment by Van Groenigen et al. (2011) found that rising atmospheric CO₂ levels increased N₂O emissions by 25 % in upland natural ecosystems. This increase was attributed to enhanced plant fine-root biomass and soil moisture, which favored the carbon availability for denitrifying bacteria under anaerobic conditions. Conversely, a recent meta-analysis argued that elevated CO₂ concentrations significantly improved plant N use efficiency by as much as 32 %, subsequently decreasing hydrological N loss by 33 % and N₂O emissions by 5 % in global forests (Cui et al., 2024). It remains unclear which of these two opposing mechanisms might play a more dominant role in field measurements, but the negative CO₂ effect on N₂O emissions, at least in our simulations (−87 % during 2011–2020 on global natural lands; Fig. 8), can be explained by enhanced vegetation N uptake and reduced soil mineral N surplus under increased CO₂ conditions (Fig. S5). Studies using other models have shown similar results (e.g., Huang and Gerber, 2015; Tian et al., 2019; Zaehle et al., 2011). However, the situation differs for N-applied agricultural lands, where the high abundance of soil mineral N is typically sufficient for crop uptake and nitrifier–denitrifier use, thus favoring N₂O production even under rising CO₂ levels. The low CO₂ effect on cropland N₂O emissions compared with natural-vegetation emissions in our simulations (−6 % during 2011–2020; Fig. 8) is in line with previous modeling findings that reactive N addition (e.g., N fixation and deposition) may diminish the negative influence of elevated CO₂ on soil N₂O emissions through mitigating N limitation to both vegetation and soil microbes (Kanter et al., 2016; Xu-Ri et al., 2012).

In this study, climate change was simulated to increase global N₂O emissions by 1.3 Tg N yr⁻¹ (+24 %; Fig. 8) during 2011–2020, which was slightly higher than the mean increase of 1.0 Tg N yr⁻¹ found by NMIP models (Tian et al., 2019) and below the 33 % warming-induced enhancement reported in a meta-analysis (Li et al., 2020). This positive effect in both simulations and observations can be explained by the increased soil temperature, which can (a) speed up the N mineralization process, resulting in N-rich substrates for microbe use, and (b) significantly stimulate the activity and population of denitrifying bacteria, which thrive more in warmer environments than nitrifiers (Pärn et al., 2018). Regionally, an N₂O decline due to climate change was seen in parts of India, Malaysia, and Indonesia (Fig. 9), likely resulting from increased rainfall over the 2 most recent decades (Xu et al., 2020). As discussed in Sect. 4.1, high soil water content linked to ample rainfall can suppress N₂O production in the tropics once soil WFPS exceeds 0.7 (Davidson et al., 2000; Pilegaard, 2013). By contrast, the N₂O decreases in semi-arid climates, such as in western Asia, could be attributed to reduced precipitation. In these regions, nitrification is the dominant process determining N₂O production under aerobic conditions and is usually proportional to the simulated soil WFPS (Fig. S3; Davidson et al., 2000).

4.3 Modeling limitations and implications

Implementing and evaluating soil N₂O emissions in models remains challenging due to the short timescales and high spatial heterogeneity of microbial processes in soils, as well as uncertainties from model input forcings, parameterization, and structure. In this study, we incorporated the influence of major controlling factors – soil moisture and temperature, carbon supply, soil texture and pH, and reactive N availability – into nitrification and denitrification processes in LPJ-GUESS, following empirical findings from previous studies. Although the model can represent the cumulative N₂O emissions satisfactorily, the daily magnitude and seasonal variation did not match experimental observations well. This discrepancy probably reflects the differences between highly controlled and/or specific local field conditions and the general protocols adopted for our simulations related to land use history, N interaction between plants and soils, initial SOM levels, and assumptions made regarding other N-related fluxes (e.g., NH₃ volatilization, N leaching, and partitioning ratio between N₂ and N₂O). Moreover, some key processes that regulate N₂O emissions in field trials – such as the life cycle of nitrifiers and denitrifiers, N₂O uptake, heterotrophic nitrification–aerobic denitrification, and freeze–thaw cycles – have not been implemented in the model.

The current LPJ-GUESS version used in this study only models two soil layers (i.e., 0–50 and 50–150 cm; see Smith et al., 2014; Wårlind et al., 2014). Soil organic matter has no defined placement within the soil column, and the N transformation processes implemented in this study are influenced solely by the temperature and moisture of the top layer (see Sect. 2.2). Given the 50 cm thickness, soil water variability during rain events is minimal, which explains the absence of any distinct WFPS peaks in our modeling results (Figs. 5 and S3). In reality, however, most N₂O-related hydrological processes occur very close to the surface, where soil saturation and drying take place rapidly over short periods. Additionally, soil water content is currently simulated between the wilting point and field capacity (Smith et al., 2014), indicating that fully saturated or completely dry soils cannot be represented in the model. This constraint would consequently impact the simulated WFPS values across different soil textures, which, in turn, affects all other N transformation processes relying on WFPS (see Fig. 1). To better account for gaseous N losses from the soil, the improvement in the soil hydrological scheme within LPJ-GUESS remains to be taken into account in future model development.

Both modeling and EF-based approaches indicated that managed grasslands were strong N₂O sources, with emissions of 1.5–2.2 Tg N yr⁻¹ between 2000–2006, including contributions from livestock excreta deposition, manure use, and N fertilizer application (Dangal et al., 2019; Oenema et al., 2005). In this study, LPJ-GUESS did not account for the effects of N management on pasture, which may have resulted in an underestimation of soil N₂O emissions across global land. However, a long-term gridded dataset for fertilizer and manure applications on pasture has recently become available (HaNi; Tian et al., 2022), offering the possibility of incorporating fertilized grasslands into our future model simulations. In addition, concentrations of soil NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$, as essential substrates for nitrifying and denitrifying bacteria, are dominant factors controlling N₂O production on agricultural lands (Fig. 1). The ratio of NH${}_{\mathrm{4}}^{+}$ to NO${}_{\mathrm{3}}^{-}$, which typically varies with fertilizer types, has been reported to significantly influence soil N₂O emissions in field experiments (e.g., Nelissen et al., 2014; Shcherbak et al., 2014). Globally, Nishina et al. (2017) and Tian et al. (2022) pointed out that the NH${}_{\mathrm{4}}^{+}$ : NO${}_{\mathrm{3}}^{-}$ ratio from N fertilizer application gradually increased from 2.0 to 7.0 during 1961–2010 as a result of the increased consumption of urea on croplands. In contrast, LPJ-GUESS assumed this ratio to be a constant of 1.0 when agricultural soils received N fertilizer inputs. Using this fixed parameter in our simulations cannot reflect the variability of fertilizer-type-induced N₂O emissions in reality, particularly in highly fertilized regions. Considering the model's sensitivity to N fertilization (Fig. 3) and the contribution of synthetic N fertilizer to global N₂O increases (Fig. 8), a step forward could be to harmonize the existing fertilizer-species-dependent datasets to reduce the simulated uncertainties of cropland N₂O emissions.

In addition to N fertilizer and manure management, conservation agriculture – like reduced tillage, residue retention, and legume cover crops – has, for many years, been recommended as a promising climate mitigation practice because of its ability to enhance soil carbon sequestration (Poeplau and Don, 2015; Smith et al., 2020). However, much experimental evidence showed that these conservation practices had the potential to offset the CO₂ mitigation effect due to the increased N₂O emissions (e.g., Lugato et al., 2018; Mei et al., 2018; Quemada et al., 2020; Yangjin et al., 2021). In this study, we only simulated conventional management because of the limited adoption of conservation agriculture in current global croplands (see Sect. 2.1). Whether such a trade-off between CO₂ uptake and N₂O emissions due to conservation practices would also emerge on a large scale – particularly in regions with high N fertilizer applications – needs to be investigated in future work.

5 Conclusions

In this study, we implemented mechanistic representations of nitrification and denitrification in LPJ-GUESS to account for soil N₂O emissions in global terrestrial ecosystems. The simulated N₂O fluxes from natural soils and croplands were compared against observations ranging from the site level to the global scale. Our results showed that the N₂O scheme implemented in the model realistically responded to changes in soil moisture, temperature, and reactive N inputs. It produced cumulative N₂O emission rates comparable with measured data despite some deviations in seasonal patterns.

In our simulations, global soil N₂O emissions from land ecosystems showed a rapid increase between 1960–2020, rising from 5.6±0.2 Tg N yr⁻¹ in the 1960s to 9.9±0.3 Tg N yr⁻¹ in the 2010s. While natural vegetation was the predominant N₂O source in the 1960s, its emissions were gradually surpassed by croplands over the study period. During 2011–2020, East Asia emerged as the largest regional source of N₂O, with N fertilization (including synthetic fertilizer and manure use) being identified as the major contributor. On average, global N₂O increases due to N fertilization increased from 0.22±0.13 Tg N yr⁻¹ in the 1960s to 3.2±0.2 Tg N yr⁻¹ in the 2010s. We also found that atmospheric N deposition and climate change have both contributed to the rise in global N₂O emissions, although effects varied significantly among different vegetation types. Conversely, rising CO₂ levels were found to reduce simulated N₂O emissions through increased plant N uptake, whereas land use change had varied spatial effects on emissions depending on how nitrogen was managed after land cover conversion.

Incorporating key transformations of soil mineral N into LPJ-GUESS offers the opportunity to evaluate total N loss from the soil to the atmosphere, which is essential for accurately quantifying the global terrestrial N cycle in response to changing environmental conditions. This representation also facilitates the assessment of climate mitigation potential in global agricultural soils by examining how CO₂ uptake and N₂O emissions may interact under various conservation practices.