Gaseous reactive nitrogen losses of agricultural systems in China influenced by crop trade

China plays an important role in the international trade of agricultural commodities. Considering the dynamic reactive nitrogen (Nr) losses of agricultural systems in China, a hypothesis was proposed that crop conversion in China would be correlated with the extent of crop trade, influencing Nr losses in agricultural systems. The objective of this study was to verify the hypothesis based on a hybrid approach, which incorporated life cycle analysis (LCA), copula–Markov chain Monte Carlo (MCMC) simulation, and copula sampling. The approach was proven to be of benefit in (a) evaluating Nr losses in crop planting based on a LCA framework, (b) identifying dependencies and co-movements of the correlated variables in planting structures and crop trade using copula–MCMC simulations, and (c) recognizing fluctuations in Nr losses of crop planting in the future using copula-based sampling method. The planting structures and international trade of four types of crops (i.e. wheat, soybeans, maize, and rice) in 20 provinces of China indicated significant correlations, thus supporting the initial hypothesis. With the improvement of self-sufficiency in crop production, especially soybeans, Nr losses from the crop production of China in 2025 and 2030 would decrease by 8.43% and 4.26%, compared with those in 2018 (i.e. 1916.74 kt N).


Introduction
Direct reactive nitrogen (Nr) losses in crop production systems traditionally arise from fertilizer application. Normally, 36%-40% of nitrogen (N) in fertilizer is taken up by the plant and a substantial percentage of N is released into the environment through ammonia volatilization, denitrification, leaching, and runoff (Cui et al 2016). Globalization of agriculture not only aids optimization of land, water, and energy resources to meet growing food demands, but also leads to redistribution of Nr losses through food supply chains (Dalin andRodríguez-Iturbe 2016, Lassaletta et al 2016). For example, abundant Nr losses are transferred through global food and trade supply chains, which are larger than the proportion of direct Nr losses (Cui et al 2016). Also, international trade of crop has intensified perturbation of the global biogeochemical N cycle in both importing and exporting regions (Lassaletta et al 2016). Influenced by the import of N-fixing crops (e.g. soybeans), farmland was normally used for the production of high N-demanding crops (Sun et al 2018). Thus, intensified by international trade of crops, complexity of Nr losses in local farmland should be identified (Lassaletta et al 2014).
Previous studies have predominantly focused on virtual water, land resources, and N input of agricultural activities in the context of food trade (Hu et al 2020). Dalin and Rodríguez-Iturbe (2016) analyzed groundwater depletion for irrigation and incorporated it into global virtual water in relation to global food trade. Huang et al (2019) quantified the indirect inputs of N and land resources in China together with the supply chains of food trade. Based on the input-output analysis of food trade, the complexity of the biogeochemical N cycle was aggravated under globalization of food production and consumption (Oita et al 2016). Field experiments demonstrated that multiple production structures influenced Nr losses in crop production through varied N-fertilizer application . For example, 17.7% of N-fertilizer input was lost through NH 3 emission in rice field, being markedly higher than Nr losses in soybean production . Thus, crop import, except for soybeans, would obviously relieve Nr losses in importing countries through the decrease in N-fertilizer application (Huang et al 2019). Conversely, with consideration of variations in Nr losses between N-fixing and other crops, crop import tended to increase Nr losses in importing countries when the original cropland used for soybean production was converted to other crops (Sun et al 2018). Therefore, crop trade is an important factor that influences crop structures and Nr losses in importing countries.
Multiple environmental impacts (e.g., global warming potential, resource depletion, human toxicity, and land occupation) from agricultural systems have been evaluated and compared within the framework of life cycle assessment (LCA) (Prechsl et al 2017, Brunklaus et al 2018, Lienhardt et al 2019. Specific environmental performance in agricultural activities (e.g. carbon, N, or water footprint) was widely analyzed (Chojnacka et al 2019). Eranki et al (2019) developed a LCA framework to account for the carbon footprint in a crop rotation system. Inputoutput analysis has been incorporated with consideration of inputs of nutrients under supply chains of agricultural activities. For example, Hamilton et al (2018) applied an environmentally extended multi-regional input-output analysis to quantify the summed supply chain emissions of N and phosphorus. However, dynamic variations in crop production and crop trade, influencing the final results of environmental impacts, should be considered by means of network and correlation analyses (Schaffer-Smith et al 2018).
In agricultural systems, Nr losses are not static, and are influenced by dynamic supply chains under multiple time scales. Specifically, given variations in planting structures and fertilizer application, it is difficult to predict definitively the features of Nr losses of agricultural systems in the future. Thus, dynamic LCA methods must be developed to indicate the non-static characteristics of Nr losses in crop production (Zhang et al 2014, Shimako et al 2018. Based on numerical and dynamic simulation methods, many previous studies focused on the variations in arising from a multi-tiered structure and interactions along the life-cycle stages of products or services (Adhitya et al 2011, Pigné et al 2020. For example, Shimako et al (2018) explored variations in carbon emissions with consideration of time dependence on energy consumption under future strategies of carbon capture, storage, and utilization. Incorporating LCA indicators in dynamic simulation models, Levasseur et al (2010) identified the noisy and random features of a life cycle inventory as a function of time. Although previous research measured dynamic Nr losses in agricultural activities, co-movement of planting structures and crop trade were rarely considered. Among the most flexible statistical methods to examine dependencies and co-movements in a complex system, copula functions were applied to analyze variations in crop yields under dynamic climatic conditions (Segura et al 2018, Cai et al 2019, Albulescu et al 2020. For example, farm business profit was predicted based on correlations between agricultural production and the El Niño-Southern Oscillation using a statistical copula method (Nguyen-Huy et al 2020). Copula functions were also used to reveal interactions between water resources availability and electricity demand (Yu et al 2020). In addition to simulating the dependence structures of crop yields and drought hazards, copula functions can be used to estimate the probability of crop residues in the future (Leng andHall 2019, Ribeiro et al 2019). Thus, statistical methods were required to support dynamic analysis of the Nr losses in an agricultural system.
Considering the variation of Nr losses in an agricultural system, a hypothesis that crop conversion in China would be correlated with the extent of crop trade was proposed. The objective of this study was to verify correlated features of farming structures and crop import based on LCA, copula-Markov chain Monte Carlo (MCMC) simulation, and a copulabased sampling method. The approach can (a) evaluate the life cycle inventory of Nr losses in agricultural systems based on a LCA framework, (b) identify dependencies and co-movements of the correlated variables in planting structures and crop import using copula-MCMC simulations, and (c) recognize fluctuations of Nr losses in agricultural systems using a copula-based sampling method.

Methodology
Crop production and its Nr losses fluctuate under the influence of crop trade in China. The varied features of Nr losses may arise from variations in the planting structure of N-fixing and high N-demanding crops. Variation in Nr losses could be verified by potential correlation between farming structures and crop trade to apply adaptive management in future decision-making. Thus, dynamic LCA of Nr losses in crop production was conducted to reveal the potential Specifically, a dynamic LCA framework for evaluating Nr losses in crop production was proposed based on the integration of LCA, copula-MCMC simulations, and copula sampling (figure 1). Four steps were adopted in the approach, composed by (a) setting the system boundary and functional unit of Nr losses in crop production, and (b) analyzing the life cycle inventory of the Nr losses, (c) assessing correlated features in planting structures and crop trade based on copula-MCMC simulations, and (d) deriving distributions of Nr losses using a copula-based sampling method to predict the distributions of Nr losses in the future.

LCA of Nr losses in crop planting system 2.1.1. System boundary and functional unit
The following stages in crop production, associated with Nr losses, were included in the system boundary (figure 1): (a) combustion of fossil fuels in agricultural machines, (b) fertilizer production and application, (c) food waste disposal, and (d) biomass reuse from food waste. Functional unit is 1000 ha of crop fields.

Life cycle inventory analysis 2.1.2.1. Nr losses from fertilized soil
In this study, the Nr volatilization (i.e. N 2 O, NH 3 , and NO x ) was considered to be the main cause of gaseous Nr losses. The Nr losses in maize, soybean, wheat, and rice fields were evaluated based on Zhang et al (2018) and Calvo Buendia (2019). Nr losses in fertilized cropland are caused by NH 3 volatilization from fertilized soil into the atmosphere and its conversion into NO and N 2 O. Based on Calvo Buendia (2019) and Bouwman et al (2002), Nr emissions from managed soils for maize, wheat, rice, and soybean production were estimated with equation (1): where N is the annual amount of Nr produced from atmospheric deposition of Nr volatilized from fertilized soils, F SN i is the annual amount of synthetic N fertilizer applied to the i th crop field (kg N per year), and EF N2O,i and EF NH3−NOx,i are the fractions of synthetic fertilizer N that volatilizes as N 2 O and NH 3 -NO x , respectively, for the i th crop field. With to Calvo Buendia (2019), EF N2O,i and EF NH3−NOx,i in equation (1) for wheat, soybean, and maize production were estimated as 0.01 kg N 2 O-N/kg and 0.11 kg NH 3 -NO x , and EF N2O,i in rice field was estimated as 0.004 kg N 2 O-N/kg N. NH 3 emission of rice field, being different from that of arid cropland, was estimated based on experimental data derived from Wang et al (2018) (table 1). Province NH3-N volatilization

Nr losses from agricultural machinery
Nr losses from agricultural machinery were evaluated based on the emission factor of N 2 O by agricultural machinery. The emission factor is the default value for off-road mobile sources (IPCC 2006) [equation (2)]: where F l is the amount of the l th energy consumption by agricultural machinery, EF l is the emission factor for N 2 O in the l th energy consumption, and E is the Nr losses from agricultural machinery. The N 2 O loss [i.e. EF l in equation (2)] from diesel oil in agricultural machinery was 28.6 kg TJ −1 according to IPCC (2006).

Nr losses from biomass reuse
Percentage of crop residues was estimated based on the average data from Gustavsson et al (2011). Biomass reuse from crop residues was evaluated based on Yue et al (2021). The emission factors for N 2 O associated with fertilizer production were 28.6 kg TJ −1 from fossil energy and 41 kg TJ −1 from ethanol according to IPCC (2006). Percentages of cereals, oilseeds, and pulses in food residues and waste through agricultural processes were 22% and 14%, respectively, based on average data from industrialized Asia (Gustavsson et al 2011). It was assumed that 50% of the crop residues would be composted. N 2 O emission in the process of compost production from crop residues is described in equation (3): where EF i is N 2 O emission factor (i.e. 0.3 g N 2 O/kg to wet mass of crop residues and 0.6 g N 2 O/kg to dry mass of crop residues), and M is mass of crop residues (unit: kg).

Nr losses from fertilizer production
Three types of chemical fertilizers (i.e. nitrogen, potash, and phosphate fertilizers) were considered. Related Nr losses in the production stage of chemical fertilizers are mainly caused by energy consumption in the production process. The amount of energy consumption in the production stage of chemical N fertilizers was referred to Yue et al (2018).

Uncertainty analysis based on copula functions 2.2.1. Correlated features analysis based on copula-MCMC simulations
As planting structure may be influenced by crop import, the co-relationship of the two indicators can be verified based on copula functions. According to Sklar's theorem (Sklar 1959), copula functions can be used to describe the joint distributions of correlated random variables. The joint probability density function (PDF) of the two random variables [i.e. f X,Y (x, y)] can be expressed as equation (4) (Gaupp et al 2019): where X i and Y j are two correlated random variables, X i is the area proportion for the i th crop field, Y j is degree of the self-sufficiency in the j th crop which negatively correlates with the j th crop trade, estimated by ratio of local crop production and crop consumption, F X i (x i ) and F Y j y j are their respective marginal cumulative distribution functions, c(•) is the fitted copula function, and f X i (x i ) and f Y j y j represent the PDFs of X i and Y j , respectively. In this study, Gaussian, Student-t, Clayton, Frank, and Gumbel copulas were used for correlation analysis (table A1) (Nelsen 2006). Bayesian analysis was employed to model uncertainties of copula parameters through identifying their posterior distribution (Ali et al 2018, Liu et al 2020, Zhao et al 2022 [equation (5)]: where p θ|Ỹ and p (θ) indicate the posterior and prior distributions of parameters, and p Ỹ = θ p Ỹ |θ dθ is a normalization constant which denotes marginal likelihood acts. p Ỹ |θ ≡ L θ|Ỹ can be evaluated by the following likelihood function (i.e. L θ|Ỹ ) (Sadegh et al 2017): MCMC algorithms can provide samples from high-dimensional complex distributions for estimating posterior distributions of parameters (Naseri and Hummel 2022). Thus, the posterior distributions of parameters in copula functions can be estimated through hybrid-evolution MCMC algorithm, based on Multivariate Copula Analysis Toolbox (Sadegh et al 2017). Details of hybrid-evolution MCMC algorithm can be found in Sadegh et al (2017).
To select the best-fitted copula functions, two indicators [i.e. root mean square error (RMSE) and Nash-Sutcliffe efficiency coefficient (NSE)] were used to examine the performance of the copula functions [equations (7) and (8)] (Sadegh et al 2017): where RMSE ∈ [0, +∞) and NSE ∈ (−∞, 1] indicate the extents of variation in observed and predicted values, respectively,f i represents the joint probability in observed variables of X i and Y, and f i is their joint probability predicted by copula functions.

Copula sampling
In this study, copula sampling was introduced to estimate the probability distribution of the unknown parameter (i.e. crop planting structures) through the probability distribution known parameter (i.e. crop self-sufficiency) and their joint probability distribution. Based on the best-fitted copula functions in the section 2.2.1, joint probability distribution can be obtained. Thus, probability distribution of crop planting structures can be determined by copula sampling. Copula sampling algorithm was described by the following steps (Cai et al 2014): Step 1 Select suitable copula function according to the methods in the section 2.2.1.
Step 3 Determine the cumulative distribution functions (CDFs) [i.e. w x = F −1 (x) and w y = F −1 (y)] of correlated parameters (i.e. x and y) based on historical data.
Step 4 Generate random number matrix W N×2 = [x, y], based on matrix U N×2 and the CDFs generated in Steps 2 and 3.
Step 5 Obtain random number matrix x ′ based on future range of y (e.g. [y 1 , y 2 ]).

Overview of research area and data source
Maize, wheat, rice, and soybeans were chosen as staple crops for Nr losses analysis in relation to crop trade. The research area was chosen to cover the majority of regions used for the production of the staple crops in China. Thus, 20 provinces were selected for

Nr losses from fertilized soil
As indicated in table A4 of supplementary material, per unit area of Nr losses from rice, wheat, soybean, and maize fields were evaluated. In 2018, the maximum amount of N 2 O was volatilized from fertilized soil for wheat production in Inner Mongolia, and the maximum amount of NH 3 was volatilized from fertilized soil under rice production in Sichuan. The minimum amount of N 2 O and NH 3 was emitted from fertilized soil for soybean production in Anhui. Compared with Nr losses from other crops, the per unit area of Nr emitted from soybean field was the least, while the per unit area of Nr emitted from maize, rice, and wheat fields accounted for more than 97.3% of total Nr losses.

Nr losses from agricultural machinery
The per unit area estimated of Nr losses from agricultural machinery are summarized in table A4 of supplementary material. The maximum amount of N 2 O was emitted from agricultural machinery for rice production in Liaoning in 2018.

Nr from biomass reuses
It was assumed that 50% of crop residues would be composted. The amount of N input in crop residues reuse was allocated to fertilizer application and production. As indicated in table A4 of supplementary material, the maximum amount of N reuse was derived from crop residues in per unit area of rice field in Jiangsu in 2018.

Nr losses from fertilizer production
The per unit area of Nr losses for fertilizer production are summarized in table A4 of supplementary material. The maximum amount of N 2 O was emitted from fertilizer production supporting wheat cultivation in Inner Mongolia, whereas the minimum amount was emitted from fertilizer production supporting soybean cultivation in Anhui in 2018.

Life cycle inventory of Nr losses
Based on statistic data for crop structures in the research area, the total Nr losses in crop production varied with the production structure and sown area of the crop in the research area ( The fertilized soil contributed more than 90% of Nr losses in crop planting. Previously, Nr losses from agricultural production were mainly evaluated based on emission factors and field experiments (Schroeck et al 2019). For example, as indicated by the field experiments from Jiang et al (2018), Nr losses in crop planting varied with fertilizer applications (e.g. conventional fertilizer, urea-ammonium mixed nitrogen fertilizer, stabilized urea, and organic fertilizer). Compared with the Nr losses from crop fields (i.e. 36.94-41.20 and 32.41-38.11 kg N ha −1 in wheat and rice fields, respectively) under urea-ammonium mixed nitrogen fertilizer application from Jiang et al (2018), (a) Nr losses from wheat field in 2018 were 17.92%-68.79% lower in most provinces; and (b) Nr losses from rice field in 2018 were 1.38%-77.20% lower in most provinces.

Correlation analysis of planting structures and crop trade
As indicated by the indicators of RMSE and NSE, the best-fitting association parameters between planting structures and crop trade were described by the copula functions, listed in table A6 of supplementary material. As indicated by their association parameters (i.e. θ), the correlation between production structures and crop trade varied with crop and region. Given that a low RMSE and a high NSE indicate a strong correlation, in most provinces the proportions of maize field were correlated with wheat, maize, and soybean trade, as well as the proportions of soybean, wheat, and rice fields were correlated with wheat trade. Thus, these results supported the original hypothesis that planting structures would be correlated with the extent of crop trade in China.  57%-7.19% in Shanxi, Heilongjiang, and Shaanxi provinces; (c) planting structures of soybeans in 2025 and 2030 would increase by 0.99%-4.1% and 1.23%-3.79% in most northern provinces; (d) planting structures of rice field in 2025 and 2030 would decrease by 1.72%-6.67% and 1.37%-5.42% in Jilin, Jiangsu, Anhui, Hubei, and Jiangxi Provinces. Planting structures of rice field in 2025 and 2030 would increase by 1.38%-4.98% and 2.42%-6.81% in Liaoning, Zhejiang, Guangdong, Guangxi, Sichuan, Yunnan, and Fujian Provinces; and (e) the maximum increase in planting structures of maize, wheat, soybean, and rice fields would observe in Shaanxi, Sichuan, Jilin, and Heilongjiang.

Future trends
Influenced by improvement of self-sufficiency and decrease of crop import in 2025 and 2030, the planting structures of crops would be varied (tables A2 and A7 in supplementary material). For example, compared with the crop structures in 2018, (a) soybean production of Hebei and Liaoning Provinces would increase more than two times in 2025 and 2030, in order to fulfill self-sufficiency of soybeans (i.e. 21% to 40%); (b) wheat planting of Heilongjiang Province would increase more than eight times in 2025 and 2030, in order to fulfill self-sufficiency of wheat (i.e. 95% to 105%); and (c) slight variations (less than 0.5 times) of maize and rice production would occur in 2025 and 2030.  12.15%-12.88%, and 4.52%-6.14% in 2025 and 2030. Conversely, Nr losses of rice production in other regions would increase by 0.06%-43.87% in 2025 and 2030.

Total Nr losses in crop planting
Multiple distributions of Nr losses in crop production were collated to generate distributions of total Nr losses using Latin hypercube sampling (Yue et al 2018). The total Nr losses revealed the following trends (figure 5 and table A9 of   Previous studies demonstrated domestic dynamics of soybean production influenced by international soybean trade (Sun et al 2020). Indicated by the correlation between planting structures and crop trade in table A6 of supplementary material, other crop import (e.g. wheat, maize, and rice) would also influence domestic farming structures. Similarly, food import or food self-sufficiency influenced N flows and other environmental impacts (Shi et al 2016). Thus, under the background of the improvement in grain self-sufficiency in China, dynamics in farming patterns and Nr losses of the staple crops should be revealed. The results of this study indicated that the amount of Nr losses related to the production of staple crops would decrease in 2025 and 2030 with the improvement of the grain selfsufficiency in China. The trend of Nr losses in this study coincided with related findings in some previous studies. For example, as indicated by Sun et al (2020), increase of domestic soybean planting would decrease Nr losses and their environmental pollution. Concurrently, as dynamic environmental burdens of food systems were also influenced by the interactions among food demand, production and trade (Zhao et al 2021), demand-side correlation between food consumption and socio-economic development can also be identified through copula-based methods. Differentiated from N budget accounting in agricultural systems or anthropogenic activities [e.g. Gu et al (2015) and Cui et al (2013)], the system boundary of this study did not incorporate Nr losses from other agricultural products. Thus, variations of total Nr losses from food production of China influenced by international trade will be directions of future research.

Conclusions
Although many previous studies have investigated Nr losses in agricultural systems, research on dynamic Nr losses in crop production under the influence of crop trade is lacking. A hypothesis that crops conversion in China would be correlated with the extent of crop import in China was verified in this study based on LCA, copula-MCMC simulation, and a copula-based sampling method. This approach improved on conventional LCA methods for analysis of dynamic Nr losses by (a) identifying dependencies

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request. and 41801203), Guangdong Basic and Applied Basic Research Foundation (2022A1515012023 and 2019A1515110681), National Social Science Grant (20BGL190), and Guangdong Provincial Key Laboratory Project (2019B121203011). The authors much appreciate the editor and the anonymous reviewers for their constructive comments and suggestions which are extremely helpful for improving the paper.