# Global Phosphorus Losses from Croplands under Future Precipitation Scenarios

Wenfeng Liu, Philippe Ciais, Xingcai Liu, Hong Yang, Arjen Hoekstra, Qiuhong Tang, Xuhui Wang, Xiaodong Li, Lei Cheng

► To cite this version:

Wenfeng Liu, Philippe Ciais, Xingcai Liu, Hong Yang, Arjen Hoekstra, et al. Global Phosphorus Losses from Croplands under Future Precipitation Scenarios. Environmental Science and Technology, 2020, 54 (22), pp.14761-14771. 10.1021/acs.est.0c03978. hal-03025964

HAL Id: hal-03025964

https://hal.science/hal-03025964v1

Submitted on 6 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
This document is the accepted manuscript version of the following article:
Liu, W., Ciais, P., Liu, X., Yang, H., Hoekstra, A. Y., Tang, Q., ... Cheng, L.
(2020). Global phosphorus losses from croplands under future precipitation
scenarios. Environmental Science and Technology, 54(22), 14761-14771.

https://doi.org/10.1021/acs.est.0c03978

This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/

Global phosphorus losses from croplands under future precipitation scenarios

Wenfeng Liu¹,²,*, Philippe Ciais¹, Xingcai Liu³,⁴, Hong Yang³,⁵, Arjen Y. Hoekstra⁶,⁷, Qiuhong Tang⁴,⁸, Xuhui Wang¹, Xiaodong Li⁹, Lei Cheng¹⁰,¹¹

¹Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France.

²College of Water Resources and Civil Engineering, China Agricultural University, Beijing 100083, China.

³Eawag, Swiss Federal Institute of Aquatic Science and Technology, Ueberlandstrasse 133, CH-8600 Duebendorf, Switzerland.

⁴Key Laboratory of Water Cycle and Related Land Surface Processes, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China.

⁵Department of Environmental Sciences, MGU, University of Basel, Petersplatz 1, CH-4003 Basel, Switzerland.

⁶Twente Water Centre, University of Twente, Enschede, The Netherlands.

⁷Institute of Water Policy, Lee Kuan Yew School of Public Policy, National University of Singapore, Singapore.

⁸College of Resources and Environment, University of Chinese Academy of Sciences, Beijing, China.

⁹State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, No. 24 South Section 1, Yihuan Road, 610065 Chengdu, China.

¹⁰State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China.

¹¹Hubei Provincial Collaborative Innovation Center for Water Resources Security, Wuhan 430072, China.

*Corresponding author.

Email: wenfeng.liu@lsce.ipsl.fr or wfliu2012@gmail.com

Telephone: +33 1 69 08 87 39
# 1 Abstract

Phosphorus (P) losses from fertilized croplands to inland water bodies cause serious environmental problems. During wet years, high precipitation disproportionately contributes to P losses. For the first time, we combine simulations of a gridded crop model and outputs from a number of hydrological and climate models to assess global impacts of changes in precipitation regimes on P losses during the 21st century. Under the baseline climate during 1991–2010, median P losses are 2.7 ± 0.5 kg P ha⁻¹ yr⁻¹ over global croplands of four major crops, while during wet years P losses are 3.6 ± 0.7 kg P ha⁻¹ yr⁻¹. By the end of this century, P losses in wet years would reach 4.2 ± 1.0 (RCP2.6) and 4.7 ± 1.3 (RCP8.5) kg P ha⁻¹ yr⁻¹ due to increases in high annual precipitation alone. The increases of P losses are highest (up to twofold) in arid regions of Middle East, Central Asia and northern Africa. In three quarters of the world's river basins, representing about 40% of total global runoff and home up to seven billion people, P dilution capacity of freshwater could be exceeded due to P losses from croplands by the end of this century.

15

# 16 TOC/Abstract Art

<figure><img src="image_2.png" /></figure>

17

18

# 19 1. Introduction

Losses of phosphorus (P) from croplands to aquifers, streams, lakes and finally to oceans have resulted in severe environmental consequences around the world¹, such as emerging

2

3

4

5

6

7

8

9

10

11

12

13

14

16

18

20

21
eutrophication² and spreading anoxic dead zones in the coastal ocean³. P losses are projected to increase from croplands⁴, ⁵, with more fertilizers needed to sustain a higher food production⁶. But P losses are also sensitive to climatic factors⁷, ⁸. To keep P loads below critical regional and global boundaries⁹, it is necessary to quantify and understand the interaction between precipitation regimes and P losses from fertilized croplands.

Under global climate change, increasing frequency and magnitude of high precipitation are projected¹⁰⁻¹², which could be an important driver of P losses in intensive cultivation areas in the future. Michalak¹³ outlined the overlooked impact of (increasing) climate extremes on water quality, especially for intensified P and nitrogen (N) losses. Ockenden et al.¹⁴ found that precipitation had strong effect on P losses and future climate (mainly precipitation) alone would increase winter P loads up to 30% by 2050 in three catchments in the U.K. Using P load measurements in two tributaries of Lake Mendota in the U.S., Carpenter et al.¹⁵ found that precipitation extremes were associated with extremes of discharge and substantial increases in P loads. In the same region, Motew et al.¹⁶ applied a process-based watershed model to detect the impacts of P inputs and precipitation on P loads and found a significant and positive relationship between these two factors and P concentration in water bodies. On the one hand, an increase in the balance between precipitation and evapotranspiration causes higher runoff, and this process will dilute P pollution. On the other hand, higher precipitation will increase the intensity of P losses on croplands e.g. shown in Carpenter et al.¹⁵. It is thus the balance between these two mechanisms that determine how P pollution will change with precipitation in the future. So far, studies on influence of precipitation changes on P losses related to agriculture have only been conducted in a few specific watersheds¹⁴⁻¹⁶. Examining future P losses based on the projections of the General Circulation Models (GCMs) under different climate scenarios is still missing on a global scale. Further, previous studies did not provide quantitative estimates
of exposed populations to degraded water quality across different catchments throughout the
world.

49

The objectives of this study are to project regional P losses at global scale from croplands in
future wet years and to determine whether increased runoff in wet years is sufficient to dilute
increased P loads. We estimated global P losses through fitting a substitute linear model trained
to reproduce the simulation outputs of the process-based PEPIC crop model¹⁷, ¹⁸. Then, the
substitute model was used to detect the impacts of future precipitation in wet years on annual P
losses. To investigate the water pollution induced by P losses, the concept of grey water
footprint (GWF) was used¹⁹. We quantified the amount of freshwater required to dilute P
pollution below a critical concentration threshold of 0.02 mg P L⁻¹ based on the GWF
guideline²⁰. Then the grey water stress (GWS) was calculated to indicate water P pollution
levels²¹, ²². We focused on the detection of regions where GWS will reach above 1, meaning
that local water resources are insufficient to dilute the P loads²¹, and assessed the number of
people exposed to severe P pollution. This study, which shows the validity of a linear substitute
model to estimate the impacts of future precipitation in wet year on P losses, helps to quantify
the regional mechanisms that determine how P loads will change with precipitation under
climate change.

65

## 2. Material and methods

### 2.1. Simulation of phosphorus losses

In this study, P losses were defined as losses of P from cropland soils with surface runoff and leaching water as well as with soil erosion. The PEPIC model is used to simulate global P losses from these sources on croplands at a resolution of 0.5°¹⁸. PEPIC is a grid-based version of the Environmental Policy Integrated Climate (EPIC) model²³ coded under the Python environment. It performed well in simulating crop yields, crop water use, and N dynamics at large scales¹⁸.

72
24, 25. More recently, PEPIC was successfully used to assess global P losses by conducting an integrative crop-soil-management modelling approach¹⁷. The ratio of P losses simulated by PEPIC to P inputs is in the middle range of previous reports and regional P losses estimated by PEPIC matched well with reported ones (Tables 2 and S8 in Liu et al.¹⁷). In addition, PEPIC presented good and comparable ability in representing inter-annual variation of crop yields between 1980 and 2010 at the country level among 14 large scale crop models²⁶, ²⁷ participating in the Agricultural Model Intercomparison and Improvement Project (AgMIP)²⁸. In this study, we used the verified PEPIC model to estimate P losses. Detailed information on the PEPIC model and estimation of P losses can be found in the Supporting Information.

The PEPIC model was run 100 times for each crop using 100 sets of parameter combinations to investigate uncertainties associated with model parameters. We carefully chose 19 model parameters (related to crop growth, water balance, as well as N, P and carbon routines) and their possible ranges based on Wang et al.²⁹ and the EPIC user's manual³⁰ (Table S1). Parameter sampling was conducted by applying the Latin Hypercube Sampling method³¹, ³². These parameters and sampling method were used in our previous study¹⁷.

Four major crops, i.e. maize, rice, soybean and wheat, were included in this study. P fertilizer consumption for these four crops together accounts for 73% of total P fertilizers among 17 most commonly produced crops³³. P inputs from chemical fertilizer and manure were applied before crop planting along with tillage, which is a widely used P application method³⁴. P losses for each grid were estimated as average P losses weighted by cropland areas of the four crops. P losses at the river basin level were aggregated by using area-weighted average over cropland in the river basin. It deserves to note that the basin-scale P losses defined here are the direct losses without considering any biochemical processes of P in water bodies. The estimated P loads could be higher than delivered P from the outlet of a given basin. Here, the river basin map was
derived from the FAO GeoNetwork (http://www.fao.org/geonetwork/srv/en/metadata.show?id=38047), consisting of 228 watersheds in total. P losses in the basins with the smallest cropland areas (for a total of 0.1% of global total cropland areas of the four crops) were not considered in the analysis.

103

## 2.2. Projection of future trend in P losses

PEPIC simulates P losses given precipitation and other climate variables through a set of complex equations controlling the water balance of the soil and its P budget. As runoff is roughly proportional to P losses, and runoff is related to precipitation (and evapotranspiration), one can expect a positive correlation between P losses and precipitation. We developed a substitute linear model from the model outputs by regressing modelled P losses from each of 100 simulations during the baseline period 1991-2010 as a function of precipitation at annual level for each river basin, according to:

$$ \ln(Pl) = slope \times \ln(pr) + intercept \quad (1) $$

where Pl is time series of total annual P losses [kg P ha⁻¹ yr⁻¹] averaged from all the cropland grid cells in each river basin (188 basins considered here) over the croplands; ln is natural log transformation, which is used to consider possible nonlinear relationship between P losses and precipitation³⁵; pr is the time series of annual precipitation [mm yr⁻¹], which was estimated as an area-weighted average of precipitation according to the cropland fraction of each grid cell; slope [ln(kg P ha⁻¹ mm⁻¹)] and intercept [ln(kg P ha⁻¹ yr⁻¹)] are model parameters. The slope presents the response of ln(Pl) to one unit change of ln(pr) over the river basin considered. It also indicates the elasticity of P losses to precipitation, i.e. 1% of change in precipitation corresponding to slope% of change in P losses. The coefficient of correlation (r) and p-value by the student-t test were used to detect the significance of the relationship.

123
In order to determine the impacts of high annual precipitation on P losses, we calculated the median values of P losses during the wet years (P losses of high precipitation, referred to as $Pl_{h,2000}$) and the median P losses (referred to as $Pl_{m,2000}$) during all the years in the period 1991-2010 (referred to as subscript 2000). Wet years are defined as the years with annual precipitation above 75th percentile of the whole distribution in a given period. Then the differences between $Pl_{h,2000}$ and $Pl_{m,2000}$ were treated as the additional P losses due to high annual precipitation. We then regressed $\ln(Pl_{h,2000})-\ln(Pl_{m,2000})$ against slope×($\ln(pr_{h,2000})-\ln(pr_{m,2000})$), where $pr_{h,2000}$ and $pr_{m,2000}$ are the high and median precipitation in the baseline period, to see whether it is appropriate to use the substitute linear model to detect the impacts of high precipitation on P losses as:

$$ \ln(Pl_{h,2000}) - \ln(Pl_{m,2000}) = a_h \times \text{slope} \times [\ln(pr_{h,2000}) - \ln(pr_{m,2000})] \quad (2) $$

where $a_h$ is fitted parameter. In addition to the comparison for high P losses, we also investigated the agreement for maximum P losses and minimum P losses, as:

$$ \ln(Pl_{max,2000}) - \ln(Pl_{m,2000}) = a_{max} \times \text{slope} \times [\ln(pr_{max,2000}) - \ln(pr_{m,2000})] \quad (3) $$

$$ \ln(Pl_{m,2000}) - \ln(Pl_{min,2000}) = a_{min} \times \text{slope} \times [\ln(pr_{m,2000}) - \ln(pr_{min,2000})] \quad (4) $$

where $Pl_{max,2000}$ and $Pl_{min,2000}$ are P losses in the year with maximum ($pr_{max,2000}$) and minimum ($pr_{min,2000}$) baseline precipitation, respectively; $a_{max}$ and $a_{min}$ are model fitted parameters. The robustness of relationships expressed in eqs. (2-4) was detected by the coefficient of determination ($R^2$) and p-value. The results suggested quite good performance (Fig. 1d and Fig. S1). Besides, we found that the slopes derived from eq. 1 are mainly influenced by coefficient of variation (cv) of precipitation (Fig. 1c), which will not change much in the future. The future precipitation in wet years generally ranges between baseline minimum and maximum precipitation (Fig. S2). Therefore, we used the slope obtained from the historical regression model, assuming that the slope would not change under future climate conditions, to project the impacts of future high precipitation on P losses for four periods: 2020-2039 (referred to as
subscript 2030), 2040–2059 (referred to as subscript 2050), 2060–2079 (referred to as subscript 2070), and 2080–2099 (referred to as subscript 2090).

$$Pl_{h,future} = Pl_{m,2000} * \left( \frac{pr_{h,future}}{pr_{m,2000}} \right)^{slope} \quad (5)$$

where $Pl_{h,future}$ is future P losses of high precipitation [kg P ha⁻¹ yr⁻¹]; slope is derived from eq. 1; $pr_{h,future}$ is the high precipitation in the four future periods [mm yr⁻¹]. Eq. 5 was only applied to river basins where the numbers of 100 simulations with significant (p-value < 0.05) linear relationship between ln($Pl$) and ln($pr$) are more than 50, while for the other basins (responsible for only about 2% of global croplands of the four crops), $Pl_{h,future}$ was kept the same as $Pl_{h,2000}$.

To examine whether the linear substitute model performs well in representing the PEPIC simulated P losses under future climate conditions, we ran the PEPIC model forced by future climate data in eight major river basins around the world. For simplicity, we only ran PEPIC under two Representative Concentration Pathway scenarios (RCP2.6 and RCP8.5) for the periods 2030 and 2090. Then we compared future P losses estimated by PEPIC to those estimated by the linear substitute model.

## 2.3. Assessment of water pollution intensity associated with phosphorus losses

The concepts of GWF and GWS were used to assess the water pollution intensity due to P losses into water bodies. GWF was introduced by Hoekstra et al.¹⁹ and it measures the water requirements to dilute pollutants (P losses in this study) based on ambient water quality standard and natural background concentration. It is estimated as:

$$GWF = 100 \cdot Pl/(C_{max} - C_{nat}) \quad (6)$$

where GWF is the grey water footprint over cropland [mm yr⁻¹] associated with P losses; $C_{max}$ and $C_{nat}$ [mg P L⁻¹] are ambient water quality standard and natural background concentration
of P. According to the GWF concept, the range between C<sub>nat</sub> and C<sub>max</sub> defines the dilution capacity of freshwater. The values of C<sub>max</sub> and C<sub>nat</sub> were selected as 0.02 mg P L<sup>-1</sup> and 0.01 mg P L<sup>-1</sup> according to the GWF accounting guideline<sup>20</sup>. These values were also used by Mekonnen and Hoekstra<sup>22</sup>. The GWS was calculated as a ratio [dimensionless] of GWF to runoff simulated by large-scale hydrological models (HMs), indicating to which level the dilution capacity of freshwater has been consumed by pollutants<sup>21, 22</sup>. Here, GWF is over cropland, while runoff is over terrestrial land. Hence, GWF is scaled by a ratio of cropland to terrestrial land for each river basin. With this concept, GWS > 1 for a given type of pollutant means that freshwater has exhausted its dilution capacity of this pollutant. Therefore, we focused on the regions with GWS greater than 1 to explore to what extent land and runoff are facing severe P pollution conditions and how many people are residing in the high P pollution regions for both the baseline situation (2000) and future projection (2030, 2050, 2070, and 2090). Similar to P losses, we calculated median GWF and GWS, and GWF and GWS of high precipitation.

187

## 2.4. Data description

Input data of PEPIC include elevation, soil, climate, fertilizer, crop calendar and crop land use<sup>18</sup>. Climate data for PEPIC historical simulation were obtained from the WFDEI dataset<sup>36</sup>, including precipitation, temperature, wind speed, and relative humidity at a daily step. Nutrient inputs of P and N from mineral fertilizer and manure were derived from the EarthStat dataset<sup>33</sup>,<sup>37</sup> Crop land use data were based on the MIRCA2000 dataset<sup>38</sup>, which provides irrigated and rainfed land area for 26 crop species. Detailed information of input data of PEPIC can be found in the Supporting Information.

196

Future annual precipitation data (2020-2099) were downloaded from the data archive of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP, www.isimip.org)<sup>39</sup>. The ISIMIP climate data include five CMIP5 GCMs (GFDL-ESM2M, HadGEM2-ES, IPSL-CM5A-LR,
200 MIROC-ESM-CHEM, and NorESM1-M) forced by four RCP scenarios (RCP2.6, RCP4.5,
201 RCP6.0 and RCP8.5). These precipitation data were downscaled to a resolution of 0.5° and
202 bias-corrected with a trend-preserving method<sup>40</sup>. The five GCM models account for 55% of
203 uncertainties of the entire CMIP5 GCM models<sup>41</sup>. Annual runoff data were provided by four
204 HMs, i.e. DBH<sup>42</sup>, H08<sup>43</sup>, PCR-GLOBWB<sup>44</sup>, WaterGAP2<sup>45</sup>. These runoff data were generated
205 for the ISIMIP fast track phase<sup>39</sup>. These models performed well in representing global large-
206 scale hydrological cycles<sup>46</sup>. Grid-based population (s ?) data for the baseline and future periods
207 were also obtained from the ISIMIP archive. Future populations were estimated based on the
208 Shared Socioeconomic Pathway scenario 2.

209

210 In the analysis, we assumed that fertilizer P inputs and land use patterns remain unchanged for
211 the future at the baseline level. As it is likely that P inputs will increase in the future in response
212 to the higher food demand<sup>47</sup>, future water P pollution might be underestimated in this study.

213

## 3. Results

215 **3.1. Responses of phosphorus losses to precipitation**

216 The regression between PEPIC-simulated P losses and annual precipitation time series (both
217 are in natural log form) shows a strong positive linear relationship across most river basins of
218 the world (Figs 1a and S3). Among 188 river basins considered in this study, 177 of them (94%)
219 have more than 50 simulations out of 100 presenting statistically significant positive
220 relationships (*p*-value < 0.05) between natural log transformed P losses from croplands and
221 natural log transformed annual precipitation. The median coefficient of correlation (r) for
222 simulations with significant regression relations is higher than 0.6 for 87% of the basins. The
223 slope of the linear regression model in each river basin reflects the overall elasticity of P losses
224 to precipitation (see Materials and Methods). This slope shows high spatial variations and takes
225 the highest values in the North China, the central U.S., and the Middle East catchments (Fig.
226

1b). The values of the slope are mainly influenced by the cv of precipitation, as the partial correlation between the slope and cv of precipitation across different catchments reaches 0.7 when controlling for P inputs and median precipitation (Fig. 1c). These results indicate that a linear substitute model relating P losses to precipitation with catchment specific slope derived from the fully-fledged PEPIC simulations can explain most of the spatial variance of present-day estimated P losses.

227

228

229

230

231

232

233

The linear substitute model also has a good ability to reproduce year to year changes in P losses from precipitation, as shown by the comparison of abnormal log P losses during high annual precipitation (log of mean P losses during the 75th highest precipitation years in each catchment minus log of median P losses during all the years in the period 1991-2010) between the results from the original PEPIC model and the substitute linear model (Fig. 1d). In addition, the linear substitute model has a good performance for reproducing maximum minus median P losses (log form), and median minus minimum P losses (log form) from the fully-fledged PEPIC (Fig. S1).

234

235

236

237

238

239

240

It also performs well in representing PEPIC estimated P losses forced by future climate data in eight major river basins throughout the world. The estimated future P losses with the linear substitute model vs. with the PEPIC model were found to be comparable ($R^2$ = 0.86 between the two approaches across the eight basins; Fig. S4). We noticed large standard deviation for the East Brazil basin in the comparison of P losses estimated by both methods. This is mainly because future high precipitation in East Brazil has much larger standard deviation across GCMs than in other basins. In addition, the level of high annual precipitation in the future generally lies in the range between baseline minimum and maximum precipitation and cv of precipitation would not change much into future (Fig. S2). This suggests that the substitute linear model calibrated on the range of baseline precipitation variability is not extrapolated outside this range when forced by future precipitation fields. Besides, the slope values for the eight river basins derived from the future simulations also match well with the baseline slope

241

242

243

244

245

246

247

248

249

250

251
values (Fig. S5), further confirming robustness of the linear substitute model for estimating future P losses. However, the approach makes an assumption that the interannual baseline elasticity of P losses to precipitation is conserved for each catchment under future climate conditions. It ignores, for instance, effects of changes in precipitation frequency and intensity within a year, changes in agricultural practices and in cultivated areas from each crop type.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Relationship between natural log transformed phosphorus losses (ln(Pl)) and natural log transformed annual precipitation (ln(pr)) at the river basin level. **a:** coefficient of correlation (r) of the linear regression model between ln(Pl) and ln(pr) for the period 1991–2010. **b:** slope [ln(kg P ha⁻¹ mm⁻¹)] of the linear regression model. **c:** influential factors of the slope. **d:** comparison of differences in ln(Pl) (losses during the wet years minus the median of all years) between the fully fledged PEPIC model (x-axis) and the linear substitute model (y-axis). In the plot, basins with the smallest areas (for a total of 0.1% of global total cropland</figcaption></figure>
areas) and numbers of simulations with significant relations (*p*-value < 0.05) between ln(PI) and ln(pr) lower than 50 are discarded. In subplot **c**, p_r[slope, pr-cv|(pr, Pin)] is the spatial partial correlation coefficient across basins between slope and coefficient of variation of pr (pr-cv) when controlling pr and P inputs (Pin), p_r[slope, pr|(pr-cv, Pin)] is (?) the spatial partial correlation coefficient between slope and pr when controlling pr-cv and Pin, and p_r[slope, Pin|(pr-cv, pr)] is the spatial partial correlation coefficient between slope and Pin when controlling pr-cv and pr. In subplot **d**, the equation represents the linear relationship of the dashed red line; R² is the coefficient of determination of the equation; the dashed blue line is the 1:1 line.

275

276 **3.2. Effects of baseline phosphorus losses on freshwater pollution**

The baseline median annual P losses show quite similar spatial patterns as median annual precipitation at the river basin level (Figs 2a and S6), with high P loss values in wet countries and regions like Japan, Korea, south-eastern China, Southeast Asia, and South America, where P inputs are also generally high. In these regions, the GWF values are extremely high, reaching over 4000 mm yr⁻¹ in some hotspot regions, meaning that a large amount of freshwater is needed to dilute P losses (more than available when GWS > 1). About 71% (133/188) of the river basins considered have already today exceeded their dilution capacity of P pollution, especially in China, India, the U.S., and Argentina (Fig. 2b).

285
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Impacts of high precipitation on phosphorus losses (PI) and grey water stress (GWS) for the period 1991–2010. a: median PI [kg P ha⁻¹ yr⁻¹]. b: median GWS [dimensionless]. c: percentage differences [%] between PI during wet years and median PI across all years. d: percentage differences [%] between GWS during wet years and median GWS across all years.</figcaption></figure>

Globally, precipitation in wet years (above the 75th percentile of all the years) is about 14% higher than median precipitation, but P losses during these wet years are 34% larger than median P losses under baseline climate (3.6 ± 0.7 (mean ± s.d.) kg P ha⁻¹ yr⁻¹ of P losses during wet years compared to 2.7 ± 0.5 kg P ha⁻¹ yr⁻¹ during median years) (Fig. 3a). This disproportionate increase in global P losses relative to precipitation during wet years is mainly due to the slope values in the linear regression model are predominantly larger than 1. Spatially, the relative increases in P losses during high vs. median precipitation years are the largest in dry regions, e.g., the Middle East, Central Asia, the western U.S., northern Africa and South Africa, and southern Australia (Fig. 2c). High P losses during wet years are even two folds higher than median P losses in some hotspot regions like the Middle East, Central Asia and northern Africa. In these regions, the relative differences between high and median precipitation
values are larger than elsewhere (Fig. S6). Consequently, GWF during wet years increases strongly in these regions. While P losses increase in most basins during wet years, this does not imply that the resultant GWS in these basins always increase, because increased runoff also favours the dilution capacity of P pollution in freshwater. Although increases in GWS shown in most regions in wet years, decreased values can also be observed (Fig. 2d). We classified 140 river basins with GWS > 1 during wet years, compared to 133 during normal (median) years. This increase (5%) in the number of basins exceeding P-thresholds in wet years suggests that large P losses occur in those basins, leading to very high pollution.

About 38% ± 4% (under normal years) and 42% ± 4% (under wet years) of global land area already experiences P pollution with baseline GWS > 1 (Fig. S7a). The related catchments
314

represent about 40% of total global runoff (Fig. S7b). About 5 billion people are affected by P polluted waters (Fig. 3b).

315

316

<figure><img src="image_5.png" /></figure>

317

318

319

320

321

322

<caption>**Fig. 3.** Responses of phosphorus losses (Pl) to high precipitation (pr) and affected population with grey water stress exceeding 1. a: median Pl and high Pl during wet years (kg P ha⁻¹ yr⁻¹) for the baseline period (around the year 2000) and high Pl during wet years for four future periods (2030, 2050, 2070 and 2090) and four Representative Concentration Pathway (RCP) scenarios (RCP2.6, RCP4.5, RCP6.0 and RCP8.5). b: total affected populations (billion).</caption>
Whiskers plots give 95th, 5th, median and inter-quartile values. Red square points represent mean values. M stands for median, H for high. Background colors indicate different periods.

### 3.3. Impacts of future wet years on freshwater pollution

With the projected increases in high precipitation given by GCMs (Fig. S8), future P losses are simulated to increase gradually from 4.0 ± 0.9 kg P ha⁻¹ yr⁻¹ in 2030 to 4.5 ± 1.1 kg P ha⁻¹ yr⁻¹ in 2090 (Fig. 3a). The increases under scenario RCP2.6 are more pronounced in 2050, but less in 2070 and 2090. P losses under RCP2.6 are higher than those under the other three scenarios during the periods of 2030 and 2050, mainly due to the higher precipitation over the croplands under RCP2.6 (Fig. S8). The highest global average P losses are found for 2090 under scenario RCP8.5, with losses in wet years of 4.7 ± 1.3 kg P ha⁻¹ yr⁻¹, which is 75 ± 48% higher than baseline median P losses and 31 ± 36% higher than baseline high P losses. Percentage increases in P losses are larger than 200% in the Middle East, north-eastern Africa, and south Asia (Fig. 4a), but relative decreases are found in a few basins, e.g. northern Africa and Northwest China, where future high precipitation is lower than baseline median precipitation due to a general drying trend in these areas (Fig. S2). In many basins, GWS under future high precipitation is estimated to be smaller than under baseline high precipitation (Figs 2d and 4b), holding the fertilizer inputs constant. In some regions, the GWS during wet years in 2090 will become even lower than GWS under the baseline median precipitation, e.g. in several basins in northern Africa, southeastern Brazil, and Southeast Asia. These decreases in GWS are mainly due to the increases in runoff under future high precipitation to be larger in magnitude than the estimated increases in P losses, a process which dilutes P pollution. Nevertheless, GWS during the four future periods exceeds 1 in most intensively cultivated catchments (Fig. S9), especially in South and East Asia, West Europe and eastern U.S., accounting today for about 90% of global cropland of the four crops.
349

Although the land areas and runoff volumes affected with GWS > 1 are projected to increase only slightly during wet years in the future (Fig. S7), the population exposed to severe P pollution will increase from 5 billion people under baseline high precipitation to 7 billion people under future high precipitation mainly due to increases in future population (Fig. 3b). We also observed that increases in affected population level off after 2050 and exposed populations become lower in 2090 than 2050. This is mainly because the increases in total populations are projected to slow down after 2030 and the numbers in some populous regions start to decline after 2050, e.g. in the Ganges basin in India and the Yangtze basin in China.

350

351

352

353

354

355

356

357

<figure><img src="image_6.png" /><figcaption>**Fig. 4. Relative changes in phosphorus losses (Pl) and grey water stress (GWS) in 2090 compared to baseline situation (around the year 2000).** a: percentage differences [%] of Pl of high precipitation in 2090 in comparison to median P losses in 2000. b: percentage differences [%] of GWS of high precipitation in 2090 in comparison to median GWS in 2000.</figcaption></figure>

358

359

360

361

362

363
364 **3.4. Uncertainties analyses**

We find that uncertainties in projecting future high P losses are more related to slope values,
derived from different PEPIC parameter combinations, than different precipitation outputs from
GCMs (Fig. 5). As for affected populations, both slope values and combination of GCMs and
HMs play important role on the projection uncertainties. In addition to uncertainties associated
with slope values, GCMs and HMs, we also explicitly investigated the uncertainties due to the
selection of different thresholds of percentile for defining wet years, that is, 65th, 70th, 75th, 80th,
85th, 90th, and 95th of the distributions of both baseline and future climates (Fig. S10). A large
difference in P losses during wet years can be observed when choosing varied thresholds. P
losses will become higher when using a higher percentile for wet year definition, with the
highest P losses reaching about 6 kg P ha⁻¹ yr⁻¹. However, there is only a minor effect on the
affected populations with GWS > 1 under the different wet year definitions (?).

376

<figure><img src="image_7.png" /></figure>

377

<caption>**Fig. 5.** Uncertainties related to slopes, climate models (GCMs) and hydrological models (HMs). **a:** phosphorous losses (Pl) [kg P ha⁻¹ yr⁻¹] under future wet years averaged across different GCMs. **b:** Pl under future wet years averaged across different slope values. **c:** affected</caption>

378

379

380
populations [billion] with grey water stress > 1 averaged across different GCMs and HMs. d:
affected population with grey water stress > 1 averaged across different slope values.

383

384

4. Discussion

The impacts of high precipitation on P losses during the baseline period are the largest, even up to a twofold increase, over dry regions with high sensitivity of P losses to precipitation, like in the Middle East, Central Asia, and the western U.S. (Figs 2b and S6a). This is mainly due to high variations of precipitation, causing accumulation of P in the soil during dry years followed by high P losses during wet years in these arid regions (Fig. S2). Therefore, more attention should be paid to fertilizer management in dry lands for lessening the impacts of high precipitation on P-related pollution, e.g. determining P application rates and timing based on crop P demand and prior soil P information and avoiding P fertilization just after heavy rainfall event<sup>48</sup>. Moreover, little freshwater is available to assimilate pollution in those dry areas, which leads to further degraded water quality. On the other hand, P losses are already high (Fig. 2a) in humid regions concurrent with high P inputs<sup>17</sup>. Although future high precipitation would not increase P losses that much in relative term in these humid regions, any increase in P losses in absolute term will further deteriorate water quality there. Interestingly, we find that future high precipitation will increase P losses, but the affected land areas and runoff with GWS > 1 will be almost consistent with baseline high precipitation conditions (Fig. S7). Although future GWS in wet years would increase in many highly polluted river basins, it still stays lower than 1 in most less polluted basins, as of the baseline conditions. This implies that the dilution effect of future increased runoff is comparable to the effect of a future increase in the flow of P from cropland soils to water in these less polluted basins with baseline GWS < 1. This finding holds consistent with different percentiles for defining wet years (Fig. S10).

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405
406

It is worth noting that the GWF concept is based on the dilution effect of freshwater, which is different from the capacity of inland water P assimilation. Although assimilation effect was also referred in previous GWF studies⁴⁹, it implies to account of removal mechanisms of P in water bodies, e.g. microbial uptake, sediment burial and sorption on mineral particles before delivery to rivers. Modeling detailed future changes of P concentration in the cascade of inland-water involves complex chemical and biological processes. In this study, we did not track inland-water P concentration changes under future climate change, as it would require hydrological and biogeochemical models and detailed local data to simulate the transport and reactivity of P, which is out of the scope of this study. Addressing this question should be integrated into the future development of GWF theory.

407

408

409

410

411

412

413

414

415

416

417

Unlike nitrate pollution, which could cause methemoglobinemia and cancers in infants⁵⁰, high P concentration in water does not directly affect human health. However, with a growing number of people living in regions with degraded aquatic environments (e.g. dilution capacity for P will be exhausted in about three fourths of global total river basins shown in this study), continuous P pollution will pose significant indirect impacts on human well-being⁸. High level of precipitation will further amplify P losses under projected future climate changes, particularly in high P input regions in China, India, and eastern U.S.⁴ Sinha et al.⁵¹ reported that high precipitation will increase N losses in these regions. The concurrent additional losses in P and N triggered by high precipitation could greatly challenge ecological and humans' health, as both elements have largely transgressed their global safe planetary boundary⁵², ⁵³. However, P could cause more serious environmental problems due to its much lower acceptable concentrations in water bodies. Therefore, it is necessary to take urgent actions to reduce P losses, e.g. through precision farming and better nutrient management¹⁷, ³⁷, international food trade for improving global nutrient use efficiency⁵⁴, ⁵⁵, and better P recycling and recovery⁵⁶.

418

419

420

421

422

423

424

425

426

427

428

429

430

431
The modelling framework developed in this study provides an effective and robust tool to quantify the impacts of changing precipitation on P losses, especially in arid regions with high slope values in the P losses and precipitation regression analysis (Fig. 1b). However, the results are still subject to uncertainties. For instance, uncertainties related to fertilizer types and as well as application timing affects the estimation results⁵⁷. Due to data unavailability, we had to assume simple P application method on large scale with only one global P dataset. We considered baseline land use pattern unchanged for the future conditions, which could lead to conservative estimation of P losses because global cropland area is still growing⁵⁸. We apply the baseline relations between P losses and precipitation to estimate future P losses, which relies on future precipitation regimes only, but do not directly simulate future P losses based on crop modelling. Lacking future crop-specific P fertilizer information and difficulties on handling nutrient management in crop modelling for a long period⁵⁹ prevent us from doing a direct centurial simulation with PEPIC for the globe. The approach to translate future changes in high precipitation to P losses is verified for baseline precipitation regimes and several major river basins by running PEPIC with future climate. The future contrasts between high and median precipitation are found not far out of the range from the baseline contrast values (Fig. S2a). P fertilizer consumption is expected to continuously increase towards the end of this century⁶⁰. Although this increase in P uses may not increase the elasticities of P losses to annual precipitation (Fig. 1c), it will certainly increase P losses¹⁷. We also noticed that CV of future precipitation in few river basins is slightly different from the baseline CV (Fig. S2b). This change, to some degree, could affect the projection of P losses under the future climate. Therefore, detailed local investigation is needed in these specific regions in future studies.

In addition to the agricultural sector, P losses from other sectors, e.g. domestic and industry, also generate P loads in inland water bodies. For instance, exposed populations with GWS > 1 increase by 0.5-1 billion (Fig. S11) when considering industrial and domestic P losses obtained
from Mekonnen and Hoekstra ²². As non-agricultural P losses are with high uncertainty to project⁶⁰, we did not consider them in this study. Instead, we isolated the impacts of precipitation changes on P losses, as changes in future precipitation can hardly be avoided through local management. We focused on inter-annual patterns of precipitation in this study, while changing distribution on precipitation within a year will also affect P losses¹⁵. Therefore, investigation on the effects of P input increases and intra-annual precipitation changes on P losses should be the focus of future research. This study, for the first time to our best knowledge, explicitly estimates the impacts of future precipitation changes on P losses and highlights the hotspot regions and associated negative consequences on a global scale. It is of particular importance to inform policy decisions for controlling P-related pollutions in the context of changing climate.

469

470 Author contributions

W.L., P.C. and H.Y. designed the research; W.L. and X.L. collected the ISIMIP data; W.L. run the PEPIC model, analyzed the data and wrote the draft; all authors participated in the interpretation of results and the writing and editing process.

474

475 Competing financial interests

The authors declare no competing financial interest.

477

478 Acknowledgments

This work is dedicated to the memory of A.Y. Hoekstra. We thank Prof. Dr. M.M. Mekonnen from the University of Nebraska for sharing datasets of domestic and industrial phosphorus losses. We are grateful to the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) for providing the climate, runoff, and population data. This study was supported by the funding from the Swiss National Science Foundation (P2EZP2_175096), the French State Aid managed
484 by the ANR under the "Investissements d'avenir" programme with the reference ANR-16-CONV-0003, the European Research Council Synergy project SyG-2013-610028
485 IMBALANCE-P, and the Open Fund Research project by the State Key Laboratory of Hydraulics and Mountain River Engineering (Sichuan University) (SKHL1709). X.L. and Q.T.
486 are supported by the Strategic Priority Research Program of Chinese Academy of Sciences
487 (XDA20060402) and National Natural Science Foundation of China (41790424).
488
489

490

491 References

492 1. Chen, M. P.; Graedel, T. E., A half-century of global phosphorus flows, stocks, production, consumption, recycling, and environmental impacts. *Global Environ Chang* **2016**, *36*, 139-152.

493 2. Carpenter, S. R., Phosphorus control is critical to mitigating eutrophication. *P Natl Acad Sci USA* **108**, 105 (32), 11039-11040.

494 3. Diaz, R. J.; Rosenberg, R., Spreading dead zones and consequences for marine ecosystems. *Science* **2008**, *321*, (5891), 926-929.

495 4. Bouwman, L.; Goldewijk, K. K.; Van Der Hoek, K. W.; Beusen, A. H.; Van Vuuren, D. P.; Willems, J.; Rufino, M. C.; Stehfest, E., Exploring global changes in nitrogen and phosphorus cycles in agriculture induced by livestock production over the 1900-2050 period. *Proc Natl Acad Sci USA* **110**, (52), 20882-7.

500 5. Sutton, M., A.; Bleeker, A.; Howard, C. M.; Bekunda, M.; Grizzetti, B.; de Vries, W.; van Grinsven, H. J., M.; Abrol, Y. P.; Adhya, T. K.; Billen, G.; Davidson, E. A.; Datta, A.; Diaz, R.; Erisman, J. W.; Liu, X., J.; Oenema, O.; Palm, C.; Raghuram, N.; Reis, S.; Scholz, R. W.; Sims, T.; Westhoek, H.; Zhang, F. S., *Our Nutrient World: the challenge to produce more food and energy with less pollution*. Centre for Ecology and Hydrology (CEH): 2013.

505 6. Tilman, D.; Balzer, C.; Hill, J.; Befort, B. L., Global food demand and the sustainable intensification of agriculture. *P Natl Acad Sci USA* **2011**, *108*, (50), 20260-20264.

509 7. Tong, Y.; Zhang, W.; Wang, X.; Couture, R.-M.; Larssen, T.; Zhao, Y.; Li, J.; Liang, H.; Liu, X.; Bu, X.; He, W.; Zhang, Q.; Lin, Y., Decline in Chinese lake phosphorus concentration accompanied by shift in sources since 2006. *Nature Geoscience* **2017**, *10*, (7), 507-511.

512 8. Michalak, A. M.; Anderson, E. J.; Beletsky, D.; Boland, S.; Bosch, N. S.; Bridgeman, T. B.; Chaffin, J. D.; Cho, K.; Confesor, R.; Daloglu, I.; DePinto, J. V.; Evans, M. A.; Fahnenstiel, G. L.; He, L. L.; Ho, J. C.; Jenkins, L.; Johengen, T. H.; Kuo, K. C.; LaPorte, E.; Liu, X. J.; McWilliams, M. R.; Moore, M. R.; Posselt, D. J.; Richards, R. P.; Scavia, D.; Steiner, A. L.; Verhamme, E.; Wright, D. M.; Zagorski, M. A., Record-setting algal bloom in Lake Erie caused by agricultural and meteorological trends consistent with expected future conditions. *P Natl Acad Sci USA* **110**, (16), 6448-6452.

518 9. Elser, J.; Bennett, E., Phosphorus cycle: A broken biogeochemical cycle. *Nature* **2011**, *478*, (7367), 29-31.

520 10. Fischer, E. M.; Knutti, R., Observed heavy precipitation increase confirms theory and early models. *Nature Climate Change* **2016**, *6*, (11), 986-991.

522 11. Stott, P., How climate change affects extreme weather events. *Science* **2016**, *352*, (6293), 1517-8.

524 12. IPCC *Managing the risks of extreme events and disasters to advance climate change adaptation*; Cambridge University Press, Cambridge, UK, and New York, NY, USA, 582 pp., 2012.

526 13. Michalak, A. M., Study role of climate change in extreme threats to water quality. *Nature* **2016**, *535*, (7612), 349-350.

528 14. Ockenden, M. C.; Hollaway, M. J.; Beven, K. J.; Collins, A. L.; Evans, R.; Falloon, P. D.; Forber, K. J.; Hiscock, K. M.; Kahana, R.; Macleod, C. J. A.; Tych, W.; Villamizar, M. L.; Wearing,
530 C.; Withers, P. J. A.; Zhou, J. G.; Barker, P. A.; Burke, S.; Freer, J. E.; Johnes, P. J.; Snell, M. A.;
531 Surridge, B. W. J.; Haygarth, P. M., Major agricultural changes required to mitigate phosphorus losses
532 under climate change. *Nat Commun* **2017**, *8*, (1), 161.

533 15. Carpenter, S. R.; Booth, E. G.; Kucharik, C. J., Extreme precipitation and phosphorus loads
534 from two agricultural watersheds. *Limnology and Oceanography* **2018**, *63*, (3), 1221-1233.

535 16. Motew, M.; Booth, E. G.; Carpenter, S. R.; Chen, X.; Kucharik, C. J., The synergistic effect of
536 manure supply and extreme precipitation on surface water quality. *Environmental Research Letters*
537 **2018**, *13*, (4), 044016.

538 17. Liu, W.; Yang, H.; Ciais, P.; Stamm, C.; Zhao, X.; Williams, J. R.; Abbaspour, K. C.; Schulin,
539 R., Integrative crop-soil-management modeling to assess global phosphorus losses from major crop
540 cultivations. *Global Biogeochemical Cycles* **2018**, *32*, (7), 1074-1086.

541 18. Liu, W.; Yang, H.; Folberth, C.; Wang, X.; Luo, Q.; Schulin, R., Global investigation of
542 impacts of PET methods on simulating crop-water relations for maize. *Agr Forest Meteorol* **2016**, *221*,
543 164-175.

544 19. Hoekstra, A. Y.; Chapagain, A. K.; Aldaya, M. M.; Mekonnen, M. M., *The water footprint
545 assessment manual: Setting the global standard*. Earthscan: London, UK., 2011.

546 20. Franke, N.; Hoekstra, A.; Boyacioglu, H. *Grey water footprint accounting: Tier 1 supporting
547 guidelines*; Value of Water Research Report Series No. 65, UNESCO-IHE, Delft, The Netherlands:
548 2013.

549 21. Liu, W.; Antonelli, M.; Liu, X.; Yang, H., Towards improvement of grey water footprint
550 assessment: With an illustration for global maize cultivation. *J. Cleaner Prod.* **2017**, *147*, 1-9.

551 22. Mekonnen, M. M.; Hoekstra, A. Y., Global anthropogenic phosphorus loads to fresh water and
552 associated grey water footprints and water pollution levels: A high-resolution global study. *Water
553 Resour. Res.* **2017**, *54*, 345-358.

554 23. Williams, J. R.; Jones, C. A.; Dyke, P. T., A modeling approach to determining the
555 relationship between erosion and soil productivity. *Transactions of the Asae* **1984**, *27*, (1), 129-144.

556 24. Liu, W.; Yang, H.; Folberth, C.; Muller, C.; Ciais, P.; Abbaspour, K. C.; Schulin, R.,
557 Achieving high crop yields with low nitrogen emissions in global agricultural input intensification.
558 *Environ Sci Technol* **2018**, *52*, (23), 13782-13791.

559 25. Liu, W.; Yang, H.; Liu, J.; Azevedo, L. B.; Wang, X.; Xu, Z.; Abbaspour, K. C.; Schulin, R.,
560 Global assessment of nitrogen losses and trade-offs with yields from major crop cultivations. *Sci.
561 Total Environ.* **2016**, *572*, 526-537.

562 26. Müller, C.; Elliott, J.; Chryssanthacopoulos, J.; Arneth, A.; Balkovic, J.; Ciais, P.; Deryng, D.;
563 Folberth, C.; Glotter, M.; Hoek, S.; Iizumi, T.; Izaurralde, R. C.; Jones, C.; Khabarov, N.; Lawrence,
564 P.; Liu, W.; Olin, S.; Pugh, T. A. M.; Ray, D.; Reddy, A.; Rosenzweig, C.; Ruane, A. C.; Sakurai, G.;
565 Schmid, E.; Skalsky, R.; Song, C. X.; Wang, X.; de Wit, A.; Yang, H., Global gridded crop model
566 evaluation: benchmarking, skills, deficiencies and implications. *Geoscientific Model Development*
567 **2017**, *10*, 1403-1422.

568 27. Porwollik, V.; Müller, C.; Elliott, J.; Chryssanthacopoulos, J.; Iizumi, T.; Ray, D. K.; Ruane,
569 A. C.; Arneth, A.; Balkovič, J.; Ciais, P.; Deryng, D.; Folberth, C.; Izaurralde, R. C.; Jones, C. D.;
570 Khabarov, N.; Lawrence, P. J.; Liu, W.; Pugh, T. A. M.; Reddy, A.; Sakurai, G.; Schmid, E.; Wang,
571 X.; de Wit, A.; Wu, X., Spatial and temporal uncertainty of crop yield aggregations. *European Journal
572 of Agronomy* **2017**, *88*, 10-21.

573 28. Rosenzweig, C.; Jones, J. W.; Hatfield, J. L.; Ruane, A. C.; Boote, K. J.; Thorburn, P.; Antle,
574 J. M.; Nelson, G. C.; Porter, C.; Janssen, S.; Asseng, S.; Basso, B.; Ewert, F.; Wallach, D.; Baigorria,
575 G.; Winter, J. M., The Agricultural Model Intercomparison and Improvement Project (AgMIP):
576 Protocols and pilot studies. *Agricultural and Forest Meteorology* **2013**, *170*, 166-182.

577 29. Wang, X.; Williams, J. R.; Gassman, P. W.; Baffaut, C.; Izaurralde, R. C.; Jeong, J.; Kiniry, J.
578 R., EPIC and APEX: Model use, calibration, and validation. *T Asae* **2012**, *55*, (4), 1447-1462.

579 30. Gerik, T.; Williams, J. R.; Dagitz, S.; Magre, M.; Meinardus, A.; Steglich, E.; Taylor, R. EPIC
580 User's Manual Version 0810; Blackland research and Extension Center, Temple, Texas, US.: Temple,
581 Texas, 2015.

582 31. Abbaspour, K. C. SWAT-CUP: SWAT calibration and uncertainty programs; Eawag:
583 Dübendorf, Switzerland, 2011.
584

32. Mckay, M. D.; Beckman, R. J.; Conover, W. J., A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. *Technometrics* **1979**, *21*, (2), 239-245.

585

586

33. West, P. C.; Gerber, J. S.; Engstrom, P. M.; Mueller, N. D.; Brauman, K. A.; Carlson, K. M.; Cassidy, E. S.; Johnston, M.; MacDonald, G. K.; Ray, D. K.; Siebert, S., Leverage points for improving global food security and the environment. *Science* **2014**, *345*, (6194), 325-328.

587

588

589

34. Elliott, J.; Müller, C.; Deryng, D.; Chryssanthacopoulos, J.; Boote, K.; Büchner, M.; Foster, I.; Glotter, M.; Heinke, J.; Iizumi, T., The Global Gridded Crop Model Intercomparison: data and modeling protocols for Phase 1 (v1. 0). *Geosci. Model Dev.* **2015**, *8*, (2), 261-277.

590

591

592

35. Fu, J.; Wu, Y.; Wang, Q.; Hu, K.; Wang, S.; Zhou, M.; Hayashi, K.; Wang, H.; Zhan, X.; Jian, Y.; Cai, C.; Song, M.; Liu, K.; Wang, Y.; Zhou, F.; Zhu, J., Importance of subsurface fluxes of water, nitrogen and phosphorus from rice paddy fields relative to surface runoff. *Agricultural Water Management* **2019**, *213*, 627-635.

593

594

595

596

36. Weedon, G. P.; Balsamo, G.; Bellouin, N.; Gomes, S.; Best, M. J.; Viterbo, P., The WFDEI meteorological forcing data set: WATCH Forcing Data methodology applied to ERA-Interim reanalysis data. *Water Resour. Res.* **2014**, *50*, (9), 7505-7514.

597

598

599

37. Mueller, N. D.; Gerber, J. S.; Johnston, M.; Ray, D. K.; Ramankutty, N.; Foley, J. A., Closing yield gaps through nutrient and water management. *Nature* **2012**, *490*, (7419), 254-257.

600

601

38. Portmann, F. T.; Siebert, S.; Doll, P., MIRCA2000-Global monthly irrigated and rainfed crop areas around the year 2000: A new high-resolution data set for agricultural and hydrological modeling. *Global Biogeochem Cy* **2010**, *24*, GB1011.

602

603

604

39. Warszawski, L.; Frieler, K.; Huber, V.; Piontek, F.; Serdeczny, O.; Schewe, J., The Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP): Project framework. *Proceedings of the National Academy of Sciences of the United States of America* **2014**, *111*, (9), 3228-3232.

605

606

607

40. Hempel, S.; Frieler, K.; Warszawski, L.; Schewe, J.; Piontek, F., A trend-preserving bias correction - the ISI-MIP approach. *Earth System Dynamics* **2013**, *4*, (2), 219-236.

608

609

41. McSweeney, C. F.; Jones, R. G., How representative is the spread of climate projections from the 5 CMIP5 GCMs used in ISI-MIP? *Climate Services* **2016**, *1*, 24-29.

610

611

42. Tang, Q. H.; Oki, T.; Kanae, S.; Hu, H. P., The influence of precipitation variability and partial irrigation within grid cells on a hydrological simulation. *J Hydrometeorol* **2007**, *8*, (3), 499-512.

612

613

614

43. Hanasaki, N.; Yoshikawa, S.; Pokhrel, Y.; Kanae, S., A global hydrological simulation to specify the sources of water used by humans. *Hydrol. Earth Syst. Sci.* **2018**, *22*, (1), 789-817.

615

616

44. Wada, Y.; Wisser, D.; Bierkens, M. F. P., Global modeling of withdrawal, allocation and consumptive use of surface water and groundwater resources. *Earth System Dynamics* **2014**, *5*, (1), 15-40.

617

618

619

45. Müller Schmied, H.; Adam, L.; Eisner, S.; Fink, G.; Florke, M.; Kim, H.; Oki, T.; Portmann, F. T.; Reinecke, R.; Riedel, C.; Song, Q.; Zhang, J.; Doll, P., Variations of global and continental water balance components as impacted by climate forcing uncertainty and human water use. *Hydrol Earth Syst Sc* **2016**, *20*, (7), 2877-2898.

620

621

622

623

46. Schewe, J.; Heinke, J.; Gerten, D.; Haddeland, I.; Arnell, N. W.; Clark, D. B.; Dankers, R.; Eisner, S.; Fekete, B. M.; Colon-Gonzalez, F. J.; Gosling, S. N.; Kim, H.; Liu, X.; Masaki, Y.; Portmann, F. T.; Satoh, Y.; Stacke, T.; Tang, Q.; Wada, Y.; Wisser, D.; Albrecht, T.; Frieler, K.; Piontek, F.; Warszawski, L.; Kabat, P., Multimodel assessment of water scarcity under climate change. *Proc Natl Acad Sci U S A* **2014**, *111*, (9), 3245-50.

624

625

626

627

628

47. Mogollón, J. M.; Beusen, A. H. W.; van Grinsven, H. J. M.; Westhoek, H.; Bouwman, A. F., Future agricultural phosphorus demand according to the shared socioeconomic pathways. *Global Environmental Change* **2018**, *50*, 149-163.

629

630

631

48. Bindraban, P. S.; Dimkpa, C. O.; Pandey, R., Exploring phosphorus fertilizers and fertilization strategies for improved human and environmental health. *Biology and Fertility of Soils* **2020**, *56*, (3), 299-317.

632

633

634

49. Zhi, Y.; Yang, Z. F.; Yin, X. N.; Hamilton, P. B.; Zhang, L. J., Using gray water footprint to verify economic sectors' consumption of assimilative capacity in a river basin: model and a case study in the Haihe River Basin, China. *J Clean Prod* **2015**, *92*, 267-273.

635

636

637
638

50. Ward, M. H.; Jones, R. R.; Brender, J. D.; de Kok, T. M.; Weyer, P. J.; Nolan, B. T.;

639

Villanueva, C. M.; van Breda, S. G., Drinking water nitrate and human health: An updated review. *Int*

640

*J Environ Res Public Health* **2018**, *15*, (7), 1557.

641

51. Sinha, E.; Michalak, A. M.; Balaji, V., Eutrophication will increase during the 21st century as

642

a result of precipitation changes. *Science* **2017**, *357*, (6349), 405-408.

643

52. Carpenter, S. R.; Bennett, E. M., Reconsideration of the planetary boundary for phosphorus.

644

*Environ Res Lett* **2011**, *6*, (1), 014009.

645

53. Rockström, J.; Steffen, W.; Noone, K.; Persson, A.; Chapin, F. S.; Lambin, E. F.; Lenton, T.

646

M.; Scheffer, M.; Folke, C.; Schellnhuber, H. J.; Nykvist, B.; de Wit, C. A.; Hughes, T.; van der

647

Leeuw, S.; Rodhe, H.; Sörlin, S.; Snyder, P. K.; Costanza, R.; Svedin, U.; Falkenmark, M.; Karlberg,

648

L.; Corell, R. W.; Fabry, V. J.; Hansen, J.; Walker, B.; Liverman, D.; Richardson, K.; Crutzen, P.;

649

Foley, J. A., A safe operating space for humanity. *Nature* **2009**, *461*, 472-475.

650

54. Liu, W.; Yang, H.; Liu, Y.; Kummu, M.; Hoekstra, A. Y.; Liu, J.; Schulin, R., Water resources

651

conservation and nitrogen pollution reduction under global food trade and agricultural intensification.

652

*Sci. Total Environ.* **2018**, *633*, 1591-1601.

653

55. Lassaletta, L.; Billen, G.; Garnier, J.; Bouwman, L.; Velazquez, E.; Mueller, N. D.; Gerber, J.

654

S., Nitrogen use in the global food system: past trends and future trajectories of agronomic

655

performance, pollution, trade, and dietary demand. *Environ Res Lett* **2016**, *11*, (9), 095007.

656

56. Powers, S. M.; Chowdhury, R. B.; MacDonald, G. K.; Metson, G. S.; Beusen, A. H. W.;

657

Bouwman, A. F.; Hampton, S. E.; Mayer, B. K.; McCrackin, M. L.; Vaccari, D. A., Global

658

opportunities to increase agricultural independence through phosphorus recycling. *Earth's Future*

659

**2019**, *7*, (4), 370-383.

660

57. Wang, Q.; Zhou, F.; Shang, Z.; Ciais, P.; Winiwarter, W.; Jackson, R. B.; Tubiello, F. N.;

661

Janssens-Maenhout, G.; Tian, H.; Cui, X.; Canadell, J. G.; Piao, S.; Tao, S., Data-driven estimates of

662

global nitrous oxide emissions from croplands. *National Science Review* **2020**, *7*, (2), 441-452.

663

58. FAO *The future of food and agriculture – Alternative pathways to 2050*; Rome, 2018; p 224.

664

59. Folberth, C.; Gaiser, T.; Abbaspour, K. C.; Schulin, R.; Yang, H., Regionalization of a large-

665

scale crop growth model for sub-Saharan Africa: Model setup, evaluation, and estimation of maize

666

yields. *Agr Ecosyst Environ* **2012**, *151*, 21-33.

667

60. Van Vuuren, D. P.; Bouwman, A. F.; Beusen, A. H. W., Phosphorus demand for the 1970-

668

2100 period: A scenario analysis of resource depletion. *Global Environ Chang* **2010**, *20*, (3), 428-439.

669