Introduction

The global population surge has led to a large increase in food demand, resulting in the large and widespread use of pesticides and fertilizers1. Agricultural diffuse pollutants such as nitrogen, phosphorus, pesticides, and antibiotics, lead to eutrophication, and damage to aquatic ecosystems2. Owing to climate warming, two-thirds of the world’s land is exposed to wetter and more variable precipitation patterns3, leading to increased agricultural diffuse pollution losses and challenging watershed ecological security4. A study showed that interannual variability in nutrient discharge loads within watersheds was up to 2.3 times higher, with 76% attributed to annual precipitation variability5. In extreme precipitation years, annual nitrogen loads are 30% higher than the long-term median6. However, the response of pesticide discharge levels to precipitation variability within watersheds remains unclear7, despite two-thirds of the global agricultural land being at risk of pesticide contamination.

The pesticide transport in watersheds are associated with rainfall events, pesticide application, soil type, and slope8,9. Among these factors, soil type and slope have been relatively stable parameters over a long period, and their changes have a negligible impact on the annual variance in pesticide loss. By contrast, precipitation as a meteorological factor is a major factor driving pesticide transport in watersheds10. Most studies on pesticide diffuse pollution acknowledge that pesticide discharge is driven by precipitation and runoff events; however, due to the random occurrence of precipitation events, trends are often discussed by simply comparing precipitation and pesticide losses in parallel, ignoring the multi-year trends in precipitation and the effects on pesticide discharge loads11. All calculated loads are estimates because the amount of data necessary to calculate the actual load is always insufficient. Current estimation methods include mechanistic models, empirical models, and field observations. Although long-term trends in pesticide emissions can be reliably simulated using mechanistic models, these models require local data homogeneity and complicated parameters, and the accuracy of the simulations in large watersheds (>500,000 km2) is questionable1. The limited scale of field observations remains a challenge for monitoring high spatiotemporal frequencies over large regions12. Current research is deficient in pesticide discharge in large watersheds. Empirical models are commonly used for estimating the agricultural diffuse pollutant losses in large watersheds such as the modified export coefficient model. The modified export coefficient model reasonably estimates the spatiotemporal variances of diffuse pollutants discharge loads by employing physical processes, such as precipitation, runoff, slope, canopy interception, and anthropogenic export parameters13. Although the modified export coefficient model has been widely used in estimating diffuse pollutants, it is still not involved in estimating pesticides loads in large watersheds.

Pesticide application is a dominant factor associated with pesticide discharge loads in watersheds14. Pesticide runoff loss typically accounts for 5.9% (0.56–14%) of the applied amount11. The primary factor in assessing watershed vulnerability is pesticide application intensity15; however, temporal trends between pesticide usage and discharge loads remain unclear. Pesticide application is mainly regulated by the demand even if it is always excessive, and the trends fluctuate; thus, the temporal impact on pesticide discharge loads has largely been overlooked. Because of the increasing concern about diffuse pesticide pollution, pesticide reduction policy under policy intervention are mandatory to reduce the amount of pesticide use. A zero-growth pesticide usage policy implemented by government regulators in China in 2015 led to a 152,000 tons reduction in usage and a 1.6% increase in efficiency in 1 year16. From 2016 to 2019, a further 390,000 tons reduction was observed. Correspondingly, field monitoring campaigns in China have demonstrated decreased concentrations of pesticides in some rivers17,18. Although the pesticide reduction policy mitigates diffuse pollution, ignoring the long-term precipitation effects may overestimate the actual efforts to discharge loads, which should be addressed by conducting long-term exploration to fill in the gaps.

Under the dual impact of increased precipitation variability and pesticide reduction policy, identifying the driving factors of pesticide discharge loads in the watershed is essential. A representative herbicide, atrazine, characterized by high usage, high mobility, high persistence, high detection, and endocrine disruption, was used as the target pesticide19,20. Atrazine is a widely used herbicide in corn, broomcorn, sugarcane and other crops. The annual usage of atrazine in USA and China has been ranked number two among conventional pesticides21. Atrazine is moderately dissolved in water and is highly persistent in the environment, with half-life up to 2 years22. Therefore, atrazine is ubiquitously in aquatic environments with 100% detected frequency, and poses potential risks to aquatic life23. Based on the model simulation, field observations, and sensitive analysis, we analyzed the spatiotemporal pattern of atrazine discharge loads and the change of the driving factors of atrazine discharge loads before and after the pesticide reduction policy in the watershed.

Results

Decadal precipitation and herbicide usage variation in a watershed scale

Precipitation in the YRW has increased significantly (p < 0.01) from 2010 to 2019 (Fig. 1B). The average annual precipitation has increased from 453 ± 50 mm during 1990–2009 to 585 ± 62 mm after 2010. After 2012, the average annual precipitation in the YRW generally exceeded 600 mm, with the maximum annual precipitation of 641 mm occurring in 2017. The precipitation increased significantly (p < 0.05) in the watershed in 2017. In addition, annual precipitation showed wide swings between the wet and dry extremes. The variance in annual precipitation between 2010 and 2019 was higher than that between 1990 and 2009. From 2011 to 2012, precipitation increased by 137 mm, whereas from 2014 to 2015, precipitation decreased by 59 mm.

Fig. 1: Background information and database of the Yellow River Watershed.
figure 1

A Land use, sub-basin delineation, sampling sites location, digital elevation model (DEM), normalized difference vegetation index (NDVI), and slope data of Yellow River Watershed (YRW). B Annual precipitation of YRW from 1990 to 2019 (mean ± standard deviation). C Arable area, herbicide usage, and herbicide intensity in YRW from 1990 to 2019.

Full size image

On a decadal scale, herbicide application intensity in the YRW showed an increasing trend (Fig. 1C). The arable area in the YRW remained at 13–14 million ha, and herbicide usage increased by 4.9-fold between 1990 and 2015, resulting in a 4.6-fold increase in herbicide application intensity. Since 2015, the government has implemented a zero-growth policy for herbicides. Both herbicide usage and application intensity decreased by 24% between 2015 and 2019, with an average annual decrease rate of 3.61% and 3.69%, respectively. The herbicide usage and application intensity in 2019 was 29,900 tons and 2.14 kg/ha, respectively, equivalent to that in 2011. Reducing pesticide usage effectively mitigated the discharge load exports to the aquatic environment.

Comprehensive water inflow coefficients of diffused pollutants

Among the five factors considering the comprehensive water inflow coefficients, the terrain and the retention index were stable on a decadal scale. While the precipitation, runoff transport, and leaching indices varied with the annual precipitation. Therefore, the spatiotemporal patterns of the diffuse pollutant discharge into watersheds were calculated using precipitation as the independent variable and the comprehensive water inflow coefficient as the dependent variable. The annual spatial distributions of the five factors are shown in Supplementary Fig. 1.

The comprehensive water inflow coefficients level of diffuse pollutants showed higher level in intensive agricultural areas than others within the watershed (Fig. 2A). The Hetao and Fenwei Plains, characterized by gravity-fed irrigation and traditional irrigation, respectively, are critical discharge areas for diffuse pollutants because of their intensive agricultural activities. In the gravity-fed irrigation areas, arable land was located on both sides of the mainstream for the convenience of irrigation (Fig. 2B). Diffuse pollutants were easily discharged from the soil into the water with a comprehensive water inflow coefficient close to 1 when the mainstream served as the irrigation network. The comprehensive water inflow coefficients of diffuse pollutants in the traditional irrigation areas were slightly low because tributary confluence lengthened the transport path to the water, and arable land was mainly located beside the tributaries.

Fig. 2: Spatiotemporal distribution of the comprehensive water inflow coefficients of diffuse pollutants in the Yellow River Watershed.
figure 2

A Spatial distribution of the comprehensive water inflow coefficients of diffuse pollutants in the Yellow River Watershed. B Average level (mean ± standard deviation) of comprehensive water inflow coefficients of diffuse pollutants in Hetao Plain (HT, covering sub-basin 22), Fenwei Plain (FW, covering sub-basin 17, 18, 39, 41, 43), and other areas in the Yellow River Watershed. C Temporal pattern of the leaching index (LI), precipitation index (α), transportation index (TI), and comprehensive water inflow coefficients (λ) in the Yellow River Watershed. D The coefficient of variance (CV) of LI, α, TI, and λ during 2010–2014 (before) and 2015–2019 (after).

Full size image

The comprehensive water inflow coefficient followed the variance of precipitation but showed more variability than precipitation did (Fig. 2C, D). The coefficient of variation among these three indices suggests that the variability of the precipitation and runoff transport index was higher than that of the leaching index, showing a conservative relationship with precipitation variance. Prior to 2015, the comprehensive water inflow coefficient varied conservatively with precipitation, which mainly reflected in 2010–2012. The variability of the coefficient was temporally enhanced. Notably, the precipitation increased by 70 mm (12.3%) in 2017 compared to that in 2016, corresponding to a 1.29-fold increase in the comprehensive water inflow coefficient. Under the impact of climate change, precipitation increased the variability of diffuse pollution discharge to watersheds.

Spatiotemporal distribution of atrazine load and field evidence

Most studies on atrazine transport loss in watersheds have explored Corn Belts in the United States because of the high corn production and high atrazine application (Table 1). Compared with agricultural watersheds with similar soil properties, land use proportion, and slope, the ratio between atrazine and total nitrogen in the Walnut Creek watershed shared the same alfisol and cambisol soil types and average slopes with the YRW. Although Walnut Creek had a large proportion of cropland, various parameters in the known watersheds were most similar to those in the YRW. In addition, research on the relationship between atrazine and total nitrogen was abundant (N = 34). The ratio of total nitrogen and atrazine loads from the Walnut Creek watershed was used as the θ parameter in this study. The average annual discharge load of atrazine in the YRW was estimated using the modified export coefficient model and was found to be 101.7 ± 3.2 tons, accounting for 18% ± 2.4% of the total herbicide use. The spatial pattern of the runoff discharge loads of atrazine (Fig. 3A) matched the comprehensive water inflow coefficients, which mainly highlighted two intensive agricultural activity areas, gravity-fed irrigation and traditional irrigation. The runoff discharge loads of atrazine showed the similar pattern in different years 2010, 2015, and 2019 (Supplementary Fig. 2). In the gravity-fed irrigation area Hetao Plain, the annual runoff loads of atrazine were 14.7 tons, accounting 2.5% of the atrazine applied amount, whereas in the traditional-irrigated Fenwei Plain, annual runoff loads of atrazine were 9.7 tons, accounting 1.4% of the atrazine applied amount. The field observations also highlighted these intensive agricultural activity areas (Fig. 3B), the mainstream was dispersed into several branches between sites BYGL and SHHK, and the atrazine concentration in SHHK increased 1.4-fold than that in BYGL, indicating a high contribution of atrazine in the gravity-fed irrigation areas. The traditional irrigation areas with intensive arable land dedicate a lot of atrazine discharge loads, such as TR, FR, and WR, whose controlled sub-basins showed high atrazine discharge loads.

Table 1 The background information of soil type, land use, and slope among watersheds.
Full size table
Fig. 3: Spatiotemporal distribution of annual runoff discharge loads of atrazine and the field observation.
figure 3

A Annual runoff discharge loads of atrazine simulated by the modified export coefficient model. B Atrazine concentration in surface water from the field observation. C Temporal variation among inflow coefficients of atrazine, runoff discharge loads of atrazine, precipitation, and herbicide usage smoothed by spline model.

Full size image

The temporal discharge variations and comprehensive water inflow coefficients were typically consistent (Fig. 3C). After the launch of the pesticide reduction policy in 2015, atrazine usage decreased, reducing the runoff discharge loads of atrazine, except in 2017. The impact of the comprehensive water inflow coefficient on pesticide reduction increased the runoff discharge load of atrazine. Although pesticide usage decreased by 10% from 2016 to 2017, precipitation increased by 12.4%, increasing the atrazine water flow coefficient and discharge load (31.0% and 19.1%, respectively). In other years, fluctuations in precipitation did not affect the decreasing trend in atrazine discharge loads with reduced pesticide use.

Impact of increasing precipitation on atrazine discharge load

In the case of stable arable land areas, the application intensity of atrazine was mainly subject to variations in atrazine usage, affecting the comprehensive water inflow coefficients. A Bayesian Gaussian process regression model involving the full Bayesian Monte Carlo sensitivity analysis method was used to explore the sensitivity of the application intensity of the atrazine discharge loads and comprehensive water inflow coefficients by dividing the periods before and after the pesticide reduction policy into two groups.

The main response factors for atrazine discharge loads changed substantially after implementing the pesticide reduction policy (Fig. 4A). Whether before or after the pesticide reduction policy, the response level of comprehensive water inflow coefficients to atrazine discharge loads was stable ranging from −0.3 to 0.2. However, the response level of atrazine application intensity significantly increased after the pesticide reduction policy. The response level for atrazine application intensity was only between −0.06 and −0.02 before, but it increased to −4.75 to 4.44 after implementing the pesticide reduction policy. In addition, the overall sensitivity indices (Fig. 4B) indicated the interactions between atrazine application intensity and the comprehensive water inflow coefficients to discharge loads of atrazine. The atrazine discharge loads were controlled by the interaction of the comprehensive water inflow coefficients and atrazine application intensity due to the total sensitivity index was close before the policy. The comprehensive water inflow coefficients indicated a stronger influence to atrazine discharge loads. After implementing the policy, the discharge loads of atrazine was mainly influenced by atrazine application intensity while the effect of the comprehensive water inflow coefficients sharply decreased, indicating that the discharge loads of atrazine were mainly controlled by the application intensity after the policy intervention.

Fig. 4: Sensitivity analysis of two parameters to runoff discharge of atrazine.
figure 4

A Response of two parameters to runoff discharge simulation before reduction policy and after reduction policy. B Total effect sensitivity indices of export intensity of atrazine and the inflow coefficient (λ) to runoff discharge of atrazine before and after the reduction policy. The center line denotes the median value (50th percentile), while the box contains the 25th to 75th percentiles of dataset. The black whiskers mark the 5th and 95th percentiles. C Precipitation and atrazine discharge load in actual value and simulated values using spline model.

Full size image

However, in unusual years, when precipitation increased significantly, the discharge load of atrazine increased. To explore the impact of the unusually increased precipitation, we used data from other years to estimate the precipitation and discharge loads and then compared the model-simulated precipitation with the actual precipitation in 2017 (Fig. 4C). The simulated scenario showed a discharge load 56.8 tons in 2017 when the actual atrazine discharge load was 74.5 tons. The simulated precipitation was 12.4% less than the observed values, and the corresponding runoff discharge was 23.7% less than the observed values. This result suggested that the effect of precipitation replaced pesticide reduction as the driving factor in years with an unusual increase in precipitation, and in years with unusual increase in precipitation, it promoted discharge. This difference suggests that the unusual increase in precipitation drives the discharge loads of atrazine in the context of pesticide reduction policy.

Discussion

The atrazine concentration in surface water reflects the regional discharge level. Therefore, we collected 180 sets of atrazine concentration in surface water from the relevant literature worldwide (Supplementary Table 4). The global average level of atrazine in surface water was 430 ng/L. The atrazine concentrations in different countries varied by up to three orders of magnitude, with developed crop-farming countries, such as the United States, Australia and China, exerting high concentrations (Fig. 5A). In addition to crop farming, sugarcane is an atrazine-consuming source, and countries, such as Brazil and Myanmar, have harmful atrazine concentrations in surface water. In the subregions, the highest atrazine level, 552 ng/L, was found in the Asia-Pacific region, followed by North America and Latin America (Fig. 5B). Europe banned atrazine application in the 1990s, and its atrazine levels in surface waters have decreased19. Notably, the regions with the largest number of studies were typically the areas most polluted with atrazine, such as the Corn Belts in the United States. The United States is the largest corn producer globally with an annual production more than 350 million tons;24 the average concentration of atrazine in the surface water in the United States was 794 ng/L. Although Australia’s corn production was not high, the corn is intensively cultivated on the Great Barrier Reef, resulting in high atrazine concentrations in surface waters25. China is the second largest corn producer globally and has huge atrazine usage20. The average concentration of atrazine in the surface waters of China (535 ng/L) exceeds the global average, while the investigations about atrazine transport behavior and discharge loads in China are inadequate.

Fig. 5: Occurrence of atrazine concentration worldwide and the temporal distribution of global herbicide usage.
figure 5

A Distribution of atrazine concentration worldwide from 2000 to 2021. B Concentration level in different continents/regions. The center line denotes the median value (50th percentile), while the box contains the 25th to 75th percentiles of dataset. The whiskers mark the 5th and 95th percentiles. C Global herbicide usage in the world, China, and Yellow River Watershed (YRW) from 1990 to 2020. D Growth rate of herbicide usage in the world, China, and YRW during the pesticide reduction policy.

Full size image

Atrazine exposure in surface waters is closely related to the use of herbicides. Global herbicide usage reached 1 million tons in 2004, increased by 1.2-fold from 1990 to 201226, and leveled off between 1.37 and 1.43 million tons (Fig. 5C). Notably, agriculture in China is a crucial part of the global agricultural industry. Herbicide use in China increased to 141,244 tons in 2014, an increase of 135% from 1990. In 2015, to mitigate the risk of diffuse pollution, the Chinese government launched a zero-growth strategy for pesticides16. This policy led to a staged decline in pesticide use, especially the 10% decrease in 2017–2018 (Fig. 5D). The YRW in China also attained this goal: with a 13.6% decline in herbicide use in 2018–2019. The application of the zero-growth policy effectively reduced herbicide usage, mitigating atrazine levels in surface water. Therefore, the dataset for USA and China in the literature was divided into two groups for comparison: before (n = 21) and after (n = 17) implementation of the zero-growth policy (Supplementary Fig. 3). The average atrazine level before the implementation of the policy in China was 259 ng/L, which decreased to 105 ng/L thereafter. In the YRW from 2014 to 2019, the average concentration decreased from 154.8 to 99.1 ng/L17. However, in USA, another atrazine intensively used country, the average atrazine level slightly increased from 1250 to 1321 ng/L in this period.

Spatiotemporal variations and driving factors of atrazine were abundant in small agricultural watersheds. In small watersheds, atrazine discharge levels are affected by the atrazine application intensity, precipitation, soil type, and slope, the main factors determining the ratio of atrazine to total nitrogen loss9. Atrazine concentrations were proved to reach the highest during the first few runoff events after application and decreased rapidly thereafter, with a 78% loss within 1 month of application27. However, atrazine showed excellent persistence and can be detected in the aqueous environment 20 years after concentration decay19. Regarding spatial distribution, herbicide concentrations in tributaries with intensive agricultural activities are usually higher than those in the mainstreams28,29. The average tributary concentration in this study was 9.1 ng/L, which was significantly higher than the average concentration of 5.8 ng/L in the mainstream. However, when the scope expands to the large watershed, the discharge loads of atrazine are mainly obscured by high spatiotemporal heterogeneity30. The pesticide discharge areas are mainly located in intensive agricultural areas. In contrast with the tributary confluence in traditional irrigation areas, dispersing into an irrigation network from the mainstream in gravity-fed irrigation areas shortens the transport pathway. The dispersed river channel form evolves into a vast irrigated farmland area; thus, atrazine can be directly discharged to the mainstream31. There was a 1.4-fold increase in atrazine concentration as the Yellow River water flowed through the gravity-fed irrigation area. This finding suggests that agricultural structures with short migration distance have a great risk of pesticide diffuse pollution.

Chemical loads, products of water discharge, and chemical concentrations are good indicators of chemical loss over time32. The modified export coefficient method provides reliable discharge loads for atrazine in watersheds13. The spatial discharge loads of atrazine were generally consistent with the measured concentration patterns (Fig. 3A, B). The two major parameters of the export coefficient method, λ and E, reflect the dynamics of precipitation and herbicide application intensity, respectively. In spatial distribution, atrazine discharge loads were closely related with comprehensive inflow water coefficients. Although annual atrazine use was lower in gravity-fed irrigation areas than in traditional irrigated areas, the short migration distances resulted in a high comprehensive water inflow coefficients and discharge loads in gravity-fed irrigation areas, which contributed the largest proportion of atrazine loads (14.4%) within the watershed. Temporally, the climate change drives the precipitation variability and further impacts the comprehensive water inflow coefficient. The comprehensive water inflow coefficients and discharge loads of atrazine was greater in wet years with higher comprehensive water inflow coefficient than that in dry years33. Typically, because the atrazine comprehensive inflow coefficient is influenced by both precipitation and runoff, cumulative effect of extreme runoff events and extreme precipitation in the watershed leads to higher variability to the comprehensive inflow coefficient than precipitation34. Increased precipitation variability enhances the frequency of runoff events during wet years and farmland soil erosion, promoting atrazine discharge loads in large vulnerable watersheds. In addition, increasingly wet environments in a large watershed and on a global scale are accompanied by increasing precipitation variability3. Because the variability of comprehensive inflow water coefficients combined the interaction with precipitation, runoff and leaching process, the increasing of precipitation variability can enlarge the diffuse effects of pollutants within watersheds35. During the period of increasing precipitation variability, a 12.4% increase in precipitation from 2016 to 2017 in this study led to a 31.0% and 19.1% in atrazine water flow coefficient and discharge load, respectively, despite the pesticide reduction policy was already implemented (Fig. 3C). Therefore, instances of drought and flood increase profoundly affect pesticide migration behavior and environmental risks at the vulnerable watershed and deteriorate the process of the pesticide reduction policy.

Some anthropogenic efforts have mitigated diffuse pollution through regulatory management18. Long-term water quality monitoring of the Yellow River showed that agricultural chemical reduction made the largest contribution to improved water quality when considering precipitation, runoff, evaporation, reservoir storage, and anthropogenic discharge. Mechanistic evidence through Fourier spectra and wavelet analysis shows that anthropogenic activities exceeded most natural biogeochemical processes at all spatial scales of the watershed. Systematic changes in the spectral properties of water flow have been associated with anthropogenic activities36. Since the government launched a pesticide reduction policy in 2015, increasing pesticide efficiency and reducing pesticide usage have mitigated the intensity of diffuse pollution in the YRW region. Before the pesticide reduction policy, precipitation variation mainly controlled the interannual discharge load of atrazine within the watershed, while the effect of atrazine application intensity turned to major factor after implementing the policy (Fig. 4). Projections based on the current regulations suggest that chemical pesticides, including atrazine, will be completely replaced by biopesticides by 205014, and the diffuse pollution from chemical pesticides will be gradually eliminated. Although unusual increases in precipitation affect human efforts to optimize environmental pollution issues, we are confident in the on-going environmental optimization policies.

Notably, this study has several limitations. First of all, the spatial heterogeneity of atrazine application data affects the accuracy of the estimation37. The mismatch in spatial resolution between the atrazine application database and the watershed’s inflow coefficient increased the uncertainty of the estimation results. We acknowledge that the estimated atrazine discharge loads values remained significant uncertainty; however, this uncertainty did not affect the main points in this study, as the spatial resolution’s uncertainty did not propagate to the temporal trend variations and field observation verified the overall spatial pattern of atrazine. The development of a regular monitoring system is required to provide the basic data on the migration behavior of pollutants in the watershed. In addition, the pesticide reduction policy was launched for 8 years and the available data could only trance for 5 year, which needs continue to follow up the policy impact. Although this study is a regional simulation, the changes in precipitation caused by climate change are a global issue that affects us all3. Reducing the use of conventional chemical pesticides is also an inevitable path for sustainable agriculture14,38. Exploring the process of environmental changes under the dual influence of climate change and human policies is beneficial for adjusting current policies and contemplating future development directions.

The current status determines that pesticide reduction policy requires ongoing effort. Our comparison of the atrazine concentrations calculated from the annual atrazine loads and runoff discharge in the watershed demonstrated that the average atrazine concentration from 2014 to 2017 exceeded the surface water atrazine concentration limit of 3 μg/L (Supplementary Fig. 4). Excess atrazine concentration can damage the ecological environment, especially for algae, fish, and frogs, and has potential neuroendocrine effects on mammals20. Because of the spatial variability in the entire watershed and the high exposure to atrazine in agricultural areas29, we speculated that atrazine levels would exceed the surface water concentration limits in intensively cultivated areas during the first few runoff events after its application39. Such high levels of atrazine after runoff mainly occur in corn cultivation areas, such as the US Corn Belt35, the Great Barrier Reef region of Australia25. These regions in intensively atrazine applied areas should pay more attention to the ecological risks to aquatic atrazine under the climate change effect, and require continuously tactical efforts to reduce the atrazine usage and seek alternatives to green pesticides in the future. Current estimations of the atrazine discharge load are useful for understanding the dynamic patterns of pesticide diffuse pollution in the dual effects of climate change and policy regulations and for providing knowledge for future policy development. The scale gap in refined modeling of diffuse pollutants in watersheds needs to be bridged in the future to better support government decision-making.

Methods

Study area

The Yellow River is the second-longest river in China and extends 5,464 km from its source to the estuary, covering 795,000 km2 within the watershed (Fig. 1A). The Yellow River is one of the largest vulnerable watersheds worldwide, and it surrounds the Loess Plateau, which has the most severe soil erosion rate worldwide40. Tremendous amounts of diffuse pollutants are discharged with soil erosion in watersheds, threatening aquatic environmental quality and ecological health41. Agriculture is one of the most important economic parts in the Yellow River Watershed. The main crop types include corn and wheat. Currently, the corn acreage in the watershed is 15.6 million hectares, and corn production is 98.9 million tons24. The main location of intensive agricultural areas is on both sides of the mainstream. The Hetao Plain is the largest gravity-fed irrigation area in Asia, with the mainstream of the Yellow River serving as an irrigation channel. The Fenwei Plain is a traditional irrigation area, the fourth largest plain in China, and the largest alluvial plain in the middle reaches of the Yellow River. The map of Yellow River Watershed and its sub-basins were delineated using ArcGIS 10.2 software based on a 90 m resolution digital elevation model (DEM). All of the visualization results about the Yellow River Watershed in this study were based on this map.

Export coefficient model

The overall frame work of the research process and the methodology is provided (Fig. 6). The atrazine discharge load in the watershed was calculated for the sub-basins using a modified export coefficient model involving the comprehensive water inflow coefficient. Based on the conventional export coefficient method, the modified model incorporates water inflow coefficients and calculates the atrazine loads in each sub-basin13. The modified equation is as follows:

$$L=\theta \lambda \sum {E}_{i}{A}_{i}$$
(1)

where L represents the diffuse pollution load of atrazine (kg/a), θ represents the ratio between atrazine and nitrogen loads, λ represents the comprehensive water inflow coefficient, E represents the application intensity of atrazine for the i-th sub-basin (kg/(km2·a)), and parameter Ai represents the tillage area in the i-th sub-basin (km2). Detailed calculations of λ are provided in the Supplementary Material (Supplementary Methods).

Fig. 6
figure 6

Overall framework of the research process and the methods.

Full size image

The λ coefficient is determined by precipitation, slope, surface/underground runoff, and vegetation interception factors in the watershed. Among these factors, the terrain index is related to slope and is stable over long time in a large watershed. The retention index is related to the fractional vegetation cover, which were stable in the last two decades in over 95% of the Yellow River Watershed42. However, the precipitation, transport, and leaching index were driven by annual precipitation in the decadal scale (2.09 ± 0.53%)43. We also compared the decadal variance of precipitation in the study area from 1990 to 2019 (Supplementary Fig. 5). The average precipitation from 2010 to 2019 significantly increased compared to the previous 20 years. Based on this background, the precipitation was determined as the independent variable and the λ coefficient was determined as the dependent variable in this study.

Field sampling and quantification of the atrazine concentration

Thirteen sites were sampled in September 2019 in the mean stream and tributaries (Fig. 1A). Sampling site locations are listed in Supplementary Table 1. Sampling sites in main stream were determined with the entry and exit of Hetao Plain (BYGL and SHHK), Fenwei Plain (LM and TG) and other agricultural areas along the upstream, middle-stream and down-stream, which could represent the pesticides discharge level in these areas. In addition, four agricultural tributaries were sampled in sites joining the main stream: the Fen, Wei, Tao, and Huangshui Rivers, which could represent pesticide discharge levels in these tributaries. A 2 L sample of river water from each site was poured into a dark-glass bottle and then filtered with a glass fiber membrane (0.77 μm, Whatman, USA). The samples were stored at 4 °C prior to the next procedure.

Extraction and quantification of atrazine were performed as described by Ouyang et al.12. In sequence, a 6 mL mixture of dichloromethane-methanol (v/v, 80/20, acidized by 5 mM HCl), 2 mL methanol, and 5 mL ultrapure water were permeated via graphitized carbon black solid-phase extraction (SPE) columns to activate the SPE columns. Desethylatrazine-13C (30 ng) was used as the surrogate mixture in filtered water. Next, the water samples were pre-concentrated using SPE cartridges at a flow of 1.0 mL/min in vacuum. Then, 1.5 mL methanol and 6 mL dichloromethane methanol (v/v, 80/20, acidized with 5 mM HCl) were used to elute the extract. When the eluents were concentrated by nitrogen flow to 200 μL, the residual volume was fixed at 1 mL using methanol. Next, the residual was filtered using a 0.22 μm polytetrafluoroethylene membrane for equipment detection. Atrazine was quantified using a UPLC-MS/MS machine loaded with a liquid chromatograph and a mass spectrometer (Supplementary Table 2). Atrazine-13C (30 ng) was used as the internal standard. A 0.1% formic acid water solution was used as eluent A, and a 0.1% formic acid acetonitrile solution was used as eluent B. Extracted ion chromatograms of atrazine, desethylatrazine-13C, and atrazine-13C are shown in Supplementary Fig. 6.

Quality control

Quality control methods, namely duplicate field samples, blank field samples, laboratory blanks, and surrogates, were used to ensure measurement accuracy. The calibration equations, method detection limits, and method quantification limits for atrazine and its surrogates are presented in Supplementary Table 3. A loss control method was used for the extraction process44. The average recovery ± standard deviation of surrogates desethylatrazine-13C in surface water was 87 ± 8.3% (N = 13). The average relative percentage difference of atrazine was 7 ± 13.6% (N = 13).

Statistical analysis

R version 4.0.3 was used for data analysis and visualization45. A Fully Bayesian Monte Carlo sensitivity analysis scheme was used to identify the most sensitive factor for atrazine discharge loads in the watershed based on the Bayesian Gaussian process regression model in R package tgp46. Random Latin hypercube samples were drawn at each Markov chain Monte Carlo iteration to estimate the main effects and the total sensitivity indices. Atrazine application intensity and precipitation were set as input variables and the atrazine discharge loads was set as output responses. In the setting of variables in R program, the size of random Latin hypercube samples was set as 1000, the span was set as 0.3, the 3-vector of Monte-Carlo parameters (B)urn in, (T)otal, and (E)very were set as 4000, 8000, and 10, respectively.

The ANOVA assumptions were tested before comparing the means among sites. The spline model is a reliable tool for describing nonlinear relationships between complex variables47. The discharge load of atrazine in the unusual precipitation years was estimated using a spline model.