Increasing sensitivity of terrestrial nitrous oxide emissions to precipitation variations

Nitrous oxide (N2O), a major greenhouse gas and ozone-depleting agent, is generated over land mostly from two key biochemical processes—nitrification and denitrification. Nitrifying and denitrifying N2O production occurs preferably under alternative oxic and anoxic conditions, which are closely linked with variations in water filled soil pores, and thus indirectly with precipitation. We show here that the interannual anomalies in the annual growth rate of the global land N2O emissions are significantly (P < 0.001) correlated with precipitation anomalies, with an overall sensitivity ( αPRE , changes of land N2O emission variations per precipitation anomalies) of 2.50 ± 0.98 Tg N2O–N per 100 mm of precipitation across the global land (1998–2016). The sensitivity ( αPRE ) and precipitation-driven N2O anomalies increased during 1998–2016, partly due to increased nitrogen inputs to agricultural lands and enhanced precipitation anomalies. Spatially, we find that the αPRE increases with aridity. We predict a larger αPRE under future climate conditions (with radiative forcing levels of 4.5, 7.0 and 8.5 Wm−2) by the year 2100 if nitrogen fertilization follows the present practice.


Introduction
Nitrous oxide (N 2 O) is a major long-lived greenhouse gas and ozone-depleting substance, whose atmospheric concentration has been increasing since the mid-eighteenth century. Human induced surplus (i.e. beyond biological demands) of reactive nitrogen, through, for example, the Haber-Bosch process (the main industrial procedure to produce synthetic nitrogen fertilizers) and expansion of nitrogen-fixing crops, is one of the major drivers of the increasing trend. This trend is expected to continue into the future due to the growing demand for food, feed, fuel and fibre [1]. While there is a strong scientific consensus on the increasing trend of land N 2 O emissions and its origin, the year-to-year variations in the growth rate of emissions and their drivers at the global scale remain largely unknown.
Nitrification and denitrification are the two major biochemical processes (mostly controlled by microorganisms) that produce N 2 O in terrestrial and aquatic ecosystems. Nitrification is an aerobic process from which ammonium is oxidized to form nitrate. N 2 O is produced as a by-product during this process. Denitrification, which includes heterotrophic denitrification, nitrifier denitrification, chemodenitrification (abiotic) and codenitrification (biotic) [2], reduces oxidized forms of nitrogen (e.g. nitrate) to produce a series of gaseous nitrogen species, including N 2 O, and requires low oxygen conditions. Nitrification and denitrification, as well as their associated gaseous N 2 O emissions, are sensitive to multiple environmental drivers. Temperature, soil water, redox potential (or oxygen levels), nitrogen availability (therefore nitrogen fertilizer input), carbon supply and soil pH are common regulators of land N 2 O emissions documented across laboratory, field and modelling studies [2][3][4][5]. How those abiotic and biotic, human induced and natural factors control land N 2 O emissions are complex. N 2 O fluxes are highly heterogeneous in space and time [3,6]. The water filled pore space, sometimes represented by soil moisture, plays a particularly important role as it controls nitrification and denitrification rates through regulating the substrate availability and soil redox potential (oxygen diffusion rates are much slower in water-filled than air-filled pore space), and also alters the partitioning among gaseous nitrification and denitrification products (e.g. NO, N 2 O and N 2 ) [5,7]. Most well-known hot (with disproportionately high intensity) spots and moments of emissions, such as those at the terrestrial-aquatic interface, riparian zones, oxic-anoxic interfaces of drained organic soils, or/and triggered by rainfall and freeze-thaw events, are related to water availability [8].
Recent climate change has altered precipitation regimes across large areas of the global land, resulting in an increasing frequency of both dry and wet spells, as reflected in soil moisture anomalies [9,10]. Field observations show relatively low N 2 O flux under dry conditions and high emissions following rewetting, partly due to the accumulation of organic nitrogen or/and nitrate during the preceding dry periods that favour denitrifying N 2 O productions under wet conditions [11,12]. Using an inverse procedure combining flux tower measurements, Griffis et al [4] showed enhanced N 2 O emissions during warm and wet periods in the US Corn Belt. The notorious heterogeneity of land N 2 O emissions in space and time makes it difficult to upscale sparse local observations to derive global patterns and drivers. Regional and global studies on the response of land N 2 O emissions to precipitation are scarce. From the N 2 O Model Intercomparison Project (NMIP) [13], under idealized simulations where only climate was varied temporally, the anomaly (difference compared to the emission in 1901) of global land N 2 O emissions were simulated to increase with time (Extended Data figure 8 of [14]) in association with warming since the industrial revolution. It remains unclear from this correlation with temperature trends whether the imprint of precipitation on global land N 2 O emissions is muted due to averaging of different trends across space, or whether the signal has yet to be detected.
We hypothesize that the interannual anomalies in the growth rate of the global land N 2 O emissions are correlated with precipitation anomalies. At least two lines of arguments support our hypothesis. Theoretically, nitrifying-denitrifying N 2 O production is favoured under alternative aerobic and anaerobic conditions, which could be mediated by precipitation [3,[15][16][17]. Empirically, laboratory and field observations have revealed the strong linkages between water availability and N 2 O production processes [16][17][18]. Here and afterwards, we use N 2 OEA (nitrous oxide emission anomalies) to refer to the interannual anomalies in the growth rate of land N 2 O emissions, calculated through a 12 month moving sum window over the detrended (linear least-squares fitting) monthly growth rate of emissions [19]. The growth rate of a specific month is the difference between the N 2 O flux in that month and in the same month of the previous year (see Methods) [20]. Climate anomalies (PRE: precipitation; TEMP: temperature; RAD: solar radiation) were calculated using the same procedure.

Land N 2 O emissions
Land N 2 O emissions were derived from three independent atmospheric inversion systems, i.e. TOMCAT (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014), LMDz (1998-2016) and MIROC4-ACTM (1998-2016) [19]. These systems apply the Bayesian inversion to optimally combine observed N 2 O mixing ratios from discrete air samples and in-situ instrument networks with atmospheric transports, which account for processes such as the loss of N 2 O in the stratosphere caused by photolysis and oxidation by O ( 1 D). The inversions generate grid-level land N 2 O emissions at a monthly temporal resolution. The land here also includes the inland water bodies due to the difficulty of separating them with relatively coarse grid cell resolution (>100 km × 100 km). Detailed information for the inversion products is available in [19]. To quantify the global or regional budgets, we integrated the gridded emissions through corresponding land areas. To estimate land N 2 O emissions, we removed anthropogenic emissions from the energy, industry and waste sectors using the PRIMAP-hist historical emission database [21], despite this not having a noticeable effect on the general relationship between anomalies of N 2 O emissions and precipitation (supplementary figure 1). Keeping vs. removing biomass burning emissions using GFED-v4.1 s did not change the global pattern between land N 2 O emissions and precipitation anomalies noticeably (supplementary figure 1). For the global budget, we reported results that exclude biomass burning emissions.
We also analyzed land N 2 O emissions from process-based land models that participated in the NMIP project [13]. Three models (namely, ORCHIDEE-CNP, ORCHIDEE and OCN) provided simulations with the same temporal resolutions as the atmospheric inversions. Model simulations share the same experimental protocol, e.g. with the same climate drivers Climatic Research Unit-National Centers for Environmental Prediction (CRU-NCEP), land use (HYDE version 3.2) and anthropogenic nitrogen inputs. Nitrogen fertilizer input integrates country-level fertilizer statistics (yearly) from the International Fertilizer Industry Association and the Food and Agriculture Organization, with detailed information from [22]. The timing and splitting frequency of fertilizer application depend on model assumptions, such as one application at the beginning of the growing season.

Climate data
Gridded climate data used in this study includes aridity, monthly mean land surface air temperature, precipitation and solar radiation. We used the gridded (0.5 • × 0.5 • ) CRU time-series (TS) data version 4.04 (CRU TS4.04) [23], provided by CRU at the University of East Anglia as our climate dataset for temperature, precipitation and solar radiation. CRU TS4.04 does not include solar radiation, and we therefore used the cloud cover to approximate variations in solar radiation. We also applied the CRU-NCEP dataset (https://rda. ucar.edu/datasets/ds314.3/) as an alternative for temperature, precipitation and especially for solar radiation. CRU-NCEP combines observations from CRU with the NCEP-NCAR reanalysis data products, and was applied as meteorological forcing data in NMIP. Solar radiation in CRU-NCEP was corrected to match empirically derived monthly solar radiations based on latitude and sunshine hours, which is correlated with the CRU cloud fraction. To test the robustness of the relationship between precipitation and land N 2 O emissions and to quantify the uncertainty, we included 3 other precipitation datasets that combine different sources of rain gauge stations, satellites, sounding observations and reanalyses. The Global Precipitation Climatology Centre (GPCC) provides quality-controlled monthly data from 7,000-9,000 stations (www. dwd.de/EN/ourservices/gpcc/gpcc.html). GPCC v2018 was obtained from (https://psl.noaa.gov/data/ gridded/data.gpcc.html) (assessed Jan 2021). GPCPv2.3 (https://climatedataguide.ucar.edu/climate-data/ gpcp-monthly-global-precipitation-climatology-project, assessed Jan 2021) is the Global Precipitation Climatology Project monthly precipitation dataset from 1979-present that combines observations and satellite precipitation data into 2.5 • × 2.5 • global grids. NOAA's Precipitation Reconstruction over Land (NOAA, https://psl.noaa.gov/data/gridded/data.precl.html, assessed Jan 2021) provides 1 • × 1 • global grids of monthly precipitation reconstructed from observations over 17 000 gauge stations.
We separated the global land into hyperarid, arid, semi-arid, dry subhumid and humid regions based on the aridity index. Aridity index here is defined as the ratio between precipitation and the reference evapotranspiration, provided by [24]. Hyperarid region is classified as the land with the smallest aridity index (AI < 0.03) and covers 7.5% of the total global land area according to the United Nations Environment Programme [25]. The class with the second smallest aridity index (AI: 0.03-0.2) is the arid region, which covers 12.1% of the global land, followed by semi-arid (17.7% of land area, AI: 0.2-0.5), dry subhumid (9.9%, AI: 0.5-0.65) and humid regions (52.8%, AI > 0.65). All gridded datasets (climate and N 2 O) were regridded to 1 • × 1 • using the nearest neighbourhood method for consistency.

Anomalies, sensitivities and decompositions
The growth rate of land N 2 O emissions of a specific month was calculated as the difference between N 2 O emissions from that specific month and the month in the previous year. The monthly growth rates were detrended (linear least-squares fitting) and summed every 12 months to estimate the interannual anomalies (N 2 OEA) similarly as in [20]. Interannual anomalies of climate variables (PRE: precipitation; TEMP: temperature; RAD: solar radiation) were calculated in a similar way. Current available global datasets of anthropogenic nitrogen inputs are coarse in temporal resolution (yearly), with the signal dominated by the increasing trend in the past several decades. The direct impact of the anomaly of nitrogen inputs alone is less likely to be the major driver of N 2 OEA as the anomaly signal of nitrogen inputs is small after detrending. Instead, N 2 OEA are more likely to be driven by anomalies in climate or/and its interactions with anthropogenic nitrogen inputs. We calculated the time-lag (within 12 months) in which the highest Pearson correlation coefficient between N 2 OEA and PRE was obtained. We reported results with optimum time-lags. Calculations with vs. without time-lags yield similar results (α PRE and N 2 OEA_PRE) and patterns discovered in this study are robust.
The sensitivity (α PRE ) of N 2 OEA to PRE was quantified as the slope of PRE in the multiple regression of N 2 OEA against anomalies of precipitation, temperature and radiation. The significance of the regression slope was analyzed through the t-test. For the global total N 2 OEA, we tested incorporating total nitrogen fertilizer (both detrended and not detrended) as an independent predictor, or/and as a predictor interacting with different combinations of climate variables. In the tests, we assumed an equal application rate of nitrogen fertilizer each month with only annual nitrogen fertilization statistics and no comprehensive information on the timing of applications at the global scale. We did not find a noticeable improvement in the performance of the regression (R 2 did not increase) and precipitation remained as the dominant predictor. We therefore focused on the first order effects and reported results on the partial linear regression with only precipitation, temperature and radiation as predictors (supplementary figure 2). We calculated α PRE over a moving time window of 96 months and tested the significance of the trends through the Mann-Kendall method. To verify that the trends are not specific to the moving window length, we varied the window length from 48 to 144 months. And we found similar trends. N 2 OEA has an autocorrelation of around 6 months. To test whether our analyses were impacted by the autocorrelation, we calculated α PRE through the generalized linear mixed model (GLMM: allow for flexible nonnormal distributions of data) that explicitly considers the temporal autocorrelation utilizing the first-order autocorrelation function (i.e. accounts for correlations among the consecutive data points.) Considering autocorrelation had a minor impact on α PRE (supplementary figure 3). In addition, to test whether our results are specific to the detrending method, we applied the singular spectrum analysis, a nonparametric spectral estimation method for TS analysis, and found a similar increase in α PRE (supplementary figure 4).
For the contributions of precipitation to N 2 OEA, in addition to α PRE , variations in precipitation also play a role. We decomposed variations of climate-driven N 2 OEA into contributions from precipitation (N 2 OEA_PRE ), temperature (N 2 OEA_TEMP) and radiation (N 2 OEA_RAD) similarly as in [27]: where α d with the superscript is the sensitivity of N 2 OEA to the corresponding climate anomalies (PRE,TEMP and RAD To estimate the uncertainties in α PRE and climate-driven N 2 OEA, we applied the bootstrapping method with 1000 random samples combining different data sources (five precipitation, two radiation and two temperature datasets) for each of the N 2 O datasets. Uncertainties were expressed as the standard deviations of these samples.
To understand the spatial patterns and mechanisms that lead to the enhanced influence of precipitation on N 2 OEA, we calculated α PRE and N 2 OEA_PRE STD across spatially aggregated groups. We first divided the global land into five categories by aridity. For each category, we further separated it into two groups by their magnitude (i.e. standard deviation) of precipitation anomalies (PRE). One group has an enhanced precipitation anomaly, shown by a higher standard deviation of PRE in the recent time window (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) than in the earliest time window (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). And the other group has reduced precipitation anomalies. We also looked into contributions from agricultural lands vs. natural lands. The atmospheric inversion estimates mix agricultural and non-agricultural emissions within a grid cell. As a coarse approximation, we used the fractional area of cropland [28] as the criterion to separate agricultural and non-agricultural lands. We define non-agricultural grid cells as locations with cropland area < 1%, whereas agricultural grid cells have cropland areas bigger than a threshold value. Agricultural grid cells have an average mineral N input rate of 7 gN (m 2 cropland) −1 yr −1 (1998-2013) [22].
We tested the threshold value from 20% to 50%, and found that the results of a generally higher mean ∆α PRE and ∆N 2 OEA_PRE STD in agricultural than non-agricultural lands are robust (supplementary figure 5). ∆ represents the difference between the recent vs. the earliest time window. For process-based modelling estimates, agricultural N 2 O fluxes due to anthropogenic nitrogen inputs were directly calculated from simulations with vs. without anthropogenic nitrogen inputs.
To derive the scale dependence of the relative importance of different climate components, we define the relative importance of a climate component as the mean variance (across space) of the anomalies of this component divided by the mean variance of the climate driven N 2 OEA, following the method from [27]. For example, the relative importance of precipitation anomalies (I PRE ) is quantified as, where s is the index for locations (or grid cells), µ stands for the mean and σ represents the standard deviation.

Results and discussions
We firstly focus on the top-down estimates of land N 2 O emissions from three independent atmospheric inversion systems (TOMCAT, LMDz and MIROC4-ACTM) [19]. We then quantified the sensitivity (α PRE ) of the globally integrated N 2 OEA toPRE as the partial derivative of N 2 OEA in response to PRE in the multiple regression of N 2 OEA against climate variables (PRE, TEMP and RAD). We included radiation as a climate predictor considering especially its role in vegetation dynamics, which regulates the energy (carbon) and substrate (nitrogen) availability for nitrification and denitrification. With an ensemble of estimates combining the three N 2 O inversions and different datasets of climate variables (Methods; PRE: 5 datasets; TEMP: 2 datasets; RAD: 2 datasets), the overall global α PRE is 2.50 ± 0.98 (mean ± standard deviation) Tg N 2 O-N (100 mm) −1 during 1998-2016. That is, a 100 mm rise in precipitation over the global land increases N 2 O emissions by 2.50 ± 0.98 (t-test, P < 0.01) Tg N 2 O-N.
To translate sensitivities to contributions to N 2 OEA during 1998-2016, we decomposed N 2 OEA into the contributions from precipitation (N 2 OEA_PRE), temperature (N 2 OEA_TEMP), and radiation (N 2 OEA_RAD) (see Methods, equations (1) and (2). Climate anomalies explain up to 59% of the globally integrated N 2 OEA. More than 90% of the climate-induced variations in globally integrated N 2 OEA are attributable to precipitation (supplementary figure 2). Note here we focus on the anomaly and removed the trend signal from temperature, which has been hypothesized to be an important regulator of historical land N 2 O emissions. Adding nitrogen fertilization rate as a global explanatory variable (with or without interactions with climate variables) does not improve the regression performance partly due to large uncertainties in the timing and rate of fertilizer additions at the global scale (methods). To facilitate comparisons among locations, time-windows and climate variables, we use the standard deviation of N 2 OEA_PRE (N 2 OEA_TEMP and N 2 OEA_RAD) over specific time periods to indicate the strength of the climate-driven anomalies, which we denote as N 2 OEA_PRE STD (N 2 OEA_TEMP STD and N 2 OEA_RAD STD ) for simplicity.
Spatially, there was a large variation in α PRE in different places. Roughly, we found that α PRE increases with aridity (defined as the ratio between precipitation and the reference evapotranspiration [24], figure 2(b)). Hyperarid regions (e.g. Sahara and Arabian Deserts) have a high α PRE averaged among pixels falling in this aridity class, but this high sensitivity does not translate into a high mean N 2 OEA_PRE STD , due to small changes in the amount and therefore in the absolute anomalies of precipitation in those regions. The standard deviation of α PRE is also high in hyperarid regions likely due to limited observations and large uncertainties in estimating N 2 O emissions. Humid regions, covering 52.8% of the global land area, contribute to the largest portion (54.3%) of the globally integrated N 2 OEA_PRE STD (sum across land area), followed by the semi-arid regions ( figure 2(d)). Semi-arid regions stand out for their highest We further verified the robustness of the enhanced α PRE in the study period through (a) applying the nonparametric spectral estimation method (singular spectrum analysis) to separate the trend and other components instead of the linear least-squares detrending; and (b) varying the moving window length for calculating α PRE from 48 to 144 months. Note the time-span and testing window length of our study are limited by use of atmospheric N 2 O measurements from the year 1998 onwards. Earlier atmospheric data have lower calibration accuracy and measurement precision [19]. We also verified that considering autocorrelation only had a minor impact on the inferred increase of α PRE (Methods, supplementary figure 3). Note α PRE at the beginning of the study is not statistically significant, which may be related to a larger uncertainty before 2000 [19]. After removing data from the first two years, we obtained an initial α PRE of 2.05 ± 1.52 Tg N 2 O-N (100 mm) −1 and N 2 OEA_PRE STD of 0.38 ± 0.25 Tg N 2 O-N yr −1 . The increasing trend for all of the three atmospheric N 2 O inversions still holds (Mann-Kendall trend test, P < 0.001).
To understand mechanisms underlying this enhanced α PRE and the contribution of precipitation to N 2 O variations, we calculated the difference of these quantities between the most recent and the earliest time windows (denoted as ∆α PRE and ∆N 2 OEA_PRE STD ) by region. Hyperarid, arid, semi-arid, dry subhumid and humid regions all show positive mean ∆α PRE and ∆N 2 OEA_PRE STD across grid cells (figure 3) despite large variations in ∆α PRE and ∆N 2 OEA_PRE STD locally (supplementary figure 7). We further separated the global grid cells into two groups: one group with enhanced variations in precipitation anomalies (∆PRE STD >0, STD stands for standard deviation, ∆ indicates the difference between the recent vs. the initial time window) and the other with reduced anomalies in precipitation (∆PRE STD <0). In arid, semi-arid and dry subhumid regions, the grid cells that experienced intensified anomalies of precipitation show a higher mean ∆α PRE than those with weakened precipitation anomalies ( figure 3(b)). In humid regions, although the mean ∆α PRE is smaller, the mean ∆N 2 OEA_PRE STD is nine times greater in grid cells with ∆PRE STD > 0 than those with ∆PRE STD < 0 ( figure 3(e)). We also found that the standard deviation of the anomalies of the globally integrated precipitation over land increased over time during our study period (supplementary figures 6 (d-f)). This increased precipitation variability, driven largely by global warming [10], is one of the drivers enhancing the response of land N 2 O emission to precipitation.
In addition, we found that the mean ∆N 2 OEA_PRE STD is 2 (arid), 6 (semi arid), 4 (dry subhumid) and 3 (humid) times higher in grid cells with croplands (see Methods) compared to those without within each aridity class ( figure 3(f): supplementary figure 5). With intensive nitrogen fertilizer additions, croplands generally have higher nitrogen substrate that is susceptible to both gaseous and leaching losses. Regions with high emission rates are thus also likely to have high emission anomalies. We then estimated ∆N 2 OEA_PRE STD as a percentage change relative to N 2 OEA_PRE STD from the earliest time window. We still found a higher mean ∆N 2 OEA_PRE STD from croplands than non-croplands in arid, semi-arid, dry subhumid and humid regions (supplementary figure 8). This disproportionately higher ∆N 2 OEA_PRE STD from croplands (supplementary figure 5) is most likely linked with more nitrogen surplus from croplands than natural lands. Furthermore, broadly speaking, places with higher cropland fractions or higher PRE STD × Nfert (precipitation anomaly multiplying total mineral nitrogen fertilizer per grid cell) on average have higher α PRE and especially N 2 OEA_PRE STD (supplementary table 2), which indicates that nitrogen fertilizer additions interact with precipitation in shaping the pattern of N 2 OEA.
Next we investigated whether process-based land N 2 O models capture the top-down N 2 O-precipitation relationship. To tackle this issue, we calculated α PRE from three global land N 2 O models (ORCHIDEE-CNP, ORCHIDEE and OCN) that participated in NMIP and quantified land N 2 O emissions with the same temporal resolution as the inversion estimates (Methods). These models however, are bottom-up calculations on land fluxes that are different from observed atmospheric growth rate anomalies. The correlations between globally integrated land N 2 OEA and PRE are weaker from land models (r = 0.32-0.38, P < 0.05) than from top-down inversions (supplementary figure 9). α PRE increase with time in 2 out of 3 models, but the magnitude of α PRE and its changes are more than three times smaller than from top-down inversions ( figure 1 vs. supplementary figure 9). Compared to the globally integrated N 2 OEA, agricultural N 2 OEA shows a stronger correlation with PRE (r = 0.42-0.56, P < 0.05) and higher α PRE (in unit of per area) with the overall increasing trend despite large differences among the three models (supplementary figure 10). The factorial model simulations from NMIP enable us to separate the effect of nitrogen additions by comparing simulations with vs. without external nitrogen inputs (e.g. nitrogen fertilizer, manure and atmospheric nitrogen deposition). Enhanced N 2 O emissions from agricultural land are driven largely by external nitrogen inputs, especially nitrogen fertilizer additions [13,29]. This is supported by process-based model simulations that show the important contribution of agricultural land and therefore anthropogenic nitrogen additions to the intensified N 2 O-precipitation relationship. Both process-based models and inverse-estimations show the same tendency for N 2 O-precipitation relationship.
The lower sensitivity of N 2 OEA to PRE from NMIP models indicates uncertainties in modelling the impacts of precipitation on N 2 O emissions, which may stem from difficulties in accurately representing the hydrological processes (from precipitation to soil moisture and water filled pore space) [30], or the soil moisture responses of N 2 O emissions. Process-based models generally represent soil moisture effects on N 2 O emissions through empirical soil moisture (or water filled pore space) response functions [3,5,13,31]. These functions are typically derived from limited site-level studies, and may not be universal or adequate in capturing the multifaceted soil moisture triggered changes, for example, the rapid microbial responses following rewetting of dry soils [32]. Due to the complexity and relatively high computation requirements, these global N 2 O models typically do not simulate nitrifying and denitrifying microbial communities explicitly. Whether the lack of microbial physiology leads to biased simulation of moisture responses remains to be explored. Note there are limited numbers of process-based global N 2 O models [29], and a comprehensive assessment on the capacity to model soil moisture responses requires more effort from both modelers and data collectors. The average correlation between satellite derived soil moisture anomalies (ESACCI v05.2, the combined product) and N 2 OEA from atmospheric inversions is 0.27. Satellite derived soil moisture datasets normally only capture first top centimetres of soil, which is not adequate in representing biogeochemical active soil zones. This correlation value reached 0.42 from model estimated soil moisture (CPC soil moisture, one-layer hydrological model) or 0.40 from a reanalysis dataset (ERA5, 0-28 cm). The global scale soil moisture datasets are likely to be associated with higher uncertainties than observation-derived precipitation datasets. While water filled pore space (approximated by soil moisture in models) [5,13] is a more direct regulator of soil N 2 O emissions than precipitation, how to realistically capture its dynamics and impacts on large-scale soil N 2 O emissions requires more efforts in the future.
The top-down (atmospheric inversion) and process-based modelling estimates are consistent in terms of the decadal global N 2 O budget [14], while the climate sensitivities of land N 2 O emissions may need more studies considering the uncertainties in both atmospheric inversions and process-based modelling, and regional scale observations related to this topic are especially in need. N 2 O emissions from inland water bodies are not included in NMIP models, but are incorporated in the top-down inversion estimates. We did not remove N 2 O emissions from inland water bodies in the top-down estimates due to lack of reliable observations. Estimations of N 2 O emissions from inland waters and estuaries remain largely uncertain (Rivers: 0.03-2; Estuaries: 0.06-5.7 Tg N 2 O-N yr −1 ) [33,34], with one recent modelling study arguing for an overestimation of previous values (Rivers: 0.05; Estuaries: 0.06-0.15 Tg N 2 O-N yr −1 ) [34]. Nitrogen lost through leaching is one of the major sources that contribute to aquatic N 2 O emissions [33]. Leaching losses of agricultural nitrogen inputs are closely linked with precipitation. ORCHIDEE-CNP and ORCHIDEE show increased sensitivity of anomalies of nitrogen leaching from agricultural nitrogen additions to PRE (supplementary figure 11), while simulations without agricultural nitrogen additions do not show increasing trends despite the strong correlation between anomalies in nitrogen leaching and PRE (supplementary figure  12, r = 0.84, ORCHIDEE-CNP; r = 0.58, ORCHIDEE). Agricultural N 2 O additions that end up in inland water systems may also play a role in the increasing trend of α PRE and N 2 O PRE from atmospheric inversions, although the magnitude is likely to be small considering the relatively small N 2 O fluxes from inland waters.
The lack of temperature signal on the global scale N 2 O emission anomalies is partly due to spatial aggregation (supplementary figure 13). At the local scale, from gridded inversions fluxes of N 2 O, anomalies in temperature and precipitation play approximately equal roles in N 2 OEA. The relative importances of temperature and radiation decrease at large spatial scales, while the importance of precipitation remains high. Locally, temperature (or radiation) could have both positive and negative impacts on N 2 O emissions (supplementary figure 14). For example, high temperature stimulates nitrifiers and denitrifiers that produce N 2 O, while exacerbated soil drying under warming decreases N 2 O emissions [35] despite the effect also depends on the interaction with precipitation. Warming reduces substrate availability for nitrification and denitrification through enhancing plant growth and nitrogen uptake, but also boosts nitrogen supply through increased nitrogen mineralization [35]. These positive and negative effects, typically involving complex interactions among spatially diverse climatic, biological, edaphic and anthropogenic factors, compensate in space and result in a small temperature (or radiation) signal at the global scale. In addition, locations with high temperature anomalies are from the northern high latitudes (especially Siberia). In Siberia, the N 2 O emission rate is relatively low and inverse estimations are largely uncertain due to very sparse observations. Note here we did not quantify precipitation-temperature interactions, which might be important locally.
Global N 2 O emissions have been accelerated over the last two decades, largely driven by anthropogenic activities that intensify the global nitrogen cycle [14,19]. Interactions between the expedited nitrogen fertilization and intensified precipitation anomalies (associated with global warming) shift land N 2 O emissions. We predict land N 2 O emissions to be more sensitive to precipitation anomalies in the future under the business-as-usual fertilization practice and high emission climate scenarios. To support 10 billion people in 2050, projected nitrogen fertilizer use will have to increase to 141 Tg N yr −1 (from 107 TgN yr −1 in 2015) under the business-as-usual practice [1]. The standard deviation of PRE is 42%, 36% and 26% higher by the end of 21st century under the SSP5-85 (Shared Socioeconomic Pathways with a radiative forcing of 8.5 Wm −2 in 2100), SSP3-70 (radiative forcing: 7.0 Wm −2 ) and SSP2-45 (radiative forcing: 4.5 Wm −2 ), respectively, from an ensemble of Earth system simulations (supplementary figure 15). Only under the sustainable pathway, SSP1-26 (radiative forcing: 2.6 Wm −2 ), the standard deviation of PRE is not increased by the end of this century. The tendency of global aridification and the expansion of the world's dryland under future warming climate [36,37] also indicate a higher sensitivity of future land N 2 O emissions to precipitation. We call for special attention on N 2 O emissions in farming of the semi-arid land, which could contribute to disproportionately high N 2 O emissions (per unit area) under its precipitation variabilities.
Despite the robustness of the overall global trends, our understanding of the mechanism(s) driving enhancements of α PRE and N 2 OEA_PRE is still limited. We demonstrated the intensifying agricultural nitrogen inputs as an important mechanism, but precise data on large-scale mineral fertilizer and manure addition (especially the timing of applications) are lacking. Coordinated efforts on monitoring agricultural nitrogen inputs are in great need. Many factors and types of disturbance, such as crop harvesting, soil tillage and earthworm activity, could trigger local or temporal nitrogen surplus and N 2 O emission anomalies. Whether these factors are regionally or globally important remain to be studied. Here as a first step, we revealed the general relationship between land N 2 O emissions and precipitation anomalies, and focused on the most critical factors at the global scale, while acknowledging many questions (e.g. differences among countries/biomes/land uses) in this subject remain to be answered. We need more attention and community efforts dedicated to monitor land N 2 O emission (with its covariates) and to improve process-based understanding of N 2 O dynamics.

Conclusions
Through the top-down estimates of land N 2 O emissions from three independent atmospheric inversion systems, we found a significant correlation (P < 0.001) between the interannual anomalies in the annual growth rate of the global land N 2 O emissions and precipitation anomalies. The average sensitivity (α PRE ) was 2.50 ± 0.98 Tg N 2 O-N per 100 mm of precipitation across the global land between 1998 and 2016. α PRE increased with time during during 1998-2016, likely due to enhanced precipitation anomalies interacting with increased nitrogen inputs to agricultural lands. α PRE generally increases with aridity and under future climate conditions (with radiative forcing levels of 4.5, 7.0 and 8.5 Wm −2 by the year 2100).

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI [38]: https:// figshare.com/s/2deaae350b24a7c1bcd8.

Code availability
Calculations were conducted through Python 3.7.3, R 3.6.2 and ferret 6.72. Data processing code and code used to generate figures are provided at https://figshare.com/s/2deaae350b24a7c1bcd8.