Virtual nitrogen and phosphorus flow associated with interprovincial crop trade and its effect on grey water stress in China

The grey water footprint (GWF) is defined as freshwater requirements for diluting pollutants in receiving water bodies. It is widely used to measure the impact of pollutant loads on water resources. GWF can be transferred from one area to another through trade. Although pollution flow has previously been investigated at the national level, there has been no explicit study on the extent to which crop trade affects GWF across regions and the associated changes in grey water stress (GWS). This study analyzes pollution flow associated with interprovincial crop trade based on nitrogen (N) and phosphorus (P) loss intensity of three major crops, namely, maize, rice and wheat, which is simulated by a grid-based crop model for the period 2008–2012, and evaluates the spatial patterns of GWS across China. The results indicate that the integrated national GWF for N and P was 1271 billion m3 yr−1, with maize, rice, and wheat contributing 39%, 37%, and 24%, respectively. Through interprovincial crop trade, southern China outsourced substantial N and P losses to the north, leading to a 30% GWS increase in northern China and 66% GWS mitigation in southern China. Specifically, Jilin, Henan, and Heilongjiang Provinces in the northern China showed increases in GWS by 161%, 114%, and 55%, respectively, while Fujian, Shanghai, and Zhejiang in the south had GWS reductions of 83%, 85%, and 80%, respectively. It was found that the interprovincial crop trade led to reduced national GWF and GWS. Insights into GWF and GWS can form the basis for policy developments on N and P pollution mitigation across regions in China.


Introduction
Nitrogen (N) and phosphorus (P) fertilizer applications are crucial for crop production. They boost production to facilitate feeding of the growing world population (West et al 2014). Globally, N and P fertilizer application in agricultural fields in 1960-2010 increased nine-and three-fold, respectively (Sutton et al 2013). However, the global N and P use efficiencies decreased from 0.58 to 0.44 and from 0.82 to 0.61, respectively, over the period 1950-2000 (Bouwman et al 2013). This has been associated with continuous increases in N and P losses (Liu et al 2016a(Liu et al , 2018, which reached 120 Tg yr -1 for N and 2-7 Tg yr -1 for P in agricultural soils over the period 2000-2010 (Sutton et al 2013). Substantial N and P losses have resulted in notable environmental issues such as water quality degradation and groundwater contamination, which has caused eutrophication worldwide (Galloway andCowling 2002, Vitousek et al 2009).
According to Hoekstra (2003), the water footprint (WF) is a multidimensional indicator of water use. It includes green water (soil moisture in the unsaturated zone), blue water (surface and groundwater), and grey water (polluted water). The grey water footprint (GWF) measures the volume of water required to dilute pollutants in water bodies to an acceptable degree based on natural background concentrations and ambient water quality standards Chapagain 2008, Hoekstra andMekonnen 2012). Accordingly, grey water stress (GWS) is defined as the ratio of GWF to total water resources and quantifies the degree to which water available to dilute pollutants has been consumed by GWF (Liu et al 2017). GWF and GWS link water quantity to quality, and therefore, they have been widely used to evaluate the impact of agricultural production on water resources at various spatial levels Hoekstra 2015, 2017).
Crop trade involves the virtual flow of resources, such as water and land (Hoekstra and Hung 2005) as well as the redistribution of pollution, such as N and P loads in water (O'Bannon et al 2014, Lun et al 2018. Since the inception of the virtual water concept (Allan 1996), substantial work has been conducted on the quantification of virtual water transfer at various spatial and temporal scales (Dalin et al 2012, Hoekstra and Mekonnen 2012, Suweis et al 2013, Hoekstra 2017. These studies have helped improve our understanding of the role of crop trade in compensating for local water deficiency (Yang and Zehnder 2007). In comparison to intensive studies on virtual water trade, studies on virtual pollution flow and the related effects on local GWS are limited (Liao et al 2019). O'Bannon et al (2014) constructed a global network of virtual grey water flow for 1986−2010 and found that international food trade has led to pollution that is more heavily concentrated in a few regions, such as the United States and Canada. Zhao et al (2019) quantified global GWF for 40 countries/regions with 35 economic sectors using a multiregion input-output approach and noted that agriculture accounted for the largest portion of virtual grey water flow. A major limitation of these studies is the lack of information on the extent to which GWF flow changes trade partners in the GWS.
The GWF is particularly high in China (Liu et al 2017), primarily because of the large amount of food production and intensive use of both N and P fertilizers in the country (Wang et al 2018). GWS is further worsened by the mismatch of major agricultural production regions in northern China and the large water resource availability in southern China. Currently, northern China exports large amounts of food to southern China. How N and P flows are associated with crop trade, pollution, and GWS in China remains unclear.
To fill this knowledge gap, the objectives of the present study are (a) to estimate GWF and GWS in each province/municipality (hereafter referred to as province) related to the production of three major crops (maize, rice, and wheat) in China, (b) to calculate virtual N and P flows embodied in interprovincial crop trade, and (c) to investigate the effects of interprovincial crop trade on the spatial distribution of GWS in China. Through virtual N and P flows, food-importing provinces outsource N and P pollution, thereby mitigating local GWS. Conversely, food exporting provinces retain extra N and P loads, which leads to additional GWS in the local environment. This study provides useful information for agricultural pollution management in China and beyond.

Methods and data description
We first estimated N and P losses into water bodies at the grid level from the production of maize, rice, and wheat for the period 2008-2012. These three crops account for 79% of the sown area of grain crops and 54% of the total agricultural sown area in China (NBSC 2012). Then gridded N and P losses were aggregated to the provincial level with an areaweighted method. Combining interprovincial trade data (averages of the period 2008-2012) of these three crops, we investigated the virtual N and P flows and the associated effects on GWS in China. Only 31 provinces in mainland China were included. Taiwan, Hong Kong, and Macao were not included in this study because of a lack of data.

Estimation of nitrogen and phosphorus losses and the related uncertainties
In this study, gridded N and P losses into water bodies from the three crops were simulated using the Python-based Environmental Policy Integrated Climate (PEPIC) model (Liu et al 2016a(Liu et al , 2018) at a spatial resolution of 0.5 arc degrees. PEPIC is based on the field crop model Environmental Policy Integrated Climate (EPIC) (Williams et al 1984), which simulates crop growth, nutrient dynamics, soil erosion, and hydrological processes on a daily basis. Potential biomass was estimated by multiplying the intercepted solar radiation by the crop-specific biomass energy coefficient. To derive the actual biomass, potential biomass was then scaled by the major plant stress, including water, nutrients, temperature, aeration, and salinity. At harvest, crop yield could be estimated based on accumulated actual biomass during the growing season and the crop-specific harvest index.
PEPIC applied N and P inputs as chemical fertilizers and manure. N inputs were applied twice (before planting and one and a half months after planting), and P inputs were applied once just before planting as base fertilization. PEPIC simulated the losses of N and P in water via three pathways, namely, surface runoff, percolation, and eroded soil particles. N and P losses through runoff and percolation were simulated as a product of the respective water fluxes and N and P concentrations in runoff and percolation, while losses of N and P with eroded soil were simulated by considering soil erosion, N and P concentrations in the topsoil layer, and respective N and P enrichment ratios in the soil. These processes considered the impacts of soil types, fertilizer inputs, topography, precipitation patterns, and soil sorption capacity on P (Liu et al 2016a(Liu et al , 2018. The large-scale PEPIC model has been shown to perform well in representing crop yields, as well as N and P losses at the country or regional levels (Müller et al 2017, Franke et al 2020, Jägermeyr et al 2020, Liu et al 2020. For example, the simulated ratios of N (P) losses to N (P) inputs by PEPIC were within the ranges of previous studies, and simulated P losses in China matched well with regional reports (Liu et al 2016a(Liu et al , 2018. In this study, we used the verified large-scale PEPIC model without further model modification.
The simulations of N and P by PEPIC are also subject to considerable uncertainties, particularly those related to model parameters (Wang et al 2012, Della Peruta et al 2014. To demonstrate the impacts of the parameters on the model simulation, 28 parameters related to N and P dynamics were selected for uncertainty analyses based on Wang et al (2012) and the EPIC manual (table S1 available online at stacks.iop.org/ERL/16/124018/mmedia). The Latin hypercube sampling (LHS) method (Mckay et al 1979) was used to resample the parameters within their specific ranges. Using the LHS method, PEPIC was run 100 times for each crop, with 100 parameter combinations. The 100 outputs of simulated N and P losses were then used to estimate the GWF associated with N and P. According to the 100 simulations of N and P losses associated with 100 parameter combinations, variables related to GWF and GWS also had 100 outputs. The mean values of the 100 outputs were used to demonstrate the simulated results, while the coefficient of variation (CV, the ratio of mean to standard variation) was used to present the uncertainties. This uncertainty analysis method was successfully applied in previous studies by the present authors (Liu et al 2016a(Liu et al , 2020.

GWF and stress calculation
The total losses of N and P from the three crop fields (maize, rice and wheat) in each province were estimated using the following equations: where TNL (i,δ) and TPL (i,δ) are total losses of N and P (kg yr −1 ) due to crop δ production in province i; NL (i,δ) and PL (i,δ) are simulated N and P losses per hectare (kg ha −1 yr −1 ) in province i for crop δ using the PEPIC model (table S2); and SArea (i,δ) is the area (ha) under crop δ in province i. The total GWF associated with each crop in each province was estimated as: where GWFN (i,δ) and GWFP (i,δ) are the total GWF related to N and P losses from crop δ in province i (m 3 yr −1 ), NC max and PC max are water quality standards for total N and total P (mg l −1 ), and NC nat and PC nat are the natural background concentrations of total N and total P in the receiving water bodies (mg l −1 ), respectively. In this study, NC max and NC nat indicate the total N concentration of the quality standard of surface water for drinking (10 mg l −1 ) and natural water (0.4 mg l −1 ). PC max and PC nat imply the total P concentration in the quality standard of surface water at Class III (0.2 mg l −1 ) and natural water (0.01 mg l −1 ) in China. The values of NC nat (0.4 mg l −1 ) and PC nat (0.01 mg l −1 ) in the receiving water bodies were derived from previous studies with large scale applications (Franke et al 2013, Mekonnen and Hoekstra 2015. However, these values are subject to changes, for example, several studies set the natural concentration at 0 (Pellegrini et al 2016, Wu et al 2016. A previous study by the present authors demonstrated that setting PC nat (0.4 mg l −1 ) and PC nat (0.01 mg l −1 ) to 0 could underestimate the GWF (Liu et al 2017). The number '1000' was used for the unit conversion. In addition to the total GWF, the GWF per unit area (m 3 ha −1 yr −1 ) was determined by scaling up the total GWF to provincial administrative areas. For a province, the maximum GWFN and GWFP were used to represent the integrated GWF (GWFI) because the volume of water required to dilute one pollutant simultaneously dilutes other pollutants. Therefore, the integrated GWF was estimated as: where GWFI (i,δ) is the integrated GWF of N and P associated with crop δ in province i; max is a mathematical expression that finds the maximum value. GWFP and GWFI, GWSN, GWSP and GWSI (integrated grey water stress of N and P) were calculated in a similar way to GWFN. In this study, the focus was on the GWSI, estimated as: where GWSI (i,δ) (dimensionless) refers to GWS due to GWFI (i,δ) and TWR i represents the total water resources (including surface water and groundwater, m 3 yr −1 ) in province i.

Virtual N and P flow estimation via crop trade
Virtual N and P flows embodied in interprovincial crop trade were estimated as: where VNF i,j and VPF i,j indicate virtual N and P (ton yr −1 ) embodied in (maize, rice, and wheat) trade flow from province i to j, and NI (i,δ) and PI (i,δ) indicate N and P loss intensity per unit production of crop δ in province i (kg ton −1 ) (table S2), defined as the ratio of total N and P losses to crop production.

Interprovincial crop trade impact on GWS
To detect the possible impact of interprovincial crop trade on GWS, additional/mitigated GWS due to virtual N and P flows in a province were calculated as: (9) where TE GWS is the effect of interprovincial crop trade on GWS (%), GWS cur is GWS under current conditions, and GWS nt is GWS under hypothetical conditions without interprovincial crop trade. Therefore, food-importing provinces need to produce extra food originally bought from food exporting provinces, whereas food exporting provinces reduce food production for export.

Data description
The inputs for the PEPIC model included geographical data (latitude, longitude, elevation, and slope), climate data, soil properties, N and P fertilizer application, land use data, and phenology data. N and P inputs, including chemical fertilizer and manure data, were based on Zuo et al (2018) and Tong et al (2017). The N and P chemical fertilizer inputs for the three crops in 2010 were derived from Zuo et al (2018), and the total N and P inputs from livestock manure in 2010 were derived from Tong et al (2017). The manure data were mapped to all of the agricultural land in each province. Table S3 details the chemical fertilizer and manure application rates for N and P at the provincial level in China. Soil properties were derived from the ISRIC-WISE dataset (Batjes 2006), which provided the information necessary for PEPIC, including organic carbon, percentages of sand and silt, pH, bulk density, and layer depth. Seven climatic variables, namely, maximum air temperature, minimum air temperature, relative humidity, wind speed, precipitation, solar radiation, and CO 2 concentration, were obtained from the WATCH Forcing Data methodology applied to ERA-Interim data (WFDEI) dataset (Weedon et al 2014). Crop-specific land use data were obtained from the SPAM2010 dataset at a spatial resolution of 0.5 arc degrees (IFPRI (International Food Policy Research Institute) 2019). SPAM2010 is currently the most updated grid-based and crop-specific cropland distribution dataset and was used to aggregate the provincial N and P losses. Data on interprovincial crop trade (Trade (i,j,δ) ) were estimated by Zhuo et al (2019) by analyzing the differences in food production and food consumption of each province. Provinces with surplus food were considered to be exporting regions. Otherwise, they were considered to be importing regions. The trade matrix was then simulated using a linear optimization approach (Dalin et al 2014) to minimize the trade cost. Statistical data on provincial crop production, yield, cultivated area, total water resources, water use, and administrative area were obtained from the National Bureau of Statistics of China (http://data.stats.gov.cn).
All abbreviations and their descriptions used in this study are provided in table 1.

Crop GWF
The GWF per unit area due to maize production is higher in northern China than in southern China (figures 1(a) and (d)). For example, the GWFN and GWFP per unit area in northern China are 465 and 518 m 3 ha −1 yr −1 , respectively, and those for southern China are 213 and 326 m 3 ha −1 yr −1 , respectively. The largest GWFN per unit area is in Shandong Province (2698 m 3 ha −1 yr −1 ), and the largest GWFP is in Liaoning Province (3206 m 3 ha −1 yr −1 ). The total volume of GWF due to maize production in China is 354.2, 427.5 and 493.0 billion m 3 yr −1 for GWFN, GWFP, and GWFI, respectively (table S4). Generally, GWF is mostly concentrated in northern China, accounting for 76% of the total GWFI related In contrast to maize, the GWF per unit area due to rice production is higher in southern China than in northern China (figures 1(b) and (e)). The GWFN and GWFP per unit area in northern China are 103 and 42 m 3 ha −1 yr −1 , respectively, and those for southern China are 965 and 875 m 3 ha −1 yr −1 , respectively. The total volumes of GWF related to rice production in China for GWFN, GWFP, and GWFI are 417.5, 348.9 and 465.3 billion m 3 yr −1 , respectively (table S4). GWF due to rice production occurs mostly in southern China, with GWFN, GWFP, and GWFI of 356.7, 323.7 and 404.5 billion m 3 yr −1 , respectively, accounting for 85%, 93%, and 87% of the total volume of rice production, respectively, in China. The GWFI is relatively high in Jiangxi, Hunan, Guangdong, Guangxi, and Jiangsu Provinces, together contributing to 56% of the total GWFI due to rice production.
For GWF per unit area associated with wheat production, GWFN and GWFP in northern China are 286 and 315 m 3 ha −1 yr −1 , respectively, and those for southern China are 97 and 282 m 3 ha −1 yr −1 , respectively (figures 1(c) and (f)). The largest GWFN and GWFP per unit area existed in Shandong and Henan Provinces, which were 2537 and 3329 m 3 ha −1 yr −1 , respectively. For the total volume of GWF related to wheat production in China, GWFN,GWFP,and GWFI are 205.4,290.8,and 313.2 billion m 3 yr −1 , respectively (table S4). In northern China, GWFI associated with wheat production occurs mainly in Shandong, Henan, Hebei, and Anhui Provinces, together accounting for 48% of the total in China. In southern China, it is mainly observed in Jiangsu and Hubei Provinces, together accounting for 18% of the total.
The GWF per unit area associated with the three crops (maize, rice, and wheat) is higher in southern China than in northern China (figure 2). The GWFN, GWFP, and GWFI values per unit area are 854, 876, and 1084 m 3 ha −1 yr −1 for northern China and 1275, 1483, and 1702 m 3 ha −1 yr −1 for southern China, respectively. However, the highest GWFN (5315 m 3 ha −1 yr −1 ) and GWFP (4637 m 3 ha −1 yr −1 ) per unit area occur in Shandong and Henan Provinces, respectively. The total volumes of GWFN, GWFP, and GWFI for the three crops across the country are 977.1, 1067.1, and 1271.5 billion m 3 yr −1 , respectively (table S4).
A comparison of GWFN and GWFP for all three crops indicates that GWFI is generally dominated by GWFP. In a few provinces with severe water shortages, such as Beijing, Tianjin, Inner Mongolia, Shanxi, and the provinces in the northwest region, GWFN has a greater effect than GWFP. This is because P is more bondable to soil and, therefore, less vulnerable to loss when precipitation is lower. In northern China, GWF is mainly related to the production of maize and wheat because of the prodigious amount of cultivation. In southern China, however, GWF is mostly related to rice production. Maize, rice, and wheat contribute 39%, 37%, and 24% to the total GWFI in China, respectively.

Spatial distribution of integrated GWS
When estimated as a ratio of GWF to total water resources in a region, the GWSI shows different spatial patterns for the three crops across China (figure 3). In accordance with the distribution of GWFI because of maize production, the provinces with GWSI larger than one are primarily found in  3(a)). Although the GWF is mostly concentrated in southern China, it does not have much influence on the local water pollution of the GWSI of rice because of the abundant water resources in the region ( figure 3(b)). For the GWSI due to wheat production, a high GWSI generally exists in Ningxia, Henan, Hebei, and Shandong Provinces ( figure 3(c)). For northern China, the GWSIs due to maize, rice, and wheat production are 0.62, 0.10, and 0.35, respectively. The values for southern China are 0.06, 0.19, and 0.05, respectively.
For the GWSI of the three crops together, a high GWSI is generally observed in northern China, including Ningxia, Henan, Hebei, Shanxi, and Shandong Provinces ( figure 3(d)). Ten provinces have a GWS larger than one, indicating that the dilution capacity of local water resources is fully consumed, and water quality is deteriorating in those provinces. For the northern and southern China, the GWSIs associated with the three crops are 1.06 and 0.29, respectively.

Virtual N and P flow embodied in crop trade
Virtual N flow occurs mainly by exportation from northern China to southern China of maize and  The virtual P flow embodied in the interprovincial crop trade has a similar pattern to the virtual N flow, that is, from north to south for maize and wheat, and from south to north for rice. The largest virtual P flows embodied in maize, rice, and wheat trade are from Jilin to Hunan, Hunan to Henan, and Henan to Guangdong, accounting for 1.6, 2.1 and 2.7000 tons yr −1 , respectively. The total virtual P flow amounts to 63.8000 tons yr −1 , among which 39.3000 tons yr −1 are embodied in crop trade from north to south, and 10.5000 tons yr −1 in the reverse direction.

Provincial N and P outsourcing and retention
Estimated as the product of local nutrient loss intensity per unit production (NI and PI) and inflow and outflow of crop trade, the provincial outsourcing and retention of N and P are presented in table 2. In northern China, the total outsourced N and P are 297 and 1.7000 tons yr −1 , respectively, and the total retained N and P are 1024 and 36.9000 tons yr −1 , respectively. Consequently, northern China has additional pollutant retention associated with domestic crop trade. Henan Province has the largest retention of both N and P, reaching 262 and 8.6000 tons yr −1 , respectively. Conversely, N and P outsourcing in southern China is 3268 and 195.8000 tons yr −1 , and the retention is 78 and 2.3000 tons yr −1 , respectively. Therefore, pollution is mitigated in southern China through domestic crop trade. Hunan has the largest outsourcing of N, and Jiangxi has the highest outsourcing of P. Total outsourcing of N and P is significantly higher than retention for all of China. This implies that interprovincial crop trade is efficient in reducing national N and P pollution in China.

Interprovincial crop trade effect on integrated GWS
The effect of crop trade on the GWSI varies with crop type (figure 5). The GWSI of maize increased by 57% in northern China relative to the GWSI without crop trade, while mitigating the GWSI by 79% in southern China ( figure 5(a)). Jilin and Heilongjiang have the largest increases in the GWSI of 502% and 381%, respectively, but in many other provinces in southern China, for example, Zhejiang, Shanghai, Fujian, Jiangxi, and Hainan, the GWSI significantly decreased. In contrast to maize, the GWSI driven by rice production dropped by 46% because of crop trade in northern China but increased by 23% in southern China ( figure 5(b)). However, in Heilongjiang and Jilin, the extra production of rice for  export increased the GWSI by 83% and 20%, respectively. For the GWSI driven by wheat production, interprovincial trade led to an additional increase in the GWSI of 47% in northern China but an 89% decrease in southern China (figure 5(c)). In Henan Province, wheat exports caused a 232% increase in the GWSI. Although Xinjiang Province has a high impact because of extra wheat production (increased by 142%), the GWSI is only 0.1 because of wheat production in the province. The interprovincial crop trade of the combined three crops led to an additional increase in the GWSI in northern China ( figure 5(d)). In Jilin Province, extra food production led to an additional increase in the GWSI of 161%. However, in some metropolises such as Beijing and Tianjin and some less developed provinces such as Gansu and Qinghai, GWS is mitigated to various extents. In contrast, crop trade mitigates the GWSI in most provinces in southern China, including Fujian, Shanghai, and Zhejiang. In productive provinces such as Jiangsu, however, the GWS increases. The interprovincial crop trade of the combined three crops led to an additional increase of 30% in the GWSI in northern China, while it mitigated the GWSI in southern China by 66%, suggesting that the GWS shifted from southern China to northern China via crop trade.

Findings from uncertainty analyses
Uncertainties in GWF and GWS as well as trade flow derived from PPEIC parameters (100 parameter combinations) are presented in figures S1 and S2 and table 2 as coefficients of variation (CV) and standard deviation. The CV values of GWFN and GWFP for each crop are mainly (88% of total provinces) within 40%, and most (97% of the total provinces) of them are lower than 60%, except for Qinghai Province (figure S1). The CVs of GWFN are lower than those of GWFP, and the CVs of GWF of rice are lower than those of maize and wheat. The CVs for both GWFN and GWFP in northern China are higher than those in southern China. The CVs of GWSI for the three crops together range from 9% to 40%, with an average of 23% (figure S2). In nearly 60% of the provinces, the CVs of the GWSI of the three crops together are between 20% and 30%. Large uncertainties can also

Interprovincial crop trade effect on water pollution
Currently, most provinces in northern China produce surplus staple food, particularly maize and wheat, to support consumption in southern China (figures S3-S5), thereby adding an extra burden on the local environment. In a global study, O'Bannon et al (2014) noted that pollution was becoming more concentrated in a few exporting countries. Without stringent regulations on water and fertilizer use, crop trade can result in a further deterioration of local water resources and environmental conditions in exporting regions. Importing regions can consider compensating for exacerbated water pollution issues in exporting regions. Otherwise, their own GWSI would exceed sustainable levels, that is, higher than 1 in some provinces.
Interprovincial food trade in China is generally criticized because of inefficient blue water use. Dry and irrigation-intensive northern provinces produce extra food for wetter and less irrigated southern provinces (Dalin et al 2014. However, the findings of this study suggest that interprovincial trade does not always have negative effects on environmental sustainability. It reduces N and P losses across the entire country. This is primarily because N and P are easily lost in wet southern provinces (Liu et al 2016a(Liu et al , 2018, where crop yields are generally low. This has resulted in higher N and P loss intensity per unit production in southern China, particularly for maize and wheat. Therefore, exporting maize and wheat from the north to the south has benefited N and P pollution in China. Under the current conditions, GWS should be considered in a comprehensive investigation of the impacts of interprovincial crop trade on the environment. For example, although crop trade led to an additional increase in the GWSI in Heilongjiang (114%) and Inner Mongolia (56%), the GWSI in these provinces remains less than one. This suggests that the capacity of local water resources to dilute pollutants in the region is sufficient. For Beijing and Tianjin, although the GWSI is mitigated to various degrees by crop trade, the GWSI value exceeds one. This implies that the pollution level in these regions is far beyond the dilution capacity of local water resources in the region, and it is necessary to increase crop imports to mitigate environmental degradation.

Challenges in safeguarding water resources in northern China
As an indicator of water resource appropriation in pollution, GWF facilitates an assessment of the sustainable, efficient, and equitable use of water resources (Franke et al 2013). According to the present study, 1271 billion m 3 yr −1 of freshwater was needed to dilute both N and P pollutants due to the cultivation of the three staple crops from 2008 to 2012 in China. The GWF is mainly distributed in the northeast, north-central, and Yangtze regions. For GWS, the most severe areas were mainly located in the north-central region (except for Ningxia Province). This is the most extensive agricultural production zone in China, yet it has the most severe water scarcity (figure S6) because of diminishing streamflow and groundwater depletion (Xiao et al 2017, Rodell et al 2018. In the Hutuo River Basin, which is a typical headwater basin in the northern China Plain and the main source of groundwater for the central region of the plain, runoff has decreased to only onethird of the value in the 1950s in nonflood seasons (Cheng et al 2014, Han et al 2019. The coincidence of large amounts of N and P loads with limited water resources worsens the GWS situation in the northcentral region of the plain.
Water endowment in northern China accounts for less than 20% of the national total water resources and arable cropland in the region accounts for 65% of China's total, with grain production reaching 57% of the national total in 2016. This mismatch between water resources and food production not only worsens water scarcity, but also causes water pollution. Northern China has significantly less water to dilute the pollution caused by virtual water transfer through food trade. The combined effect of water scarcity and pollution poses a significant challenge to water security in northern China.

Water pollution mitigation in northern China
The application of N and P fertilizers is critical for increasing crop production. However, low fertilizer use efficiency and fertilizer overdose can result in substantial nutrient loss from croplands into water bodies. This is the case in China, where the GWF is particularly high, contributing 30% to the global GWF driven by the P load (Mekonnen and Hoekstra 2017). In addition, the P input per unit area in China is significantly higher than that in many other countries in the world, and the P use efficiency is, therefore, far below the global average (Lun et al 2018).
The northeast and north-central regions are the main food production bases in northern China. They export a large amount of food and, therefore, cause additional environmental stress. The overuse of fertilizers in this region is a common strategy to improve crop production. For example, N fertilizer application per unit of cropland in the north-central region (approximately 369 kg N ha −1 for maize and wheat) is approximately 1.7 times greater than that in the High Plains (approximately 217 kg N ha −1 for the two crops) in the United States (Pei et al 2015). Therefore, reducing fertilizer inputs and improving the efficiency of fertilizer use have the potential to mitigate GWS in northern China. This agrees with the current policy. In 2015 the Zero Increase Action Plan for fertilizer use by the 2020 government plan was announced by the Ministry of Agriculture to restrict the overuse of fertilizers (Liu et al 2016b).
In contrast, adjusting to crops with less N balance is a promising way to reduce GWF driven by crop production in northern China. For example, per hectare, the N balance in cropland under soybean is negative because of biological N fixation, while per hectare N balance in cropland under wheat, corn, rice, and vegetables is positive (Sun et al 2018). The total N load will be significantly reduced by replacing other crops with a higher N balance with soybeans in major N loss regions.

Conclusions
This study estimated GWF and GWS driven by N and P fertilizer applications in the production of three major crops (maize, rice, and wheat) over the period 2008-2012 in China. Virtual N and P flows embodied in interprovincial crop trade and their effects on GWS were evaluated. The results show that the GWS in northern China is substantially higher than that in southern China, mainly because of the high crop cultivation intensity and limited water resources in northern China. The provinces in southern China outsource substantial N and P loads to northern China via crop trade, which adds to the already worsening environmental burden in northern China and in turn challenges local water resource management in the region. Crop trade has largely influenced the spatial distribution of GWS in China and shifted it from south to north. The uncertainty analyses indicated that model parameters could have a substantial impact on the estimation, and hence should be carefully considered. The findings of this study suggest that the water pollution shift via crop trade should be considered when formulating robust water management policies aimed at alleviating the existing imbalance between the water-scarce and water-rich regions of China.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).