Modeling biological nitrogen fixation in global natural terrestrial ecosystems

Biological nitrogen fixation plays an important role in the global nitrogen cycle. However, the fixation rate has been usually measured or estimated at a particular observational site. To quantify the fixation amount at the global scale, process-based models are needed. This study develops a biological nitrogen fixation model to quantitatively estimate the nitrogen fixation rate by plants in a natural environment. The revised nitrogen module better simulates the nitrogen cycle in comparison with our previous model that has not considered the fixation effects. The new model estimates that tropical forests have the highest fixation rate among all ecosystem types, which decreases from the Equator to the polar region. The estimated nitrogen fixation in global terrestrial ecosystems is 61.5 Tg N yr−1 with a range of 19.8–107.9 Tg N yr−1 in the 1990s. Our estimates are relatively low compared to some early estimates using empirical approaches but comparable to more recent estimates that involve more detailed processes in their modeling. Furthermore, the contribution of nitrogen made by biological nitrogen fixation depends on ecosystem type and climatic conditions. This study highlights that there are relatively large effects of biological nitrogen fixation on ecosystem nitrogen cycling. and the large uncertainty of the estimation calls for more comprehensive understanding of biological nitrogen fixation. More direct observational data for different ecosystems are in need to improve future quantification of fixation and its impacts.

N2N fixation by legumes. The model schematic and other calculations including carbon cycle and nitrogen cycle are inherited from an earlier version of TEM (Zhuang et al., 2003;Yu and Zhuang, 2019).
BNF is the most significant process in either symbiotic or non-symbiotic forms converting stable molecular N2 95 into N chemical compounds that are available to plants. For most terrestrial ecosystems, N2N fixers could be in many forms, such as free-living bacteria, lichens, and blue algae. But among them, symbiotic BNF is a dominant process to provide biologically accessible N, and most systematical BNF is regulated by legume plants, especially in croplands and semi-natural environments (Mus et al., 2016). In natural environments, contributions from legumes can be significant but with large uncertainties, which is greatly determined by various environmental conditions 100 (Lindemann and Glover, 1996). In this study, the N2N fixation via legume plants is modeled considering (1) the accessible N concentration in soils, (2) the limitation of temperature, (3) soil water status, (4) the carbon demand for N2N fixation, and (5) the percentage of N2N fixing plants for each ecosystem type as: where is the nitrogen fixation rate, fixpot N is the potential N2N fixation rate (g N day -1 ) , t f is the influence 105 function of soil temperature,  Table   4 for value range of related parameters.
The potential N2N fixation is highly related to the total N demand of plants and the available nitrogen in soils. Theoretically, the definition of potential N2N fixation rate should be the difference between the demand and 110 supply of N. Both of them vary with plant types, stages of growth and soil conditions. For large spatial-scale simulations for various ecosystem types, it is impossible to derive potential N2N fixation because of data availability. fixpot N can be estimated based on dry matter of root, nodule or plant dry matter (Voisin et al, 2003(Voisin et al, , 2007. However, root biomass is also difficult to measure directly. In most published studies, the potential nitrogen fixation rate was measured using an acetylene reduction array (ARA) method (Hardy et al, 1968(Hardy et al, , 1973, and some 115 used 15 N methods (Shearer and Kohl, 1986). In our simulation, fixpot N is assumed to be a constant for each ecosystem type. The fixpot N range is determined from literature and specific values for various ecosystem types are obtained through model parameterization.
Soil temperature is a controlling factor for both microbial activities and plant growth. A large number of studies show that different plants have slightly different preferences for temperature (Montanez et al, 1995; 120 Breitbarth et al., 2007;Gundale et al., 2012). For soybean, 20-35 °C is optimal (Boote et al., 2008), and for white clover the optimal temperature can be 13-26 °C (Wu and McGechan, 1999). The activity of microbes responds slightly differently to temperature among species. For most of them, the optimum temperature is 20-25 °C, and at 12-35 °C the activity is not limited. Generally, the relation between the factor and temperature is not exactly a Gaussian distribution. BNF increases as the temperature rises from minimum temperature (0-5 °C) for N fixation to optimal temperature, maximum rate occurs within an optimal range (15-25°C), and decreases from optimal to maximum temperature above which BNF will stop at 35-40 °C: where the upper limit (tmax) is set to 45 °C. There is no lower limit, but when t is low enough, ft will be close to zero 130 (Wu and McGechan, 1999;Boote et al., 2008;Holzworth et al., 2014) (Table 1). For the convenience in computing, a lower limit is set in our model. When the temperature goes beyond its upper or lower limit, ft is assumed to be 0.
Water stress has a direct effect on nitrogen fixing system (Sprent, 1972). With proper temperature, soil moisture condition is the major factor controlling nitrogen fixation rate (Srivastava and Ambasht, 1994). Soil water deficit and flood dramatically inhibits N2N fixation because of drought stress and oxygen deficit, respectively 135 (Omari et al., 2004;Mario et al., 2007). In our model, the water factor is linearly related with soil water content (Williams, 1990;Wu and McGachan, 1999): where f W (J kg -1 ) is the available soil water, which is defined as the ratio of water content to that at the field capacity. In soils, water potential generally includes osmotic and matrix potentials, ranging from -0.1 to -0.3 bar for 140 typical soils, which has little effects on the N fixation. But when the soil gets very dry, the potential can be up to -100 to -200 bar and increases rapidly. a W is the bottom threshold below which N2N fixation is totally restricted by soil moisture. b W is the upper threshold above which nitrogen fixation is not limited by soil moisture. 1  and 2  are parameters representing the linear relationship between soil water content and its effect on N2N fixation, respectively (Table 1).

145
It is generally thought that more substrate N in soils will slow down the N2N fixation, because plants can uptake N directly from soil with less energy (Vitousek and Field, 1999). By comparison, N2N fixation needs more energy and consumes more carbon than plant N uptake does. Thus, the N2N fixation is only considered to occur when the direct N uptake from soil cannot meet the plant N demand. In our model, the inhibition effect of N is defined as (Wu and McGehan, 1999): Where Nup f is a parameter related to legume biological N2N fixation and soil N. S N is the soil mineral N (g N m -N2N fixers get photosynthetic carbohydrate support from plants. Because the product of every unit of nitrogen fixed consumes a certain amount of carbon, the lack of carbon supply will inhibit the N2N fixation. The carbon cost for per unit of fixed N2N varies widely depending on environmental conditions and ecosystem types. For example, the consumption of carbon is only 1.54 times of fixed N2N for cowpea (Layzell et al., 1979), and it can 160 be 6.3 to 6.8 times for soybeans (Ryle et al., 1979). It is also related to the life cycle of plants. The carbon effect is modeled following a Michaelis-Menten equation (Boote et al., 1998): where Cr is the soil carbon content (g C m -2 ) to represent carbon availability from plants to N2N fixers. Kc is the Michaelis-Menten constant, which is plant species dependent.

Data
The classification of land cover and leguminous biomes were derived from the combination of the International Geosphere and Biosphere (IGP) land-cover classification system and the study of Schrire et al (2005). The experimental N2 fixation data for model calibration were collected for 7 major ecosystem types. Nitrogen fixation rates were determined with acetylene reduction assay (ARA) method in most published studies ( The parameters for N2 fixation module were initialized with a priori values (Table 2). Ecosystem-specific and microbe guild-specific parameters were inherited from previous TEM model (Zhuang et al., 2003;Yu and Zhuang, 2019). The global simulations were conducted at a spatial resolution of 0.5 by 0.5 degree and at a monthly time step.

175
Historical climate data including temperature, precipitation, cloudiness and water vapor pressure were derived from the Climate Research Unit (CRU) (Mitchell and Jones, 2005). Soil texture data were from Melillo et al. (1993) and calculated to show the mean difference between simulated data and observational values. The model is iterated with changing parameters until the RMSE reached a certain value for each site. Most parameters in the model driving nitrogen cycle in the soil have been defined and calibrated in previous studies (Yu and Zhuang, 2019). The 190 calibrated model is evaluated at the site level and then extrapolated to the global terrestrial ecosystems.

Model sensitivity and uncertainty Analysis
The response of N2N fixation of different biomes to input data and variation of parameters was analyzed using sensitivity testing. Four major input variables were selected, including air temperature, precipitation, soil nitrogen

Model evaluation
To evaluate the model, thirty-five observational sites were selected for 7 major ecosystem types across the globe, representing different climate and soil conditions. The experimental data of N2N fixation have a mean value of 12.9 kg N ha -1 yr -1 , with a standard deviation of 17.7 kg N ha -1 yr -1 . The maximum observed fixation occurred in 205 temperate forest in New Zealand, while the minimum rate was also for temperate forest in the state of Idaho State of in the US. Our simulations are comparable with the observed data for all major ecosystem types with the coefficient of determination (R 2 ) of 0.44 and with a slope of 0.46 ( Figure 2). The regression results are mainly influenced by some observed data greater than 30 kg N ha -1 yr -1 . By removing the outliers of observational data, the slope of regression increases to 0.72. Observational data for temperate forests show the greatest variation among all major 210 ecosystem types, with a maximum value reaching 800 times of the minimum one. Simulations are closer to the observations across sites in temperate forests with R 2 of 0.26 and slope of 0.42. Our model underestimated nitrogen fixation rate in temperate forests. The large variation in observations may be due to the distribution of legume plants, different sampling time periods (e.g., growing and non-growing seasons), and varying climate conditions.
For tropical forests, our model estimates of N2N fixation are higher than observations with the slope of 0.75 and R 2 215 of 0.44.

Model sensitivity analysis
The model sensitivity analysis quantifies the impact of changes in forcing data on nitrogen fixation rate.
Climate conditions including air temperature and precipitation, and soil characteristics of nitrogen content and carbon content varied at 3 levels to examine the sensitivity. The response of nitrogen fixation rate emissions is 220 quantified for each ecosystem type. The sensitivity test was conducted for all observational sites (Table 2).
Temperature is the most sensitive variable ( Figure 1). Nitrogen fixation is more sensitive to the change of all forcing conditions. Increasing soil nitrogen results in a lower N2N fixation. Abundant soil nitrogen content inhibits BNF activity, but stimulates nitrification and denitrification processes.

225
Tropical forests in South America, Central Africa and South Asia show a wide range of N2N fixation rates between 1 and 200 kg N ha -1 yr -1 (Bruijnzeel et al, 1991). Here all plants in tropical rainforest are assumed to fix nitrogen and one set of parameters are applied for all tropical forests. The coverage for tropical forests in the landscape was assumed to be 15% (Cleveland et al., 1999), ranging from 5% to 25%. The N2 fixation rate was estimated to be 18.2 kg N ha -1 yr -1 , which is the highest among all vegetation types. Our simulations show that the 230 total fixed nitrogen ranges from 10.8 Tg N yr -1 to 54 Tg N yr -1 , with the average value of 32.5 Tg N yr -1 (Table 3).
Nitrogen fixation in tropical forests is almost half of the global total amount and a principal contributor of BNF in natural ecosystems. Tropical forests have the largest potential to fix nitrogen given that the optimal temperature and soil moisture for BNF is relatively easy to have under tropical climatic conditions. kg N ha -1 yr -1 . The estimates of the total nitrogen fixation were between 1.9 and 19.14 Tg N yr -1 (  (Woodmansee et al., 1981). The coverage of nitrogen fixers was assumed to be from 5% to 25%, 250 (Cleveland et al., 1999). Our simulation assumed that nitrogen fixers cover 15% of the land, resulting in 1.9 kg N ha -1 yr -1 fixation, representing a much smaller fraction compared to forest ecosystems. Total fixed nitrogen in grasslands appeared to range from 0.62 to 3.1 Tg N yr -1 , with an average of 1.86 Tg N yr -1 . For savanna, the total contribution was less due to its relatively small area. The minimum, average and maximum values were estimated to be 0.45, 1.34 and 2.23 Tg N yr -1 , respectively.

255
In tundra and boreal forest regions, both host plants and their rhizobia are adapted to the environment with low temperature. Nitrogen fixation rate is extremely variable for boreal ecosystems. For tundra, the coverage was assumed to be 3-15%, and for boreal forest, the coverage was 4-18%. But in general, the low temperature and permafrost conditions limit the activity of nitrogen fixers (Alexander, 1981). We estimated that tundra ecosystems fix nitrogen at 3.2 kg N ha -1 yr -1 . Their total BNF was between 0.51 to 2.55 Tg N yr -1 with average of 1.54 Tg N yr -1 .

260
In boreal forests, the fixation rate was much lower (2.1 kg N ha -1 yr -1 ) compared to temperate forests.
The fixation could be neglected in deserts because of the extremely dry conditions. Only few legumes may exist in deserts and their growth is highly depended on precipitation events. Even in semi-arid areas, the N2 fixation rate is much lower than that in tropical and temperate forests (5.7 kg N ha-1 yr-1).
In deserts, although some legumes may exist in extremely dry conditions, and some species may grow rapidly 265 after rainfall, their fixation could be neglected. However, in semi-arid areas, legumes are common plants with several species, their N fixation is lower than tropical and temperate forests (5.7 kg N ha -1 yr -1 ).
Mediterranean ecosystems such as in southern California and some areas in southern Australia are characterized with mild rainy winter and hot dry summer, containing both evergreen and deciduous shrublands, in which nodulated legumes are prominent (Sprent et al., 2017). These legumes are more active in comparatively wet season 270 than in dry season (Sá nchez-Diaz, 2001). The ability to fix nitrogen is considered to be one of the most important features that enable legumes and plants to survive under severe environments (Crisp et al., 2004). We estimated that the N2N fixation rate of these legume species is similar to that in grasslands (2.7 kg N ha -1 yr -1 ).
Spatially, the highest rate of N2N fixation occurred in the tropical and sub-tropical areas, as a result of proper climate and soil characteristics for fixers ( Figure 3). N fixation from tropical forests and xeric shrubland contributes 275 to nearly half of the global terrestrial amount ( anthropogenic N input to terrestrial ecosystems is expected to inhibit the natural BNF and might lead to less BNF in the future.

Comparison with other estimates of biological nitrogen fixation (BNF)
There is a large uncertainty in estimating the N input into terrestrial ecosystems, especially from BNF (Sutton et  the mean annual global BNF was estimated to be 89-100 Tg N yr -1 . By assuming a steady state between N input to and loss from ecosystems, Vitousek et al. (2013) estimated the BNF to be 58 Tg N yr -1 with a plausible range of 40 -310 100 Tg N yr -1 , which is similar to our estimates. However, Xu-Ri and Prentice (2017) estimated that the N2N fixation was about 340 Tg N yr -1 which is almost 5 times larger than our estimates. In their study, BNF was determined by plant N requirement across all biome types.
In our estimation, tropical forests significantly contribute to the total BNF, which is up to 18 kg N ha -1 yr -1 . This result is highly related to the density of leguminous plants, and the physical conditions in tropical areas (Crews, 315 1999). Our simulated results are comparable to the estimates of symbiotic N2 fixation from tropical evergreen (5.5-16 kg N ha -1 yr -1 ) and deciduous forests (7.5-30 kg N ha -1 yr -1 ) (Reed et al., 2011). Barron et al. (2010) directly measured N2-fixing root nodules across lowland tropical forests and their observations also showed a large variation among individual trees. For a mature forest matrix, the average value was around 10 kg N ha -1 yr -1 , but it could be as high as 200 kg N ha -1 yr -1 for some areas. Cleveland (2013) provided a similar estimate to ours (around 12 kg N ha -1 320 yr -1 ), but higher values (20-30 kg N ha -1 yr -1 ) in their earlier studies (Cleveland et al., 1999). Sullivan et al (2014) analyzed human's impact on tropical N fixation and found, depending on forest ages, fixation was 5.7 kg N ha -1 yr -1 with a range from 1.2 to 14.4 kg N ha -1 yr -1 , which is lower than our estimates.
For temperate and boreal forests, we estimated that BNF fixation is 2.1-18 kg N ha -1 yr -1 . The existing BNF estimates from literature also show a large uncertainty for those forest ecosystems. For instance, LM3V-N model 325 (Gerber et al., 2009) suggested that the N input to forests to be less than 5 kg N ha -1 yr -1 . But their model also estimated that, in moist forests, the uptake of N could be 30-80 kg N ha -1 yr -1 . Deluca et al. (2002) reported that cyanobacterium and feather moss could act as a supplement to N2N fixation in boreal forests (0.5 kg N ha -1 yr -1 ) while the organic N accumulation could be 3 kg N ha -1 yr -1 . For the forests in northwest Rocky Mountain, N2N fixation amount is on average between 0.5 and 2 kg N ha -1 yr -1 (Clayton and Kennedy, 1985;Fahey et al., 1988) 330 while Kou-Giesbrecht and Menge's model (2019) estimated the N2N fixation rate to be 0 -10 kg N ha -1 yr -1 for temperate forests, and 0 to 6 kg N ha -1 yr -1 for boreal forests.
There could be a number of reasons for our comparatively lower estimates. The most important one is that there is a considerable uncertainty in estimating the coverage of N2N fixing plants. High diversity in the distribution of legume plants highly influences the estimation of total plant coverage, because our estimation was based on site-335 level experimental data. In order to improve our understanding, more investigation on legume plant distribution and associated data for N2N fixers is needed, especially in the Middle Asia, South America and Africa.
Large variations of BNF rates exist across terrestrial ecosystems spatially (Figure 3). The global BNF spatial pattern is similar to other estimates (Cleveland et al., 1999;Xu-Ri and Prentice, 2017). The highest N2N fixation rate in tropical regions (more than 50% of the global terrestrial N2N fixation) is primarily due to their warm and moist 340 soil conditions. Further, N2N fixed by human activities became increasingly influential in the past century (Galloway et al., 2002), especially in temperate regions due to their large human population. The anthropogenic N deposition contributed more to soil N than BNF did. As a result, soils became N rich, inhibiting BNF in temperate soils. This could explain why the potential N2N fixation rate was high in temperate ecosystems, but only contributed to 20% of the total fixation.

Major controls on biological nitrogen fixation
In our simulations, the N2N fixation was primarily influenced by soil temperature, moisture and soil nitrogen content. The highest N2N fixation rate in tropical ecosystems is consistent with our sensitivity analysis for temperature and soil moisture. The sensitivity analysis indicated that a 1-3°C increase of temperature led to 7% increase in N2N fixation rate. Nitrogen cycle responds differently between different biomes and legume types. But in 350 general, increasing temperature will accelerate processes in the N cycle. Soil moisture correlates with BNF in a similar way with temperature. A slightly increase of precipitation (10%) increased the nitrogenase activity.
However, the response of N2N fixation to soil water stress is not as sensitive as that to the change in temperature.
Xeric shrubland and savanna in dry tropical areas still contribute greatly to the global N2N fixation, while the contribution of boreal forests, with low temperature, is much lower.

355
BNF is highly regulated by soil nitrogen content. N-deficiency conditions usually favor BNF activities, for example, in xeric shrubland and savanna. Enhancing soil N content will decrease the N2N fixation rate, which is also consistent with our sensitivity analysis. It costs less energy for plants to take up N directly from soils rather than biologically fixing it from the atmosphere (Cannell and Thornley, 2000). However, there is an exception for some areas in tropical ecosystems. Many tropical soils are comparatively rich in nitrogen, but N2N-fixing plants are still 360 active to compensate the nitrogen depletion due to the rapid N cycling (Pons et al., 2007). This explains why N fertilization inhibits the BNF in temperate ecosystems, but BNF is still active in N-rich soils in tropical ecosystems.
In areas where the energetic cost succeeds the demand of N, the BNF rate will be comparatively lower. Sullivan et al. (2014) suggested that there were lower rates of BNF in undisturbed mature forests and higher rate in secondary forests, depending on the balance between N-demand and energy consumption.

Model limitation and future work
The incorporation of BNF into TEM allows us to more adequately simulate nitrogen cycle from natural terrestrial ecosystems. However, there are several limitations in this study.
First, the current model ignores the effect of free-living BNF. Although symbiotic BNF is critical for most natural and semi-natural ecosystems, asymbiotic organisms play an important role in extreme environments such as 370 waterlogged soils and deserts. The importance of symbiotic BNF or fixation by leguminous plants may not be as significant as previously thought. Elbert et al. (2012) suggested that cryptogam contributed nearly half BNF in terrestrial ecosystems, which was up to 49 Tg N yr -1 . In some tropical areas, the spatial N input from free-living bacteria even exceeds symbiotic input (Sullivan et al., 2014). In addition, legumes are not the only source of symbiotic BNF. Some fungi species have the ability to actively fix atmospheric nitrogen. But in most existing 375 models, fungi or mycorrhizae symbioses are not considered due to the limited knowledge about their mechanisms of fixing N (Fisher et al., 2010). A more comprehensive model that covers various types of nitrogen fixation is needed.
Second, the BNF process in our model is calibrated with a limited amount of data, imposing a general set of parameters to all plant species and soil conditions within an ecosystem type. More observational data from natural terrestrial ecosystems is desirable to improve our model.  Williams, 1990;Bouniols et al., 1991;Cabelguenne et al., 1999); SOILN (Wu and McGechan, 1999) Wa lower threshold of water content below which N fixation is totally restrict by rhe defict of soil water 0 APSIM, EPIC (Sharpley and Williams, 1990;Bouniols et al., 1991;Cabelguenne et al., 1999); SOILN (Wu and McGechan, 1999) Wb upper threshold of water content above which N fixation is not limited by rhe defict of soil water 0.5 APSIM, EPIC (Sharpley and Williams, 1990;Bouniols et al., 1991;Cabelguenne et al., 1999); SOILN (Wu and McGechan, 1999 Jarrell and Virginia (1990) 40 23.5 *ARA denotes the acetylene reduction assay method in determining biological N2 fixation rates.