Introduction

Over the past decades, the extensive application of pesticides has unquestionably contributed to meeting the demands of global food production1,2. However, agrochemicals, can cause adverse effects on non-target species and environments3,4 as well as the health of human beings5, and may even ultimately pose a threat to the Sustainable Development Goals (SDGs) set by the United Nations6. Nearly 75% of the world’s farmland has been estimated to be confronted with risks from pesticide use7. China has become the largest user of pesticides in the world, with about 43% of global pesticide consumption on <9% of global cropland8. Specifically, about 1.77 million tonnes of agricultural pesticides were applied in China, constituting about 42% of the total use globally in 20199, and the average annual usage in China between 1990–2019 has been ~1.42 million tonnes. To fully understand the likely effects of the pesticide risks under high usage, it is essential to quantify the concentrations that occur in different environments. Besides, ecological risks from pesticides which may hinder the achievement of SDGs set by the United Nations10 have attracted widespread attention internationally7,11,12,13. Challenges remain in China to achieve some of its SDGs, such as Target 2 (Zero hunger), Target 3 (Good health and well-being), and Target 6 (Clean water and sanitation) which are related to pesticide risk, according to the United Nations 2022 Sustainable Development Report14. However, a thorough analysis on the use, budget, fate, and associated risk effects of pesticides for informing SDGs is still lacking. Comparison of the available monitoring data on pesticides is significantly hindered by the inconsistent ingredients, analytical techniques or environmental compartments over a wide time scale. Besides, the diversity and trace concentrations of pesticides, as well as different transfer pathways and changes in amounts of pesticides over time, make it quite time and cost consuming to quantify fully at the large scale. Multimedia models are used to predict the transport and fate of chemicals for enactment of policies15, using readily acquired data of environmental attributes and chemical application. Extensive efforts have been made to quantify the specific fate of pesticides in different environmental compartments as well as the potential risks16,17,18. However, challenges include the presence of pesticide legacy, mixtures, and varying receptors and end-points. Residual pesticides retained in the environment, especially in soils and sediments, always act as the background value for the next application, and contribute to the pollution load next year, potentially leading to an increase in legacy risk19. Typically, the ecological risks of pesticides depend on their toxicity to organisms20,21. Thus, additive toxic effects caused by pesticide mixtures are getting attention, while consideration about associated risk is inadequate22. Ignoring these processes will lead to errors in assessing risks. Thus, a systematic analysis of the spatiotemporal pesticide risk is needed.

The hotspots where pesticides pose high threats to biodiversity, human health, and water resource, which are associated with SDGs, should be explored for tailored decision-making. Examples for the influence of certain pesticides on biodiversity include atrazine effects on amphibians23, and neonicotinoid insecticides on bees24,25. Vast areas with high pesticide risk may have potential influence on a variety of species, which partly depend on the number of risk receptors and exposure pathways. Environmental exposure to pesticides may increase the incidence of certain neurological disorders affecting human health5, and under the same pesticide mass, more densely populated regions would be exposed to higher health risk11. A large percentage of the China’s population is at water security risk, and poor water quality exacerbates water scarcity. Thus, those regions need to be identified to achieve China’s SDGs by alleviating biodiversity loss, human health risk, and water scarcity.

To address the knowledge gaps, here we: i. conducted a national scale and long-term fate analysis for four most widely-used pesticides, using national detailed data on usage, combined with a fully parameterized fugacity model for major river basins; (ii) addressed changing ecological risk for pesticide mixtures in three environmental compartments (namely soil, surface water, and sediment); (iii) identified vulnerable regions of the country by juxtaposing ecological risk with population, biodiversity, and water resource risk; and iv. analyzed future directions for sustainable development in China, based on several scenarios of public policies.

Results and discussion

Target pesticides in the multimedia environment and their spatiotemporal fate

Firstly, the fluxes for target pesticides in the multimedia environment over 30 years were predicted. Carbofuran, chlorpyrifos, atrazine and acetochlor, were selected as representatives considering their frequent and widespread applications, high ecological toxicity and complete ecotoxicological databases. The total amount of active ingredients (AIs) of the target pesticides applied in the 79 river basins of China (Supplementary Fig. 1) in 2017 was estimated as 2.38 × 105 tonnes. Yearly use leads to the legacy of pesticides that represents the overall amount of pesticide residues from the previous residues of last year and the newly applied amount in the current year. 2.24 × 104 tonnes of target AIs were estimated to be retained in the environment at the end of 2017, most of which were in soils (1.98 × 104 tonnes, 88%) and water bodies (2.56 × 103 tonnes, 11%) (Fig. 1). Details about the fluxes and legacy of target pesticides in major river basins are shown in Fig. 2 and Supplementary Table 1.

Fig. 1: Schematic diagram of transport and legacy of target pesticides in multimedia.
figure 1

Transmission rate is calculated as the overall rate of target pesticides in all basin units at the end of 2017. Legacy of four target pesticides per year in soil, water and sediment was shown in a, b, and c respectively.

Full size image
Fig. 2: Estimated pesticides usage, legacy and advection into seas for major river basins in a recent year.
figure 2

The number in each river basin represents the legacy in 2017.

Full size image

Overall, the major transfer pathways of the pesticides included: transfer through soil to nearby river networks (2.39 × 104 tonnes yr−1 (tonnes per year), with atrazine and carbofuran accounting for 51% and 31%, respectively), soil degradation (2.28 × 105 tonnes yr−1, with chlorpyrifos and carbofuran accounting for 39% and 23%, respectively), as well as water advection and degradation (2.30 × 104 tonnes yr−1, with atrazine and carbofuran accounting for 52% and 31%, respectively). Degradation in soil, water and sediments represents the main removal for carbofuran, chlorpyrifos, and acetochlor, while 12% of the atrazine loss was through advective flow into the ocean, due to its relatively long half-life in water phase (2064 h, Supplementary Table 2). The discharge of the target pesticides into the ocean from the six major rivers was estimated as 6.59 × 103 tonnes yr−1, with the Huaihe River accounting for the largest proportion (36%, 2.35 × 103 tonnes yr−1). Atrazine accounted for the largest proportion of the target pesticides for each basin, with regards to the water advection. The total amount of target pesticides from the Liaohe, Haihe and Yellow River into the Bohai Sea was about 2.97 × 103 tonnes yr−1 with atrazine accounting for 79%. For comparison, atrazine fluxes were 1.09 × 103 kg yr−1 from the Scheldt Estuary into the North Sea26 and 0.72 × 103 kg yr−1 from Ebro River to Western Mediterranean Sea27. The marginal seas receive a large amount of pesticides and the land-ocean interface would be the hotspots, such as the Yellow River Delta and the Bohai Sea28.

The environmental legacy of pesticides is projected to have increased rapidly from 1991 to 1999, then slowed between 1999–2015. Then the pesticide legacy began to decrease after 2015 (Fig. 1), which is consistent with pesticide historical usage (Supplementary Fig. 2). The recent downward usage in China can be attributed to the “zero growth in pesticide use” policy in China since 2015. Considering the compound-specific scenarios, chlorpyrifos in soil increased most (average 300 tonnes yr−1) from 1991 with the peak (about 11,000 tonnes) in 2014. It is estimated to have decreased by about 30% by 2020. As for atrazine, its legacy and growth rate were the largest in water and sediment, with the growth by 66 tonnes yr−1 and 0.58 tonnes yr−1 respectively from 1991 to 2014. The legacy was stable over the studied period, with the multi-year average in soil, water and sediment being 2563 tonnes, 148 tonnes, 0.26 tonnes for carbofuran and 1841 tonnes, 320 tonnes, 1.25 tonnes for acetochlor respectively. Chlorpyrifos shows the highest accumulation effect (defined as the ratio of the previous legacy to the new application in the following year), with the peak (13.5%) in 2018. Although the accumulation contributed only 7.1–9.2% to the input due to its short half-life, the concentrations would be underestimated if not considering the accumulation effects for long-term risks (Supplementary Table 2). This underestimation would become significant over time due to continous accumulation of pesticides.

The concentrations of target pesticides in each basin simulated for different compartments in 2017 are given in Supplementary Table 3. The predicted concentrations in air were negligible, as these compounds are involatile. The average concentrations in the soil phase ranked as: chlorpyrifos (10 ng g−1) > atrazine (4.6 ng g−1) > carbofuran (3.1 ng g−1) > acetochlor (2.4 ng g−1), which is mainly related to the usage and properties of pesticides and compartments. The highest chlorpyrifos concentration was found in the Tuhai-Majia River basin (ID: 18, 112 ng g-1), followed by Eastern Shandong River Basins (ID: 31, 94 ng g−1), with the lowest in the Northwest River basin. Atrazine had the highest legacy in water and sediment, largely due to its long half-life in these compartments (2064 h and 1920 h respectively). Areas with aquatic atrazine over 3 × 103 ng L−1 were Eastern Shandong Rivers (ID: 31), the Yellow River downstream of Huayuankou (ID: 25) and the Tuhai-Majia River (ID: 18), which are mainly located in coastal regions in Shandong province. The lowest concentrations was found in Northwest and Southwest River basins. In the Yangtze River, the highest concentration of Atrazine occurred in the middle and lower reaches (ID: 35, 570 ng L−1; ID: 37, 480 ng L−1; ID: 42, 384 ng L−1), indicating more attention is needed on the downstream ecological risk of pesticides over large river basins.

Pesticide risk and vulnerable regions

Secondly, the Risk Score (RS)7 was used to evaluate the ecological harm caused by mixed pesticides, representing the level of the compartment with the highest ecological risk of mixed pesticides that was calculated from the Risk Quotient (RQ)29. The RQ over 30 years (Supplementary Fig. 3, Supplementary Table 4) showed that the proportion of moderate and high ecological risk areas increased since 1991 and then decreased in recent years. Acetochlor posed the biggest threat to organisms in the aquatic environment, especially algae, but chlorpyrifos risk to Folsomia candida in soil with consistently over 50% in high-risk river basins, followed by carbofuran. Subsequently, the RS was obtained for a comprehensive assessment, considering pesticide mixtures and the exposed compartment. Soil was identified as the compartment with highest risk in most areas, while 16% of river basins had highest risk in surface waters, especially in Tibet, Xinjiang, and the loess Plateau region.

As for spatiotemporal pattern, the proportion of risk areas (RS > 0) expanded from 47% in 1991 to 79% in 2014, when the usage peaked (Supplementary Fig. 4). The proportion of high-risk areas (RS > 1) increased from 11% in 1991 to 30% in 2014, with two sub-basins of the Haihe and Huaihe basins (ID: 18, 31) classified as extremely high (RS > 2). By 2017, 79% of China were identified at ecological risk (RS > 0), characterized by an uneven distribution. Thirty-two river basins were identified as high risk (RS > 1) (29.8% in area), located mainly in Eastern, Southern, and Central China (Fig. 3a). This was apparently consistent with the distribution of the main grain-producing and populated areas (nearly 68% of population), revealing a spatial relation with the well-known geographical ‘Hu-Huanyong line’30. Chlorpyrifos in soil was the main cause of the high RS as it was mainly applied on rice and wheat in south and north of China. Regional analysis showed that high-risk areas were located in the Huaihe River and Shandong coastal basin (ID: 27, 28, 29, 30, 31), Haihe River basin (ID: 15, 16, 17, 18), parts of Songhua River basin (ID: 3, 6), Liaohe River basin (ID: 10, 11), the middle and downstream reaches of Yellow River basin and Yangtze River basin (ID: 23, 24, 25, 38, 39, 40, 41, 42, 43), Southeastern Rivers (ID: 44, 45, 49), Pearl River and the south coastal rivers (ID: 54, 55, 56, 57, 58, 59). Most were among the large crop producers in China; i.e., the Huaihe River Basin is the main spring and summer corn sowing area, as well as cotton, soybean and wheat production area. The Yangtze River Basin, especially the middle and lower reaches, is the main rice producing area in China, accounting for about 70% of the domestic yield31. What calls for special attention is the Huaihe River and Shandong coastal river basin (ID: 27, 28, 29, 30, 31), and Haihe River basin (ID: 15, 16, 17, 18), which were identified as high-risk areas, characterized by an abundance of fertile farmlands, large industrial cities and high population density. China’s agricultural population accounted for 41% in 201732, and the large-scale intensive farming, the great contributors of pesticide use, is also concentrated in economically developed areas. Of particular concern is the Tuhai-Majia River (ID: 18) located in Shandong province. Its high RS (RS = 2.06) was mainly attributed to the applications and properties of chlorpyrifos in soil (RQs = 112), which is consistent with previous studies reporting serious pollution of the Haihe River basin12,33. More attention should be paid here because of its intense human activities and proximity to the Bohai Sea.

Fig. 3: Pesticide risk in vulnerable regions considering population, biodiversity and water resource risk.
figure 3

79 river basins were the basic simulation units. a is the national map of pesticide risk score (RS) (RS < −1, negligible; −1 ≤ RS < 0, low; 0 ≤ RS < 1, medium; 1 ≤ RS < 2, high; RS ≥ 2, extremely high). RS ≤ 0 corresponds to <5% probability for random species to be affected by pesticides 7. e, f, and g are maps of concern level derived when the population density (b), biodiversity (species) (c), and water resource risk (d) are overlain against the pesticide risk map (a), respectively. The concern level classification is shown in Supplementary Table 4.

Full size image

To identify vulnerable regions considering multiple targets, pesticide risk was then integrated with the distribution of population density (Fig. 3b), biodiversity (Fig. 3c), and water scarcity (Fig. 3d), respectively. It was found that hotspots were mainly located in mid-eastern and southeastern China, such as the Huaihe River basin (ID: 28, 29, 30, 31), the Haihe River basin (ID: 16, 17, 18), lower reaches of the Yellow River, Yangtze River, and Pearl River basin (ID: 23, 25, 42, 43, 55, 56, 57, 58) etc. Specifically, the Tuhai-Majia River (ID: 18) within the Haihe River basin was classified as the key vulnerable region for sustainable development, with high pesticide risk, water scarcity, and population density.

Regions with high population or biodiversity providing essential ecosystem services are potentially more subject to negative influences from excessive pesticide use21. Results show that 31% of the national high-risk areas (RS > 1) were located in regions with high population density (>471 km−2), while 3% of which (RS > 1) also bear high biodiversity (>781 tetrapod species km−2). These hotspots are mainly concentrated in parts of Huaihe (high population-pesticide risk) (ID: 28, 29, 30, 31), Haihe (high population-pesticide risk) (ID: 17, 18) and Pearl River basins (ID: 55, 56, 57, 58), especially the Dongjiang delta (ID: 56) (high biodiversity-population-pesticide risk). It was reported that regional biodiversity was at risk globally, owing to elevated pesticide contamination in surface waters3,4. This threat may escalate taking other compartments (soil and sediments) and effects of mixture toxicity into consideration. Hotspots of biodiversity and high-risk pesticide did not overlap substantially, but the overlapping area might expand with increased farmland and population growth in the future, which may hinder the achievement of the “Life on land” target of China’s SDGs.

Transition of pesticides through runoff and infiltration might constitute a threat to water resource security with surface water being the susceptible compartment7,12. Our analysis showed that about 7% of mainland China suffered from limited water resource (water risk score > 3.5) and high pesticide risk (RS > 1), especially populated regions (16% of China’s population) in northern China. Pesticide-related deterioration may intensify the pressures between water supply and demand, owing to low few safe water in risk regions, which would impede the achievement of “Clean water & Sanitation” target of China’s SDGs.

New paradigm for pesticide management

The legacy and risk of pesticides under different scenarios were analyzed to inform future management. Under scenario 1 (S1, pesticide applications continue declining by 8% after 2020), China would be risk-free (RS < 0) at the national scale around 2075 (max RS is below 0). Overall legacy of target pesticides in soil, water and sediments would be 158 tonnes, 21.7 tonnes, and 0.16 tonnes respectively in 2075. Time series analysis showed that the maximum RS would be 0.99 (exactly high-risk-free) around 2047 (Supplementary Table 5). This scenario can be regarded as a baseline, indicating it will take a long time to become risk-free nationwide, although the pesticides usage had declined since 2014. Thus, current reductions may not be effective immediately due to the large previous legacy. This estimation may be conservative as there are emerging pesticides residues with potential risk in the future, given the agriculture’s dependence on pesticides.

Legacy would significantly reduce in 2025 if the utilization rate reaches 55% as the level of developed countries (scenario 3, S3) but slightly for S2 explained in Methods (Supplementary Table 6). The removal of pesticide legacy would be 22%, 24% and 20% in soil, water and sediments respectively in 2025 compared with 2020. Under this condition, the risk area (RS > 0) and the high-risk area (RS > 1) would reduce by 8% and 4%, respectively, compared with those in 2020.

In addition to application efficiency, efforts have also been made to develop lower toxic agrochemicals to reduce the possible risk for non-target organisms especially for humans34. With pesticide substitution (S4), the risk regions (RS > 0) would decrease immediately from 79% to 72% in 2017 and the proportion of high-risk regions would decline from 30% to 14% and high risk areas would be concentrated in the northern China (Fig. 4). Although the risk areas (RS > 0) decreased slightly from 75 to 73% after pesticide substitution, it is surprising that the ratio of high-risk regions (RS > 1) in 2020 increased from 18% to 44%, especially in northeastern, central, and southeastern China. This unexpected increase is related to characteristics of the substituting pesticides. Several selected substitutions with low toxicity to humans have a longer half-life, which would result in risk transfer among multimedia environments and their effects can be amplified after long-term accumulation. Although the alternative pesticides have lower toxicity to humans, their ecotoxicity may not be correspondingly low, such as low No Observed Effect Concentrations (NOEC) to Daphnia magna for diflubenzuron compared with carbofuran. If not considering this tradeoff, it may come at the expense of another goal, such as “Life on land” target, while meeting one human-health related goal. Thus, pesticides substitution should comprehensively consider the changes in usage, transmission in multimedia as well as accumulation and ecotoxicity of pesticides before implementation, in order to achieve coordination among several SDGs10.

Fig. 4: Spatial distribution of risk in 2017 (a) and 2020 (b) under substitution scenario (S4).
figure 4

The target pesticides were assumed to be substituted by diflubenzuron, chlorantraniliprole, metolachlor, and sethoxydim respectively from 2017.

Full size image

Our results show that a large part of China was still subject to pesticide risks, including several regions with high population, biodiversity, and water scarcity, which may require more tailored strategies with the purpose of minimal damage, although the legacy of pesticides has started decreasing over recent years. Farmland covers about 38% of the land surface globally35 and the amount of pesticides usage may increase due to further intensified agriculture with more population in a warmer future3,36,37. Additionally, the water resources would be more unevenly distributed with climatic warming; thus, there might be more vulnerable regions facing the double pressures from water scarcity and pesticide risk. Meeting the needs of China’s growing population (“Zero Hunger”) by applying more pesticides cannot be a priority, as it may influence other SDGs such as “Clean Water and Sanitation”, “Life Below Water” and “Life on land”. Approaches reducing pesticide risks while ensuring food security are required, so as to achieve sustainable agricultural development.

Source management is still of top concern for reducing risks given the strong correlation between pesticides usage and legacy. Solutions should focus on critical regions with high application, such as Shandong, Henan, and Hunan province in China. Decision support systems (DSSs)38 which may reduce the frequency of pesticide application, can be adopted in these regions to build sustainable agricultural systems. Moving from generic targets for reductions in pesticide inputs towards targeted mitigation actions requires identifying regional hotspots for pesticide-related impacts. Some vulnerable regions such as Huaihe River, Haihe River, Pearl River, the middle and lower reaches of Yellow River and Yangtze River can be set as the representative areas. Meanwhile, priority measures should be focused on highly toxic pesticides such as high-risk chlorpyrifos in soil as identified in our study. Certainly, new pesticides should be of concern, such as neonicotinoids25 and synthetic aldehyde pesticides39. To ensure a precise and effective control, it is suggested that detailed and frequent field investigation should be conducted in certain vulnerable regions before tailored mitigation measures can be implemented at the regional scale. Integrated pest management such as green prevention, instead of using pesticides, should be considered21,40. Organic farming16 and precision techniques41 can also be integrated with other disease management approaches. As pesticide risk would constitute an emerging threat to sustainable development, the results herein would be vital towards future risk assessment and mitigation efforts in China.

Limitations and future suggestions

Some limitations in this study were identified. Firstly, it is assumed that pesticides were directly applied into soil phase, without considering possible ways into surface water or air firstly through spraying, irrigating or washing, resulting in overestimation of the pesticides in soil phase. The spray drift loss influenced by the application methods, crops, pesticide types, and weather condition42, would be difficult to consider at a large scale7. The uncertainty analysis indicated that lack of knowledge of spray drift loss would cause about 21–78% overestimation of soil pesticide content, based on the spray drift ratio ranging from 20% to 50% over the past few decades42,43. Indeed, advanced technologies are being developed to reduce pesticide spray drift44 and there is a global shift toward seed coatings rather than aerial spraying45. The sampling survey in this study showed that about 89% of fields receiving carbofuran were treated with coated corn seeds which can greatly reduce drift of pesticides application. Besides, the spray drift always occurs on farmland and thus, most spray residues will end up in the treated/adjacent soils directly or indirectly. Thus, the soil pesticide content may have been overestimated, but within a reasonable range. The results in soil phase would be also underestimated because they were regarded without clear distinction between treated and untreated soils. Thus, the results of this study could be used at least as references for evaluating the pesticide risk at large spatiotemporal scale. Secondly, the ignorance of pesticide application in earlier years may lead to underestimation of the legacy stocks of pesticides, while the application of pesticides was relatively low with slight cumulative effects before 1991. Thirdly, this study assumes that there is no interception of pesticides by reservoirs and dams, and siltation of sediment in rivers. Thus, the discharge of the target pesticides into the ocean was overestimated. Fourthly, no account was taken of synergistic effects of pesticide mixtures which were difficult to quantify. It is still difficult to conduct precise dynamic simulation of chemicals multimedia fate at large scale, mainly resulting from the lack of detailed data.

Overall, a refined national database is recommended to be developed in future research, such as high precision spatiotemporal distribution of crop cultivation and pesticide application by watershed. The global database of the top 20 crop-specific pesticide application rates from 2015 to 2025, PEST-CHEMGRIDS46, was produced based on the data from the USGS Pesticide National Synthesis Project47 and the country-specific pesticide use data reported by FAOSTAT9, which contains three of the four target pesticide use rates in China relevant in this study (except carbofuran). Compared to the local fields survey, the inference of pesticide application rates from the USA to global scales and the map resolution of 10 km do not support refined national pesticide risk analysis for China. Nonetheless, yearly usage of the target pesticides from 1991 to 2020 in this study was extrapolated by the same application ratio and AIs obtained from the sampling survey in 2017. More factors should be considered during extrapolation like the method of PEST-CHEMGRIDS46, such as hydroclimatic and agricultural variables, policy and socio-economic indices. This also needs follow-up investigation at a more refined local level, which will reduce the uncertainty of large-scale assessment.

Methods

Creation of detailed pesticide database

In this study, mainland China was divided into 79 river basins according to the standard ‘Code for China river name (SL 249–1999)’ as the simulation units. Taiwan, Hong Kong, Macao and other islands, were not considered due to lack of available data (Supplementary Table 7). Spatial data, including digital elevation and land use maps, generated at a resolution of 1 km × 1 km, were acquired from the Resource and Environment Science and Data Center of the Chinese Academy of Sciences48. Statistical data such as the total amount of pesticides applied yearly in each province were compiled from the China Rural Statistical Yearbook and Statistical Yearbook from the National Bureau of Statistics of China, but without detailed usage by pesticide types and crops.

Additionally, a more detailed database was established, including information on application ratio, crops, seasons, intensity, and AIs contents of target pesticides types (Supplementary Fig. 5, Supplementary Tables 8 and 9). The database was based on a comprehensive sampling survey conducted in about 20,000 fields distributed in the districts and counties with representative planting patterns and production modes in each province across the country in 2017, in conjunction with the Academy of Agricultural Sciences of China. These data are provided in dataset file. Based on the investigations, carbofuran, chlorpyrifos, atrazine and acetochlor were selected as representative pesticides considering their frequency of application, high ecological risk, and complete ecotoxicological databases. All four pesticides usually exist in the same period due to their non-overlapping functions, application on major crops throughout the country according to field survey. Carbofuran (carbamate insecticide) was widely used on up to 80 crops for pest control, mainly in northeast, north China, the middle and lower reaches of the Yangtze River and coastal areas. Chlorpyrifos (organophosphate insecticide) was selected due to its relatively high ecotoxicity to aquatic life and bees, strong bioaccumulation and amplification, and the role of the preferred products of organophosphate insecticides in recent years. Atrazine (triazine herbicide) has a stable structure and very slowly biodegradable with long half-life of 696 h in soil and 2064 h in water phase. Acetochlor, the most widely used chloramide herbicide, has environmental hormone effects, and strong acute toxicity to lower vertebrates, plankton, and small and medium annelids; it has been limited by the European Commission and the USEPA. Details about four pesticides can be found in Supplementary Table 2. The mixture of the 4 compounds is used here to assess threats to flora and fauna, especially in aquatic environments.

This AI database was then used as a reference throughout the studied period, including the ratio of yearly-averaged application and averaged AIs contents for a specific pesticide, to retrieve and revise the missing datum in other years. Hence, yearly usage of the target pesticides in river basins from 1991 to 2020 was calculated by the national public pesticide application data integrated with the detailed AI database, as shown in Eq. (1). Interpolation methods were used where data was unavailable. Due to the lack of data in 2020, the application amount was estimated by using the average growth rate of total pesticide application for the last available three years.

$${\rm{Pe}}{{\rm{s}}_{{\rm{type}}\,{\rm{i}},\,{\rm{bas}}}} = {\rm{Are}}{{\rm{a}}_{{\rm{bas}}}}\sum\limits_{{\rm{s}} = 1}^{\rm{n}} {\frac{{{\rm{Pe}}{{\rm{s}}_{{\rm{pro}}\,{\rm{s}},\,{\rm{total}}}}}}{{{\rm{Are}}{{\rm{a}}_{{\rm{pros}}}}}}} \cdot {{\rm{m}}_{\rm{i}}} \cdot {{\rm{m}}_{\rm{s}}}$$
(1)

Where Pestype i, bas represents the amount of each pesticide i applied in each study basin (kg); Areabas and Areapro s represent the total cropland area of each basin and each province s in the basin, respectively (m2); Pespro s, total represents the total amount of pesticide applied in each province (kg); n represents the number of provinces in the basin; mi represents the proportion of each pesticide to the total amount of pesticide applied in the province, and ms represents the proportion of cropland area in province to the total cropland area in this basin that was obtained by the spatial analysis in ArcGIS 10.2.

Transfer processes and state simulation

Modeling and analyses were conducted using our modified Level III fugacity model, Matlab R2020a, ArcGIS software (v.10.2), and Origin 2018. The original fugacity model was based on a steady-state, non-equilibrium, and flowing system49, including four major bulk compartments (air, water, soil, and sediment). The mass balance equation was established as Eq. 2 and detailed transfer processes were defined in Supplementary Table 10. To achieve cross-region and long-term study, the model was modified with particular consideration on transmission among spatial units and accumulation over time. Based on the river map, transmission process of pollutants, referring to the advection transmission in water phase, is presented by inputting the upstream advection transmission results into the downstream unit as boundary conditions. Specifically, the advective flow of pesticides in water phase of the river basin (T02t) was the sum of all upstream inflows (T20t). The residues of pesticides in soil and sediment compartments in year i were calculated in the form of annual average background inputs (T03b, T04b) in year i+1 as initial conditions. These functions were implemented via MATLAB R2020a.

$$\begin{array}{*{20}{l}} {{{{\mathrm{Air}}}}:T_{01{{{{t}}}}} + T_{21d} + T_{31d} = T_{10t} + T_{12d} + T_{12p} + T_{12w} + T_{13d}}\\\qquad\; +\, T_{13p} + T_{13w} + T_{10m} \\ {Water:T_{02t} + T_{12d} + T_{12p} + T_{12w} + T_{32e} + T_{32l} + T_{42d} + T_{42r}}\\\qquad\quad\;\; = T_{20t} + T_{21d} + T_{24d} + T_{24s} + T_{20m} \\ {{{{\mathrm{Soil}}}}:T_{03e} + T_{03b} + T_{13d} + T_{13p} + T_{13w} = T_{30m} + T_{31d} + T_{32e} + T_{32l}} \\ {{{{\mathrm{Sediment}}}}:T_{04b} + T_{24d} + T_{24s} = T_{40m} + T_{42d}} \end{array}$$
(2)

To obtain the flux (T) and concentration of pesticides (C), it is important to clarify the relationship among those parameters by the fugacity method (Eq. 3).

$$T_{{{{\mathrm{in}}}} - {{{\mathrm{out}}}}} = G \times C,\,T_{{{{\mathrm{trans}}}}} = D \times f,\,C = f \times Z$$
(3)

Where G is the flow rate (m3 h−1), f is fugacity (Pa). Z-value (mol m−3 Pa−1) represents the capacity of a certain medium to contain pollutants at a given fugacity (Supplementary Table 11). Besides, D value (mol h−1 Pa−1) describes four aspects of environmental processes: advection, reaction, diffusion, and interface processes (Supplementary Table 12). The fluxes to the ocean were computed considering the advective transport of the major seaward rivers. These parameters were calculated by the environmental attributes, transport parameters (Supplementary Table 13), and pesticide properties (Supplementary Table 2) as well as application rate based on above investigation. The changing use and previous legacy of four target pesticides over time were used as the input for a dynamic model. The regionalized parameters among 79 river basins mainly included area of four compartments and organic carbon content of soil phase (Supplementary Table 14). The area of environmental media was obtained from the land use data at a resolution of 1 km × 1 km by the Data Center for Resources and Environmental Sciences, Chinese Academy of Science48. The tabulated and georeferenced soil organic carbon content was obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Science50, including 28 soil property parameters in >1600 sampling sites around the country. These data were processed in ArcGIS 10.2. to generate basin-specific data.

With this approach, the Predicted Environmental Concentrations (PECs), transfer fluxes, and mass inventories of target pesticides in each basin were obtained. Sensitivity analysis was conducted by changing each parameter by ±10%51.

$${{{\mathrm{SC}}}}_i = \frac{{\Delta Y_i/Y_i}}{{\Delta X_i/X_i}}$$
(4)

Where SCi represents the sensitivity coefficient for the ith input parameter. ∆Xi / Xi and ∆Yi / Yi represent the ith input parameter and relative rate of change for corresponding output parameter. In common, |SC| > 0.6 represents significant sensitivity, 0.2 < |SC| < 0.6 is moderate sensitivity, while |SC| < 0.2 is low sensitivity.

A total of 19 parameters with sensitivity coefficient > 0.6 were reasonably adjusted. Of these parameters, some associated with characteristics of the study region were adjusted from public national data to local data with high precision, such as the organic carbon content of soil and sediment. Those parameters with low sensitivity like diffusion coefficients were taken as the model default values. The simulated PECs were compared to the measured data from typical regions reported in literature, indicating a reliable simulation of pesticide fate (Supplementary Figs. 6 and 7). The distribution trends between predicted and measured pesticide concentrations were generally consistent, especially for the four target pesticides in water. The overall differences were not significant and within the acceptable range.

Uncertainty from spray drift during application was most likely to lead to overestimation of soil pesticide content, as shown by the following: we selected the Tuhai-Majia River basin (ID: 18) with the highest soil chlorpyrifos content in 2017 as an example for highest uncertainty analysis. The sampling survey in 2017 showed that about 84% fields applied chlorpyrifos were treated with foliar spray that caused more drift loss than ground spraying which is mostly used for the other three target pesticides. Chlorpyrifos also has the highest vapor pressure among the target pesticides (Supplementary Table 2). The drift loss ratio from the applied amount was estimated to be from 20% to 50%42,43. We assumed that 20–50% applied chlorpyrifos was lost to the air compartment and 50–80% was applied directly into the soil compartment, and tested the overestimation of soil pesticide content due to lack of knowledge of spray drift loss.

Ecological risk assessment

The RQ29 method was used for ecological risk assessment for composite chemicals13. The RQ is defined as the ratio between PEC (Predicted Environmental Concentration) and PNEC (Predicted No Effect Concentration) (Eq. 5).

$$RQ_i^k = {{{\mathrm{PEC}}}}_i^k/{{{\mathrm{PNEC}}}}_i^k$$
(5)

Where the subscripts of i, k were target pesticides and environmental compartments, respectively. The PEC in this study refers to the predicted concentration of certain AIs in each compartment after application considering the transport and accumulation through the fugacity modeling. The PNEC was derived from NOEC from toxicological data in PPDB52 and ECOTOX53 divided by the assessment factor (AF).

$${{{\mathrm{PNEC}}}}_i^k = {{{\mathrm{NOEC}}}}/{{{\mathrm{AF}}}}$$
(6)

The AF approach was derived based on the amount of available chemical toxicity data54, initially for aquatic organisms and applied to soil organisms as well (Eq. 6).

The RQ in sediments (RQsed) was calculated using Eqs. 7 and 8 according to Hu et al.55:

$${{{\mathrm{RQ}}}}_{{{{\mathrm{sed}}}}} = {{{\mathrm{PEC}}}}_{{{{\mathrm{sed}}}}}/{{{\mathrm{PNEC}}}}_{{{{\mathrm{sed}}}}}$$
(7)
$$PNEC_{sed} = PNEC_w \times K_{oc} \times f_{oc}$$
(8)

where PECsed is predicted concentration in dry sediment; PNECsed and PNECw are predicted no effect concentrations in sediment and water, respectively; Koc is organic carbon-water partition coefficient of for the pesticide under consideration, and foc is the organic carbon content in solids of bulk sediment. The above environmental and pesticide properties are given in the database (Supplementary Table 2 and 13).

In this study, the hierarchical approach of the Pesticide Use Risk Evaluation (PURE) decision-support system (DSS) was adopted56. The risk point (RPk) for compartment k is calculated as the log-transformed sum of all RQ in the corresponding compartment (Eq. 9). The RS in each basin division unit was defined as the maximum RP across four compartments (Eq. 10). RS was classified into several classes (RS ≥ 2, extremely high; 1 ≤ RS < 2, high; 0 ≤ RS < 1, medium; −1 ≤ RS < 0, low; RS < −1, negligible) according to the average species sensitivity distribution curve and related literature7,57.

$$RP^k = \log \mathop {\sum}\limits_i {RQ_i^k}$$
(9)
$$\it \it {{{\mathrm{RS}}}} = \max \left\{ {RP^k} \right\}$$
(10)

Where the \(\mathop {\sum}\nolimits_i {RQ_i^k}\) represents the total RQ of pesticide mixture in compartment k.

Identification of regions of concern

Vulnerable regions were then identified by juxtaposing pesticide risk map with population, biodiversity, and water resource risk maps. The dataset of the population in China at 1 km resolution was obtained from the Resource and Environment Science and Data Center58. We used the zonal statistics tool of ArcGIS 10.2 to calculate the average population intensity in each basin. The river basins with >471 persons km−2 (80th percentile of population intensity) were classified as highly populated areas. Moreover, the biodiversity was quantified by the richness of tetrapods (including mammals, birds, reptiles and amphibians). The biodiversity data was obtained from the Chinese Academy of Sciences59. Therefore, the average richness of tetrapod per square kilometer was obtained and the 80th percentile of richness (781 species per km2) was identified as high biodiversity. The risk to water resources in China was estimated by water supply reported in the water-risk database from AQUEDUCT-v2.160; again the 80th percentile of the water risk score (3.5) was considered as high risk.

The level representing the number of highly concerning indicators of each basin ranged from 0 to 4. The four highly concerning indicators were as follows: extremely high pesticide ecological risk (RS > 2), high population intensity (>471 persons km−2), high biodiversity (>781 species km−2), high water risk (>3.5) (Supplementary Table 4).

Scenario setting

To assess the effect of public policy and management, several scenarios (Supplementary Table 15) were designed, to consider possible source reduction and substitution of target pesticides. According to the National Bureau of Statistics, the pesticide application has decreased continuously since 2015, when the MARA formulated the principles and targets for pesticide reduction. S1 was designed such that pesticide application would continue to decline by the average rate in the last years (about 8%). The utilization rate of pesticides was 40.6% in 2020, so to achieved the goal of MARA, this rate could be improved to 43.6% in 2025 (S2) and the amount of pesticide application was designed to be reduced by 6.9%. In addition, there has been a noticeable gap of pesticide application between China and developed countries with the pesticide utilization rate ranging from 50% to 60%. Therefore, S3 was set to decrease the pesticide application by 26% in 2025 compared to 2020 when its utilization rate reached 55% (approximately the level of developed countries). Meanwhile, MARA enacted a list of main low-toxic and low-residue pesticides used in crop production. Therefore, with a consideration of pesticide function and ecotoxicity, S4 was designed to consider that the target insecticides carbofuran and chlorpyrifos were replaced by diflubenzuron and chlorantraniliprole, and the herbicides atrazine and acetochlor were substituted with metolachlor and sethoxydim, respectively, from 2017 to 2020, with the same usages as those before substitution. By changing the input of the pesticide properties, new risk map was calculated with these modified models.