Integrating Agricultural Practices into the TRIPLEX-GHG Model v2.0 for Simulating Global Cropland Nitrous Oxide Emissions: Model Development and Evaluation

1.Department of Biology Sciences, University of Quebec at Montreal, C.P. 8888, Succ. Center-Ville, Montreal H3C 3P8, Canada 2. Center for Ecological Forecasting and Global Change, College of Forestry, Northwest A&F University, Yangling, China 3. Institute of Wetland Research, Chinese Academy of Forestry, Beijing, China 4. College of Hydrology and Water Resources, Hohai University, Nanjing, China 10


25
Nitrous oxide (N2O) is a long-lived trace gas that has a global warming potential on a 100-year time horizon that is 265-298 times larger than that of carbon dioxide (CO2), and it simultaneously results in ozone depletion in the stratosphere (Ciais et al., 2014). The atmospheric concentration of N2O has significantly increased (i.e., by 20%) since the industrial revolution (Tian et al., 2016;Tian et al., 2020). Generally, N2O is produced as an intermediate product of soil microbial nitrification and denitrification processes and is regulated by multiple biotic (i.e., vegetation type, microbial biomass) and abiotic factors (i.e., 30 climate, soil temperature, humidity, nutrient content, and texture) (Bouwman et al., 2002;Stehfest and Bouwman, 2006;Li et al., 2000;Butterbach-Bahl et al., 2013;Tian et al., 2018).
Cropland is a hotspot of terrestrial N2O sources (Tian et al., 2019;Tian et al., 2020). The current larger emission rate of cropland soil comparing with natural soil (Davidson and Kanter, 2014) results from extensive agricultural practices, including N-fertilizer input (synthetic and manure) (Davidson, 2009;Zhou et al., 2017), irrigation (Li et al., 2000;Li et al., 2010), and 35 tillage (Powlson et al., 2014;Mei et al., 2018), because these agricultural practices directly and indirectly interfere with soil N flow and microbial activities (Cavigelli et al., 2012). Therefore, substantial observation studies have been conducted in croplands to understand the effects of different agricultural practices on N2O emissions in order to enable sustainable agricultural production (Carlson et al., 2017;Burney et al., 2010;Snyder et al., 2009). However, because of the characteristics As a trace N gas, the N2O emitted by nitrification and denitrification were simulated according to the anaerobic balloon concept. The anaerobic volumetric fraction (ANVF) was the key parameter, which represents the soil oxygen status and regulates the allocation rates of the substrates (e.g., dissolved soil organic carbon (DOC), NH4 + , and NO3 -) for nitrification and 120 denitrification. It was calculated using the oxygen partial pressure and the air-filled porosity listed in Supporting Information Table S1 (Equations (1-3)).
In the TRIPLEX-GHG model, nitrification is an aerobic process that occurs outside of the anaerobic balloon, converting ammonium (NH4 + ) into nitrate (NO3 -) driven by nitrifying bacteria with N2O as a by-product (Li et al., 2000;Zhang et al., 2017b;Morkved et al., 2007). The growth and death rate of nitrifiers  in Table S1) are highly dependent on 125 the DOC; therefore, the nitrification rate (Equations (7-10) in Table S1) was calculated using the Michaelis-Menten function based on the concentration of NH4 + and the microbial activity of the nitrifying bacteria. The effects of the soil properties were also simulated  in Table S1).
Denitrification is the process through which the nitrate is reduced stepwise into different nitrogen gases as a chain reaction process inside of the anaerobic balloon. Denitrification can be divided into 4 independent steps, which are linked by the 130 competition for DOC by the specific denitrifiers during each step (Betlach and Tiedje, 1981). Similarly, the growth and mortality rates of the different denitrifiers utilized a double substrate based (DOC and NOx) Michaelis-Menten equation  in Table S1). The consumption of for the growth of the different denitrifiers was calculated at an hourly time step according to previous studies as is shown in the following Eq. (1) (Leffelaar and Wessel, 1998;Li et al., 2000): (1)

N output
Except for gaseous N losses (N2O, NO, NH3), crop uptake/harvest processes, leaching, and surface runoff are the major N 150 outputs (Liu et al., 2010). We altered the plant N uptake process and integrated harvest practices into the original model.
Plant N uptake: As a plant grows, mineral N is taken up as NO3 --N and NH4 + -N, which is considered to be the dominant pathway of soil N loss (Sebilo et al., 2013). It has been reported that NO3 --N is much more easily absorbed by roots (Malhi et al., 1988;Kronzucker et al., 1997), so NO3 --N was set as having a higher priority of being taken up by plant roots in each soil layer using the following Eq (2) and Eq (3): Here, _ 3 − and _ 4 + are the plants' NO3and NH4 + requirements from soil layer ; and _ is the total plant N uptake requirement from soil layer . 3 is the coefficient of nitrate demand, which was set to 4.0 according to the model test. The comparison between the tendency of the modeled soil nitrate content and the 160 detected soil nitrate concentration proved the effectiveness of the model design (Fig. S1).
Harvest: Harvest practices significantly reduce the soil C and N inputs for cropland compared with natural soil. We systematically removed all of the litterfall from the cropland ecosystem at the end of the growing seasons to modify the harvest.

N Input
The input N flows of agricultural soil include N fertilizer, biological N fixation, atmospheric deposition, and returned straw (Liu et al., 2010). We integrated chemical N fertilization, manure application, and returned straw processes into the original model.

165
Returned residues: Returned residues are a significant C and N source for cropland soil and a recommended practice for improving N use efficiency. It has been reported that more than half of the N in crops (all of which is taken up from the soil) is removed from the ecosystem, and only 20% of the crop's N is returned to the soil N pools as returned residues globally (Liu 170 et al., 2010). Therefore, in our model, we returned 20% of the harvested plants to the cropland soil in order to modify the returning residue practices.
Chemical N fertilizer: Chemical N fertilizers were directly added to the NO3and NH4 + pools of the top layer of soil in the original model.
Manure N fertilizer: The manure-sourced N entered the different inorganic N and organic N pools separately. The organic 175 portion of the manure was added to up to 3 soil organic matter (hereafter SOM) pools (the non-protected, protected, and passive organic carbon pools) separately for further decomposition (Zhang et al., 2017b).

180
Here, 4 + and 3 − are manure-sourced NH4 + -N and NO3 --N, respectively, which are calculated using the ratio of ammonia and nitrate (i.e., 4 and 3 ) to total manure N. is the amount of manure that entered the different SOM pools; is the proportion of manure N added to the different SOM pools; and : is the C:N ratio of a particular SOM pool.

Irrigation and tillage
185 Some agricultural practices do not directly alter the N flow in cropland soil, instead they affect N2O by regulating the soil's physical properties.
Irrigation: The irrigation process used in this study adopted the idea of precipitation events from the DNDC (Li et al., 2000) and Agricultural Production Systems sIMulator (APSIM) (Thorburn et al., 2010). Similar to rainfall, irrigation provides extra water to the surface of the cropland soil, promoting the water-filled pore space (hereafter WFPS), thus stimulating the growth 190 of the anaerobic balloon. Nevertheless, it simultaneously induces leaching and runoff processes. In the current model, only the flood irrigation method was included.
Tillage: Tillage redistributes the soil profile and increases the availability of oxygen in each soil layer at the same time. We averaged each of the C and N pools of the top 3 soil layers (as a global conventional tillage depth) after every tillage event.
Because of more exposure to oxygen, the anaerobic conditions and diffusion pattern also vary with the different soil moisture 195 conditions, properties, and vegetation types (Rochette, 2008;van Kessel et al., 2013).

Model sensitivity analysis
We conducted initial sensitivity analysis experiments to obtain the most sensitive parameters before testing the model.

200
According to previous N2O modeling studies (Zhang et al., 2017b;Zhang et al., 2019), the coefficient of nitrification (hereafter COENR) is the key parameter driving the amount of emitted N2O in natural ecosystems probably because of the limited NO3input. In this study, considering the increased NO3input from fertilizers in cropland soil, it is conceivable for denitrification to become the dominant N2O source. Therefore, 13 major parameters, including COENR and parameters associated with denitrification (Table1), were compared in a site-specific manner. The sensitivity index (SI) in this study followed the method

205
of Lenhart et al. (2002) using the following Eq(7): where n is the total number of months from 1961 to 2015 (because in our model, chemical fertilizer application started in 1961); j accounts for the number of months from 1961 to 2015 (because in our model, chemical fertilizers were used after 1961); 0 represents the jth monthly N2O emissions with an initial parameter 0 ; and 2 1 are the N2O emission

Input data
All of the input information for the model simulation of the selected sites described above was directly obtained from the following datasets or was obtained from papers (see details below). These data were transformed into a spatial resolution of 0.5°×0.5° latitude/longitude using the ArcMap software (version 10.2) before the simulation. average, and maximum temperature, precipitation, specific humidity, air pressure, and wind speed, which were used to drive the model.

N fertilization data:
The historical chemical fertilizer  and manure (1860-2014) application data for croplands were derived from the datasets produced by Nishina et al. (2017) and Zhang et al. (2017), respectively.

230
The synthetic N fertilization dataset is mostly based on country-specific information from the Food and Agriculture Organization statistics (FAOSTAT) after filling data gaps (Nishina et al., 2017). Notably, the dataset provided application date and monthly input N fertilizer differentiated into NH4 + and NO3 − considering the seasonal crop calendars for the dominant crops in each grid (Sacks et al., 2010). The synthetic N application rates in 2011-2015 were assumed to be the same as that for 2010. In addition, if the amount and type of N fertilizer from Nishina et al. (2017) failed to match the site information 235 obtained from the literature, we utilized the site-specific amount of fertilizer application according to the published paper and the NH4 + -N/ NO3 --N ratio provided by Nishina et al. (2017).
The manure N dataset (Zhang et al., 2017a) included the annual manure production and annual application, which were reconstructed using the dataset from the Global Livestock Impact Mapping System (GLIMS) in conjunction with countryspecific annual livestock populations and the gridded cropland distribution map for 1860-2014 obtained from HYDE 3.2 240 (Goldewijk et al., 2017). The manure N production and application rates in 2015 were assumed to be the same as those in 2014.  (Lelieveld and Dentener, 2000), which used N emission estimates (van Aardenne et al., 2001) and projection scenario data (Houghton, 1996;Nakicenovic et al. 2000). Environment, version 3.2 (HYDE 3.2), which has reconstructed time-dependent land use using historical population and allocation algorithms with weighting maps (Goldewijk et al., 2017). Cropland can be classified into rain-fed and irrigated land, both of which were further divided into rice, generic C3 crops (except rice, e.g., wheat), and generic C4 crops (e.g., maize) based on the global crop distribution maps (Monfreda et al., 2008).
conditions, which was the multiyear mean climate data.
For the model calibration, we used the daily climate data for each site to drive the model along with other site-specific input information. The simulation started on January 1 st , 1901, and ended on December 31 st , 2015, with a daily time step. By

275
comparing the output N2O flux data with the observed data, we adjusted the most sensitive parameter of the N2O emissions based on the sensitivity analysis in order to fit the best model performance via trial and error and statistical model performance indicators. The index of agreement (D), the root mean square error (RMSE), and the coefficient of determination (R 2 ) were used to evaluate our model's performance, and the D-value and RMSE were calculated as follows: Here, is the ith simulated result corresponding to the number of observations; is the ith observed value; and is the mean of the observed values during the experimental period. D varies between 0 and 1, and is excessively sensitive to extreme values (Willmott, 1981). The model performance was considered to be perfect and unmeaningful when the D value was set to 1 and 0, respectively. The RMSE is the key value representing the difference between the simulated and observed values, and 285 is significantly affected by the data units.
Based on the calibration results and the fitting of the most sensitive parameter for the different sites, we used the continental mean parameter for the model validation.

Sensitivity analysis
The mean sensitivity index (SI) varied from -0.53 (EFFNO2) to 1.37 (COEdNO3) for the selected 13 parameters (Fig. 2). All of the parameters had a nonunique effect on the N2O emissions of the different sites. COEdNO3, COENR, MUENO3, MNO3, EFFN2O, COEdNO2, and COEdNO mostly had positive effects, while the remaining parameters either had negative effects (e.g., MUEN2O and EFFNO2) or had no evident impact (e.g., AMAX) on the N2O fluxes. The coefficient of the NO3consumption rate (COEdNO3)

295
was the most sensitive parameter in the current TRIPLEX-GHG model. The SI ranged from -0.61 to 5.39 (with a mean of 1.37) for the current model input information. We also noticed that the SIs of the selected parameters were not consistent with the different input information, especially for the variations in the amount of N fertilizer applied. The COEdNO3 slightly increased initially and then decreased as the N dose increased; and as the most sensitive parameter, it retained a large SI value (Fig. S2).
Overall, to simplify the parameter fitting processes and to evaluate the model's performance, we selected COEdNO3 as the 300 fitting parameter, while we set the other parameters to their original constant values as the default (Table 1).  (Table   3).

North American sites
The data collected for the North American cropland sites were located in the US and Canada and represented the dominant commercial crop species such as corn, wheat, barley, and tomatoes. Most of the measurements were collected over more than

Asia
Ten upland agricultural sites were selected in Asia (Table 2 and

Europe
Most of the wide-spread crop types were included in the calibration of the model simulation of the European cropland sites, which were located in the mid-high latitude region. The simulated trends and magnitudes of N2O were generally consistent with the measured data for most of the sites, but some of them had relatively low agreement indices. Based on the

South America & Africa
Unfortunately, there are insufficient observations of cropland N2O emissions conducted in the agriculturally dominated 380 regions of South America and Africa ( Fig. 7 and Table 2). Typical agricultural ecosystems in these regions (sugarcane, wheat) were selected for the model calibration. Compared with the results of the two sites with short experimental periods in Africa In summary, according to the calibration results, the trends and magnitudes of the simulated N2O flux were generally consistent with the measured field data.

395
The model validation (Fig. 8)  of the validated sites, and the results are also presented in Table S2. During the validation, the simulated daily mean emission rates during the experimental periods ranged from 0.048 to 5.21 mg N m -2 day -1 , and most of the values were less than 1 mg N m -2 day -1 . The regression result was close to the 1:1 line, indicating that the modeled results are quite consistent with the observed N2O emissions (R 2 = 0.86, p<0.001). However, the modeled results tend to slightly underestimate the N2O flux for 400 the low observation values (<1 mgN m -2 day -1 ) and to overestimate for large observed N2O flux values (>1 mgN m -2 day -1 ).
The model validation results further confirm that our model is capable of simulating the impacts of both climate and agricultural practices on N2O emissions across global cropland ecosystems.

405
It is important to calibrate the process-based model using reasonable parameters in order to simulate complex biogeochemical processes better. Adjusting the most sensitive parameter is an efficient method of improving the model performance and has been widely used in model development and parameterization (Wang and Chen 2012;Zhang et al., 2017b;Zhu et al., 2014). As was stated in a previous study, Zhang et al. (2017) tested 23 parameters and found that the COENR, the coefficient of the nitrification rate, controlled the N2O emission process. Meanwhile, our sensitivity analysis results revealed 410 that the coefficient of the nitrate consumption rate, COEdNO3, had the highest sensitivity level for the updated version of the TRIPLEX-GHG model. Such a divergence is probably due to the increased N input, especially NO3 -, in cropland ecosystems compared with natural grasslands and forests. Denitrification strongly contributes to N2O production in agricultural ecosystems, which requires NO3as a substrate (Wang et al., 2018a). Since it is controlled by this parameter, the NO3consumption (from NO3to NO2 -) dominates the denitrification rate and thus the N2O production rate of N fertilized soil. Globally, the COEdNO3 415 exhibited a large range of variation during the parameterization, which can partly be reconciled by the calibration method and the varying amounts of mineral N input. Because the NO3consumption rate for denitrification is difficult to measure directly, the limited field information strongly discourages the systematic adjustment of the COEdNO3, and thus, the potential uncertainty of the parameter affected the model's performance.
Generally, the TRIPLEX-GHG model reproduces the N2O emissions well for a daily time step and various cropland 420 ecosystems (e.g., wheat, maize, sugarcane, and cotton) on a global scale. The dominant characteristic of cropland N2O emission is the peaks associated with fertilization events, most of which were well simulated by our model and contributed to overall reasonable evaluation indices. Such advantages were derived from three features of our model. First, both the soil oxygen conditions and the soil water conditions were considered in the TRIPLEX-GHG model, i.e., represented by the size of the anaerobic balloon and the water-filled pore space, respectively. Previous studies have highlighted that the soil O2 status is the 425 proximal, direct, and most decisive environmental trigger of N2O production (Song et al., 2019;Zhu et al., 2013;Khalil et al., 2004). However, the majority of process-based models only integrated the WFPS into the nitrification and denitrification processes (e.g., Tian et al., 2010;Ito et al., 2018). It was reported that although the WFPS is a critical element containing information about the soil water and gaseous status, it still requires combination with other soil structural parameters in order to better predict the soil O2 concentration, microbial respiration, and subsequent gas diffusion (Farquharson and Baldock, 2008;430 Song et al., 2019;Hall et al., 2013;Rabot et al., 2015). Second, a detailed description of the manure also contributed to the improved model performance because manure is a predominant soil organic carbon (SOC) source for croplands, which is not considered by empirical models and several of the process-based models (e.g., DAYCENT, VISIT). The SOC serves as a key energy and carbon source for microbial growth, nitrification, and denitrification (Snyder et al., 2009;Butterbach-Bahl et al., 2013). Field observations have shown that the application of manure either promotes or reduces N2O emissions probably 435 because the added organic C compounds support microbial growth, but the increased SOC stimulates complete denitrification with the further reduction of N2O to N2 (Zhou et al., 2017;Meijide et al., 2007). Therefore, in this study, the manure sourced C was recalculated using the manure N and C:N ratio, which significantly enhanced the simulation of the SOC. Last but not least, the TRIPLEX-GHG model included a reasonable microbial growth and death description, which strongly improved the accurate modeling of the nitrification process because the soil microbial conditions are one of the primary determinants of the 440 soil nitrification rate at a global scale .
However, there are still major discrepancies between the modeled and measured N2O fluxes, including underestimated peak values, failure to capture emission peaks, and underestimated background emissions. First, although the timing of the simulated major emission pluses was well simulated, the peak values of the emitted N2O fluxes were underestimated. The incomplete description of the processes involving the interaction between the soil pH and the external mineral N input is 445 probably responsible for this phenomenon. The soil pH is one of the most important drivers of N2O production. Acidic soils are more sensitive to N input than alkaline soils, which probably enhances N2O production in croplands (Wang et al., 2018b;Morkved et al., 2007). Studies have shown that the pH values of agricultural soil tend to be significantly reduced by N deposition and N-fertilization at the global scale (Tian and Niu, 2015;Godsey et al., 2007;Guo et al., 2010). However, because the soil buffer capacity is difficult to quantify (Baron et al., 2014;Zhang et al., 2017b), the soil pH in our model was input  (Teepe et al., 2004;Groffman et al., 2006). The latter triggers microbial driven nitrification and denitrification processes (Sharma et al., 2006). The limited description of those processes, especially the simple empirical parameters and algorithms 460 we used for modeling snow-melting hydrology and nutrient release, are the primary error sources (Zhang et al., 2017b).
Furthermore, as for the underestimated background emissions, it is still a significant challenge for the process-based model to accurately quantify background N2O emissions due to the following possible reasons. First, our simulations used general crop classification (C3, C4, and rice) instead of detailed crop rotation information with different physiological parameters (Ito et al., 2018;Monfreda et al., 2008;Saikawa et al., 2013). Field observations have revealed that different crop 465 types or species have diverse impacts on the N2O fluxes of cropland (Rochette et al., 2018;Philibert et al., 2013;Petersen et al., 2006;Gelfand et al., 2016). For instance, legume species (e.g., soybean) have a stronger N fixation ability, which contributes considerably to the N pools in cropland soil (Liu et al., 2010), and they effectively promote background N2O emission even without N fertilization compared with other cereal crops (Lenka et al., 2017;Sanchez and Minamisawa, 2019;Yang and Cai, 2005). Second, in addition to climate conditions, the background emission rates from agricultural soils are also because these agricultural practices are controlled by the individual farmers and vary greatly at the local and subregional scales, without clear global distribution patterns such as those for soil and climate (Wang et al., 2018b). Third, the uncertainties in the 475 site history are also responsible for the inaccuracy of the modeled background emissions because the site history has a tremendous effect on the soil properties, especially the SOC content (Gelfand et al., 2016). Previous studies have demonstrated that agricultural practices, such as returning residues to the soil, tillage management, and fertilizer application, are important drivers of the SOC (Liu et al., 2014;Jiang et al., 2018;Zhou et al., 2017), but their effects vary with the intensity of the practices and the climatic conditions (Ogle et al., 2019;Snyder et al., 2009;van Kessel et al., 2013;Liu et al., 2014;Gattinger 480 et al., 2012). Unfortunately, only a few published papers have provided detailed historical land use and agricultural practice information, which is a barrier to the accurate estimation of the local SOC and thus N2O emissions.
Last but not least, the other reasons for the discrepancies between the modeled results and the observations may be the uncertainties in the field measurements and the driving data. For one, the lower sampling frequency of the fieldwork and the short-lived N2O emission pluses are particularly difficult to captured with traditional manual chambers, especially after base 485 fertilizer application in the fallow season (Lammirato et al., 2018;Lognoul et al., 2019). Moreover, the model's accuracy also relies on good quality data. A 0.5°×0.5° global scale daily climate input dataset was used for the model calibration and validation, but it is unlikely that every site was provided with detailed meteorological information due to the relatively coarse spatial resolution. The local climate may differ significantly from that of the grid input information (Wania et al., 2010). Specifically, the precipitation information is less accurate compared with the other climate data, which could significantly 490 jeopardize the model's performance since the anaerobic balloon is a precipitation-induced process (Zhang et al., 2017b).
Furthermore, the soil properties are also difficult to be precisely replicate at the site level using a global soil dataset. Because the soil texture served as a significant driver for the N2O emissions (Philibert et al., 2013;Gu et al., 2013), the mismatch of soil information is also a major cause of the disagreement between the model simulations and observations.
In response to the uncertainties described above, further modeling is suggested to improve the detail of the descriptions 495 of the key processes, and better quality datasets need to be collected.
First, it is recommended that more detailed soil microbial activities be considered in order to better model the features of the N2O emissions . The nitrifier-denitrification process may account for up to 100% of the N2O emissions from NH4 + in soils (Wrage et al., 2001;Wrage-Moennig et al., 2018), especially for N2O uptake, the occurrence of which has been widely observed in peatlands, boreal forests, (Saikawa et al., 2013) and occasionally in cropland ecosystems (e.g., Fig.   500 4b). In addition, ammonia oxidation has been found to be a significant process for the development of N2O compared to classical denitrification in extremely low-oxygen concentration soils (Zhu et al., 2013).
Next, in this study, only general PFTs were used for croplands without specifying crop types, for which the nutrient requirements, maximum productivities, C:N ratios of different organs, and biomass allocation patterns differ significantly from each other (Li et al., 2000;Shan and Yan, 2013), which significantly affects the N dynamics of cropland soils.

505
Furthermore, various agricultural management practices were not included in the current model. For example, different techniques of tillage (e.g. conventional tillage, minimum tillage, tillage with different instruments), irrigation (e.g., flood irrigation and drip irrigation), and fertilizer placement (e.g., top dressing and injection) can have diverse impacts on N2O emissions (Maris et al., 2015;Rochette, 2008). For instance, drip irrigation effectively promotes WFPS without surface runoff, which induces significant N2O flux (Sanchez-Martin et al., 2008). Considering the advantages of field studies, model 510 performances can be effectively improved at the site level.

Conclusions
Our study represents a successful attempt to fully integrate general agricultural activities into the current TRIPLEX-GHG framework for simulating global N2O emissions across cropland ecosystems. In this study, the COEdNO3, which controls the 515 NO3consumption rate of the denitrification process, was found to be the most sensitive parameter.