Pesticide mixtures in soil: a global outlook

In modern agriculture, pesticides are used in combinations to protect crops. The co-existence of multiple pesticides in soil can threaten the soil biodiversity that maintains ecosystem services, hence further posing a long-term risk on food security. Here, we introduce an assessment of global soil contamination by the residue of pesticide mixtures in nine cropping systems using a fully mechanistic, spatially explicit, and time-resolved model at 0.5° × 0.5° spatial resolution (approximately 55 × 55 km at the equator) fed with georeferenced agricultural quantities, soil properties, and hydroclimatic variables. We found that 8.3 million km2 of treated land have more than one detectable pesticide, with pendimethalin, glyphosate, paraquat, chlorpyrifos, and chlorothalonil being the five most frequently detected. The highest pesticide mixture content was found in the ‘orchards and grapes’ cropping system (95th percentile at 7.3 mg kg soil−1). Globally, the pesticide mixture in the topsoil of approximately 1.88 million km2 exceeded 1 mg kg soil−1 for more than 180 d in a year. We estimate that 0.2 million tonnes of pesticides leach below the root zone each year globally, with glyphosate contributing the greatest fraction. The major hotspots of soil pesticide contamination are located in South America and Asia, mainly in Brazil, Argentina, Chile, China, Malaysia, and Japan. Our study shows that the accumulation of pesticide mixtures in soil is a global environmental issue that has to be explicitly accounted for in the sustainability assessment of agricultural production. We propose the use of mechanistic modelling as a tool to aid in designing pesticide management strategies and minimise residue contamination.


Introduction
Intensive agriculture uses a variety of pesticides in high amounts to protect crops and secure food production. As of the year 2017, approximately 4.1 million tonnes of pesticides were used globally, corresponding to an annual per capita consumption of 0.5 kg person −1 yr −1 and a market value of USD 71 billion [1]. These figures include more than 500 active ingredients distributed across multiple chemical and functional classes [2] and spanned over 20 000 commercial pesticide products with different formulations [3].
The intensive use of pesticides exerts substantial stress on soil health by affecting soil biota responsible for maintaining soil functions. Many laboratory studies reported a shift in soil microbial community and a decrease in microbial growth and enzymatic activities as a result of exposure to single pesticides [4][5][6][7][8]. Ammonia-oxidising bacteria and archaea are the most sensitive to pesticides, as studies consistently observed a significant decline in their abundance in pesticide-contaminated soil samples [6,9,10]. Karas et al [6] also found that pesticides decrease the abundance of sulphur-oxidising bacteria that convert reduced sulphur compounds to sulphate and cause adverse effects on enzymes responsible for mediating the phosphorus cycle. Both ammonia and sulphur oxidisers are key drivers of nutrient cycles that transform nitrogen and sulphur into forms available for microbial and plant uptake, playing a vital role in sustaining soil fertility. In addition to the negative impacts on soil microflora, pesticides affect several earthworm functions by causing disruption to their enzymatic activity, increasing mortality, reducing fertility, altering their feeding behaviour, and decreasing the overall community biomass [11].
More dreadfully, the simultaneous exposure to multiple pesticides can lead to synergistic adverse effects that can deviate from the additive toxicity of single pesticides; though, our current understanding of the response of soil biota to pesticide mixtures is still inadequate [11,12]. For example, Yang et al [13] observed a clear synergistic effect on earthworms exposed to a quaternary mixture of chlorpyrifos, fenobucarb, clothianidin, and acetochlor, whereas van Hoesel et al [14] found that the effects of insecticides and fungicides on earthworm activity were intensified with the simultaneous presence of herbicides. It thus appears that the inputs intended to protect crop production-the pesticides-are potentially causing threats to food security in the long-term.
Despite the negative impacts observed on soil biota, pesticide contamination in soil remains loosely regulated. A review by Jennings and Li (2014) [15] reported that only 174 jurisdictions within 54 United Nations member states promulgated some regulatory guidance values to specify the maximum allowable pesticide residues in soil. Less than 30% of those jurisdictions regulate for more than 100 pesticides, and the guideline values for the same pesticide can vary from 2 to 10 orders of magnitude [15]. In addition, studies also emphasised that the current environmental risk assessment framework used by the European Food Safety Authority (one of the most stringent regulatory bodies) to authorise pesticide commercialisation is insufficient to protect soil ecosystem services [7,16].
Owing to the lack of regulation, large-scale quantification and monitoring of pesticide residue in the soil are limited. Most studies quantified only a few pesticide substances and these are confined to national-or field-scales [17][18][19], with Silva et al (2019) [20] being one of the most comprehensive studies that covers 76 pesticide residues across 11 European Union member states. While the lack of regulation discourages investments in monitoring campaigns, a comprehensive overview of the current state of global pesticide soil residue is crucial to convey the urgency of adopting sustainable pesticide use, increase public awareness, and substantiate the need for more stringent regulation.
This study puts forth the first global estimate of pesticide residue in agricultural soil. We target to identify the hotspots for pesticide contamination in the topsoil (TS) and leaching below the root zone (BRZ). To this end, we estimate the soil residues of the 92 most used pesticides in nine cropping systems using a mechanistic dynamic model fed with georeferenced databases that include agricultural practice, soil properties, and hydroclimatic variables. The model is solved at 0.5 • × 0.5 • spatial resolution (approximately 55 × 55 km at the equator) along a four-layer vertical soil profile down to BRZ, and spanned over a 48 year time window. This study covers 11.85 million km 2 of croplands within 155 countries and excludes pastures. We benchmark our model outputs against previously published field measurements and georeferenced global datasets. We analyse the long-term annual average pesticide residue and leaching rate of cumulative and individual substances. Model benchmark is provided in the supplementary information (available online at stacks.iop.org/ERL/16/044051/mmedia) while the interpretation of our analyses is provided in section 3.

Pesticide application rates
We used PEST-CHEMGRIDS v1.0 [21] to obtain the global georeferenced annual pesticide application rates in year 2015 for nine cropping systems that include six dominant (i.e. alfalfa, corn, cotton, rice, soybean, and wheat) and three aggregated crops (i.e. vegetable and fruits, orchards and grapes, and other crops). Following the classification in the USGS Pesticide National Synthesis Project [22], the 'vegetable and fruits' cropping system includes legumes, roots and tubers, bush fruits, and herbaceous crops; the 'orchards and grapes' includes nuts, fruit trees, and vines; while the 'other crops' includes other cereals, oil crops, and fibre crops (refer to table 2 in Maggi et al [21] for a detailed classification).
PEST-CHEMGRIDS v1.0 provides the high and low annual application rates of 95 active substances that represent approximately 84% of the global pesticide mass used in 2015. In this study, we excluded calcium polysulfide, Bacillus amyloliquefaciens, and petroleum oil from our assessment due to insufficient data on their physicochemical properties. Hence, we assessed 92 substances using their median application rates.

Model description
We estimated the pesticide residue in soil using a mechanistic, spatially explicit and time-resolved model, solved within the general-purpose multiphase and multi-component bioreactive transport simulator (BRTSim v4.0e [23]). The model explicitly describes the water, gas, and heat flow along a one-dimensional variably-saturated soil column, the transport of aqueous and gaseous chemical species, pesticide volatilisation, adsorption, and degradation.
BRTSim solves for the non-isothermal continuity and conservation laws using hybrid explicit-implicit numerical techniques within finite volumes (refer to the User Manual and Technical Guide for detailed descriptions [24]). The water flow is modelled using the Richards equation [25] in conjunction with the relative permeability-water potential-saturation relationships of the Brooks-Corey model [26]. The transport of aqueous species is modelled by the Darcy's advection and the Fick's diffusion equations. We excluded the advection of gaseous compounds, but their diffusion is explicitly described using the Fick's law.
Pesticide volatilisation and adsorption are modelled using the mass action law [23], i.e.
where [X R ] and [X P ] are the reactant and product concentrations, respectively, while x R and x P are their corresponding stoichiometric numbers. K is the Henry's constant for gas volatilisation. For adsorption, K is expressed as K = K OC × SOC × ρ soil θ , where K OC is the soil organic carbon partition coefficient in (m 3 kg-C −1 ), SOC is the soil organic carbon content in (kg-C kg-soil −1 ), ρ soil is the soil bulk density in (kgsoil m −3 ), and θ is the soil water content in (m 3 m -3 ).
We assumed that the degradation rate R of pesticides occurs by first-order kinetics and that this is moderated by biological activity, soil moisture content, temperature, pH, and organic carbon content, i.e.
where [X] and DT 50 are the pesticide concentration and half-life, respectively. SOC ref is the organic carbon content of the soil at which DT 50 was estimated.
Here, we used SOC ref = 0.02 kg-C kg-soil −1 , which is the global median. The control of pH on pesticide degradation is described by with as the concentration of H + , and K pH,L = 1×10 −9 mol l −1 and K pH,U = 1 × 10 −5 mol l −1 representing the parameters of inhibition around the optimal value of pH 7. The function f B describes the biological response to changes in soil saturation S and temperature T, expressed as where T L = 288.15 K (i.e. 15 • C) and T U = 313.15 K (i.e. 40 • C) are the lower and upper response temperatures that give an optimal rate in the range T = 21 • C-32 • C for mesophiles, while S L = S U = 0.46 are scalar parameters estimated from experimental data in Wickland and Neff (2008) [27] that give an optimal biological response at about S = 0.5. The model was deployed on a three-dimensional grid resolved at 0.5 • × 0.5 • resolution horizontally and extended vertically over two atmospheric and four soil layers. The atmospheric layers allow for water ponding and the exchange of heat and gas between soil and the atmosphere. The first three soil layers are variably saturated and extend down to 1 m depth in the root zone (RZ) with a thickness of 30 cm, 30 cm, and 40 cm, respectively. The last layer has a variable thickness that depends on the depth of either the equilibrium water table or the bedrock and has a constant soil water saturation.
Pesticides were applied at the first soil layer with the annual applied mass distributed over the year following the crop calendar as in Maggi et al [28]. We considered a 40% loss of applied mass due to crop interception and wind drift after Trevisan et al [29]. Water and heat enter and leave the soil column through boundary fluxes that include rainfall, irrigation, crop evapotranspiration, and solar radiation (short and long waves). Incoming water fluxes and solar radiation were applied at the first soil layer, while the evapotranspiration was allocated over the soil profile according to the crop root distribution. The root distribution was calculated based on a negative exponential distribution function assuming that 99% of the root mass was above the maximum root depth.
We modelled 32 768 grid cells containing the nine cropping systems. The cells were selected based on the harvested area maps distributed along with PEST-CHEMGRIDS, which were originally produced by Monfreda et al [30]. Hence, we analysed in total 11.85 million km 2 of croplands.

Model input data
We fed the model with global georeferenced datasets described in table S1 of the supplementary information and we briefly list them below.
Soil physical properties, including soil texture (sand, silt, and clay fractions), bulk density, pH, and organic carbon content were obtained from the Soil-Grids2.0 [31], while the soil porosity was obtained from SoilGrids1.0 [32]. The soil permeability, the pore volume distribution index and air-entry suction of the Brooks-Corey model, soil heat capacity, and heat conductivity were sourced from Dai et al [33]. The equilibrium water table was obtained from Fan et al [34], soil thickness from Pelletier et al [35], and soil residual liquid saturation from Zhang et al [36]. The model was initialised with the long-term average soil saturation and land surface temperature obtained from the CPC Soil Moisture dataset by NOAA/OAR/ESRL PSD, Boulder, Colorado, USA [37] and Menne et al [38], respectively.
We sourced the crop calendars from Sacks et al [39], which were used to calculate the daily pesticide application rates and irrigation volumes following a similar approach as in Maggi et al [28]. The crop water security indicator in Thenkabail et al [40] was used to determine if the agricultural land in a grid cell is irrigated.
We used the time series (from 1970 to 2017) of rainfall, atmospheric temperature, longwave and shortwave solar radiation, and potential evapotranspiration from the Climatic Research Unit time-series datasets [41] as the hydroclimatic forcing. The crop evapotranspiration was determined by multiplying the potential evapotranspiration with crop coefficient k c obtained from Allen et al [42]. The maximum root depth of each crop was also obtained from Allen et al [42].
The physicochemical properties of the 92 pesticides, including the molecular mass, Henry's law constant, half-life, adsorption affinity to organic carbon, and LC50 of earthworms were obtained from the Pesticide Properties DataBase [43].

Data analyses
We considered the top 30 cm of the soil (i.e. the first soil layer) as the TS, the top 1 m as the RZ, and the last soil layer as the BRZ. All simulations were run for 48 years and the outputs of the last 5 years were used for analyses described below.
The instantaneous residue M of pesticide X in the TS of cropping system j at time t was calculated as where C X and m X are the concentration (mol l −1 ) and molecular mass of pesticide X, respectively, V w is the water volume, and M s is the soil mass.
The total pesticide residue in cropping system j is the cumulative of all the 92 substances, i.e. The mass of pesticides leached BRZ within a time window per unit area was calculated as Similarly, the per-unit-area total leaching in a crop- In our analyses, major watersheds were identified based on AQUEDUCT-v2.1 [44]. The predicted no effect concentration of each pesticide is estimated as the LC50 of earthworms divided by a safety factor of 1000 as in the European Commission Risk Assessment guideline [45].

Model benchmark
We benchmarked our assessments against key environmental variables and the pesticide residues observed in European agricultural soils. Specifically, we benchmarked the modelled soil wetness, temperature, and pH in the TS (i.e. the top 30 cm) against values obtained from the CPC Soil Moisture dataset [37], the NOAA/NCEI land surface temperature dataset [38], and the SoilGrids2.0 dataset [31]. We compared the modelled pesticide soil residues and the number of detectable substances in the TS against field observations reported in Silva et al [20]. Results of the model benchmarking are reported in text S1 of the supplementary information.

Sensitivity of variables and quality of estimation
We used the AMAP index developed in la Cecilia et al [46] to assess the model sensitivity. AMAP quantifies the effect of a variable on the probability of a target output to exceed a threshold, with higher values signifying greater influences. The variables assessed include four boundary fluxes (i.e. pesticide application rate, rainfall, potential evapotranspiration, and net solar radiation) and seven soil properties (i.e. porosity, bulk density, permeability, pore volume distribution index, air-entry suction, pH, and organic carbon content). The target model outputs are the pesticide residue in the TS and leaching BRZ with corresponding thresholds of 0.01 mg kg-soil −1 and 1 mg m −2 yr −1 , respectively.
To measure the quality of our estimates (QI), we integrated the data quality QI APR of PEST-CHEMGRIDS v1.0 and the quality of benchmarking, i.e. QI = ( 3 · QI APR + QI S + QI T + QI pH ) /6. The overall quality of the application rate data QI APR was calculated as the average over all crops and all substances. The qualities of soil saturation (QI S ), temperature (QI T ), and pH (QI pH ) were calculated as the normalised absolute errors in soil saturation, temperature, and pH against global datasets described above. All quality indices are bounded between 0 and 1 such that QI = 0 and QI = 1 signify low and high data quality, respectively.

Pesticide residue in topsoil
Globally, the TS in croplands has an average total pesticide content of 0.62 mg kg-soil −1 and a 95th percentile of 2.12 mg kg-soil −1 , with residue heterogeneously distributed across geographic regions (figure 1(a) [47]). Specifically, the five countries with the highest average pesticide residue are New Zealand (30.7 mg kg-soil −1 ), United Kingdom (14.4 mg kgsoil −1 ), Cyprus (5.2 mg kg-soil −1 ), The Netherlands (3.6 mg kg-soil −1 ), and Japan (2.7 mg kg-soil −1 , table S2, supplementary information). The majority of croplands in Northern, Southern, and Western Europe have pesticide residues higher than the global average (i.e. about 76.8%), especially in Italy, Spain, Germany, and United Kingdom, while residues in Eastern Europe are relatively lower ( figure 1(a)). We estimated that pesticide residue in about 1.88 million km 2 of croplands exceeded 1 mg kg-soil −1 for more than 180 d in a year ( figure 1(c)). These regions are mainly located in China (0.83 million km 2 ), Brazil a. b.
c. (0.23 million km 2 ), and the United States (0.17 million km 2, table S2, supplementary information). Our analysis shows that more than 90% of pesticide residue remains in the adsorbed phase, especially in regions with high residue contents including China, Europe, and North America (figure S6, supplementary information), suggesting that sorption kinetics are important mechanisms that determined the pesticide accumulation in soil. Croplands in Africa, South and Southeast Asia, and Australia have relatively low pesticide residue, with an exception in a few countries such as Malaysia and the Philippines ( figure 1(a)).
We counted the number of pesticides with a residue greater than the typical laboratory measuring limit (i.e. 0.01 mg kg-soil −1 ) and found that 70% of global croplands (i.e. about 8.31 million km 2 ) have multiple pesticide residues in the TS ( figure 1(b) [47]). More than ten pesticide residues are detected in 16 major watersheds (figure S7, supplementary information), including Amur and Huang He (China), Colorado and Negro (Argentina), Rajang (Malaysia), Amazonas and Esmeraldas (Ecuador), Nile (Egypt), Ishikari (Japan), and Clutha (New Zealand). Croplands in Southern and Western Europe have, on average, four to five detectable pesticide residues despite the high total pesticide contents found.
Among the nine cropping systems, 'orchards and grapes' has the highest average total residue (i.e. 2.35 mg kg-soil −1 ) with the highest number of detectable pesticides (figures S8(a) and (b), supplementary information). Specifically, 'orchards and grapes' in Japan and Chile have more than 14 pesticide residues that summed up to more than 26 mg kgsoil −1 . 'Vegetable and fruits' also records high content of pesticide mixtures, especially in Malaysia and New Zealand, which have more than 24 mg kg-soil −1 of residue. Croplands for corn and soybean have comparable pesticide residue (i.e. a global average of about 1 mg kg-soil −1 ), while wheat fields have the lowest (i.e. a global average of 0.12 mg kg-soil −1 ).
The five most frequently detected pesticides in the TS are glyphosate, pendimethalin, chlorpyrifos, paraquat, and chlorothalonil (figure S7, supplementary information). Although their global average residue is relatively low (i.e. <0.4 mg kg-soil −1 ), some major watersheds are exposed to residue exceeding the predicted no-effect concentrations for earthworms ( figure 2(a)). In particular, watersheds in Argentina and Chile are exposed to chlorpyrifos residue of one order of magnitude higher than the noeffect threshold, while several watersheds in Asia and South America have average chlorothalonil residue exceeding the no-effect threshold by more than three times, such as in Lake Titicaca, Brahmaputra, Majes, Salween, and Rajang ( figure 2(a)). Among all pesticides assessed, the global 95th percentiles of 13 pesticide residues exceeded the no-effect threshold ( figure S9(a), supplementary information). Bromoxynil, phorate, and thiophanate-methyl are among the most concerning pesticides as their global average residues are more than two factors higher than the no-effect threshold ( figure S9(a), supplementary information).

Pesticides leached below the root zone
Every year, approximately 0.2 million tonnes of pesticides leach BRZ globally, corresponding to about 5.6% of the applied mass. The geographic b.
a. distribution of high leaching rate areas are substantially different from regions with high pesticide residue (figure 3(a) [47]). The leaching rate in equatorial regions is especially high, such as in Amazonas, Niger, and Congo. In particular, approximately 10 600 km 2 in the Amazonas watershed leach more than 100 mg m −2 yr −1 and, in total, 2800 tonnes leach BRZ of the watershed every year. In Asia, leaching is particularly high in parts of China, Japan, and the Philippines. Leaching greater than 200 mg m −2 yr −1 are observed in the southern regions of China, including Xi Jiang, Min Jiang, and Dong Jiang. A total of about 80 000 tonnes leach BRZs in China every year, contributing to 40% of the global pesticide leachate. Although European croplands see high pesticide residue in the TS, their leaching is generally low, with high rates mainly observed in the alpine regions between North Italy and Austria ( figure 3(a)). In figure 3(b) [47], we counted the number of pesticides that exceeded a leaching rate of 1 mg m −2 yr −1 , which corresponds to the regulatory limit of groundwater contamination stated by the European Commission (i.e. 0.1 µg l −1 [48]) after assuming that a shallow aquifer is located BRZ with a thickness of 23 m and a soil porosity of 0.43 (which are the global median values derived from SoilGrids1.0 [32], Fan et al [34], and Pelletier et al [35]). Leaching of more than two pesticides is mainly observed in the south-eastern of the United States, Brazil, Nigeria, Congo, China, Japan, Malaysia, and the Philippines. In particular, approximately 9% of croplands in China (i.e. about 0.13 million km 2 , mainly in Yangtze River, Xi Jiang, and Min Jiang) are subject to leaching of more than ten pesticides ( figure 3(b)).
Globally, 'vegetable and fruits' cropping system has the highest average leaching rate (i.e. 53 mg m −2 yr −1 ) and the greatest number of pesticides that exceeded 1 mg m −2 yr −1 , especially in Guatemala where more than ten pesticides exceeded the threshold (figures S8(c) and (d), supplementary information). On average, more than 10% of the applied mass in soybean and rice fields leaches BRZ. Wheat croplands have the lowest leaching rate.

Sensitivity of variables and quality of estimates
Pesticide application rate is the most influential variable for both the TS residue and leaching BRZ to exceed the defined thresholds of 0.01 mg kg-soil −1 and 1 mg m −2 yr −1 , respectively (figure S10, supplementary information). In general, model outputs are also sensitivity to variations in boundary fluxes. Relative to TS residues, potential evapotranspiration and net solar radiation have a higher control over the probability of exceeding the threshold than soil properties (figure S10(a), supplementary information), whereas the leaching rate is impacted by the rainfall amount and soil pH ( figure S10(b), supplementary information).
Since pesticide application rate has a major impact on the model outputs, we propagated the uncertainties in PEST-CHEMGRIDS toward quantifying the data quality of our model estimates as described in section 2.6. A high QI index signifies good quality of estimates. The distribution of the QI index is right-skewed with a median at 0.7 (figure S11(b), supplementary information). Less than 8% of the grid cells have an QI index smaller than 0.5 and are mainly located in Brazil, Saudi Arabia, Iran, and Western regions of China (figure S11(a), supplementary information).

Discussion
Our assessment extends the study in Silva et al [20] and addresses the soil contamination by pesticide mixtures at the global scale. We showed that the simultaneous existence of multiple pesticides is very much a normality in agricultural soil worldwide and, in some regions, the cumulative residues can exceed 10 mg kg-soil −1 .
More strikingly, residues of many pesticides (individually) exceeded the no-effect concentrations for earthworms, but we do not exactly know how synergistic effects of these mixtures affect earthworms. Among pesticides most frequently detected in the TS, chlorpyrifos, chlorothalonil, and dichloropropene have particularly high toxicities to soil biota and their residue exceeded the no-effect threshold in some areas. Previous studies show that the inhibitory effects of chlorpyrifos on soil microbial activity and population were significantly amplified in the presence of chlorothalonil [49,50], while noting that chlorpyrifos and chlorothalonil residues frequently coexist, especially in orchards and grapes cropping system. Yet, only less than 102 jurisdictions (not countries) worldwide have promulgated regulatory guidance values relative to soil residues for these compounds (i.e. 102 for chlorpyrifos, 75 for dichloropropene, and 67 for chlorothalonil [15]) with some guideline values being a few times higher than the LC50 of earthworms. For example, the guideline value of chlorpyrifos in Michigan, USA, is 11 000 mg kg-soil −1 [15], which is nearly two orders of magnitude greater than the LC50 of earthworms (i.e. 129 mg kg-soil −1 [43]). In addition to earthworms, pesticides can have adverse effects on soil microbial communities such as nitrifying bacteria and archaea that are responsible for transforming nitrogen into forms available for microbial and plant uptake [10]. Some pesticides can also increase the rate of ammonification [51], i.e. the fast conversion of amino acids into ammonium ions that are more prone to leaching, leading to nitrogen loss. More dreadfully, pesticides can reduce the symbiotic efficiency of nitrogen-fixing rhizobia bacteria [52], decreasing the ability for legume crops to replenish nitrogen in the soil. Hence, pesticide contamination may result in soil fertility loss. However, a comprehensive assessment of pesticide contamination risk on soil biota is currently unachievable owing to the lack of ecotoxicological data specific to each pesticide substance and soil biota species.
We also estimated that nearly 6% of the applied pesticides leach BRZ, where their degradation rates are substantially low. Pesticides may persist and accumulate in the deep soil and potentially contaminate the groundwater. Over the years, numerous pesticides were detected in the groundwater across continents, including the United States [53], Canada [54], Europe [55], and China [56], posing threats to the safety of drinking water supplies.
Our assessment suggests the need to establish a coherent regulatory guideline for pesticide contamination in soil that is in concordance with environmental and ecotoxicological studies to protect soil functions while acknowledging that there are many knowledge gaps yet to be filled. The combinatorial effects of multiple pesticides on soil biota, nutrient cycles, and crop uptake and metabolism are still not well understood, leading to high uncertainties on how pesticide residue affects soil fertility and crop yield. In a warming world, climate change adds on another layer of complexity in illustrating the role of pesticides in crop production, whether they are the protection or threats to the food security of the future.
Our model assessment provides a starting point to address the knowledge gaps. It identifies the pesticides, cropping systems, and geographic regions that require research investment and tailored strategies. Yet, we acknowledge the numerical and model limitations embedded in our estimates. We did not account for the transport and fate of the metabolites, which can retain the toxicity of parent substances and persist in the environment, and therefore we may underestimate the extent of contamination. Our modelling did not explicitly include surface runoff and soil erosion, the processes that can enhance contaminant dispersion. We conducted a parameter sensitivity screening by treating each grid cell on the map as an individual replica. This may not be the most ideal method for uncertainty quantification, but a comprehensive global sensitivity analysis is not feasible due to high computational cost, i.e. one full model run for the work presented here consumed more than 30 000 CPU hours.

Conclusions
We presented an estimate of global soil contamination by pesticide mixtures. Globally, the pesticide residue in the TS of about 1.88 million km 2 of croplands exceeded a concentration of 1 mg kg-soil −1 for more than 180 d in a year, while approximately 5.6% of the applied pesticides leach BRZ. High residue in the TS and leaching BRZ are estimated for 'orchards and grapes' and 'vegetable and fruits' cropping systems, whereas wheat fields have the lowest pesticide contamination. Croplands in parts of Europe, South America, China, and Japan have high pesticide residue with more than 10 pesticides detected in 16 major watersheds. Leaching is particularly high in equatorial and southern regions of China. Among the 92 pesticides assessed, glyphosate and chlorothalonil are the most frequently detected in the TS and to leach BRZ. Our assessment urges for more research efforts in understanding pesticide impacts on soil health and functions.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.6084/m9.figshare.12966323. Data will be available from 01 February 2021.