The locations of the 8 and 13 soil and groundwater sampling points, respectively, are shown in Figure 1 and Supplementary Table S1. This study focused on the agricultural area, so soil and groundwater in other land use in the eastern area were not densely collected. Moreover, the elevation of the groundwater wells conforms to the geographic profile. For groundwater sampling, the groundwater level was first measured using a water level meter; then, a bailer was gently dropped on the top of water column in the shallow groundwater wells until it was full. At this point, the water was transferred into an appropriate sample container. In addition, for deep groundwater wells, there was a pumping system installed for groundwater consumption and it was pumped out before collecting groundwater samples as it is a selected purging method that does not alter the geochemical and physical properties. Additional parameters measured on-site were pH, oxidation–reduction potential (ORP), dissolved oxygen (DO), electrical conductivity (EC), and temperature. The samples were stored at a temperature lower than 4°C in an ice box during transportation.

For the soil sampling, samples were collected in a sugarcane field in the study area based on different soil textures contained in the study area with different collecting methods, i.e., bulk soil sampling and soil core sampling. For the bulk soil sampling, each sampling point was collected at a depth of 15-cm from 5 different spots around the point. These five locations were approximately 10 m from each other. The sub-samples were then mixed together for representative soil at the sampling point. For the soil core sampling, surface soil to 15-cm depth was removed, and the core was collected using a 100-cm³ soil core sampler with duplicate samples. The core sampler was then hammered down to preserve the depositional sequence. Finally, bulk soil samples were air-dried for 1 week and then passed through a 2-mm sieve. Only soil particles ≤2 mm were retained for further adsorption experiments.

Soil core samples were used to determine the bulk density of soil at each sampling point after being oven-dried at 105°C for 3 days. The samples were also used to define hydraulic conductivity through the Falling head permeability test, indicating the potential for water flowing through the soil. Additionally, 2-mm diameter soil bulk samples were used to identify the texture of each sample, and pH and organic matter were measured from the bulk soil samples. 20 g of soil were added to 20 ml of distilled water (1:1 w/w) in the 60 ml PE bottle and were stirred regularly for 30 min. Then, the soil samples were left for 30 min until settled; pH of the water above the soil was determined by pH electrode (Pansu and Gautheyrou, 2006). Soil organic matter (%OM), in this study, was measured using the Walkley and Black Method (Mylavarapu, 2014), while soil texture was defined by Robinson’s pipette analysis (Augustin and Cihacek, 2016; de Oliveira Morais et al., 2019). Moreover, %OM can be converted into soil organic carbon (%OC) by the following equation: % OM = 1.72 % OC .

To define the adsorption coefficient of atrazine in soil (K_d), a batch adsorption experiment was carried out (Yue et al., 2017). Firstly, a 15-ml centrifuge tube was filled with 1 g of each 2-mm diameter bulk soil sample from the different points in the study area along with 10 ml of atrazine solution (standard atrazine with >97.0% purity purchased from Tokyo Chemical Industry Co., Ltd. in a background solution of acetonitrile and 0.01 mol/l CaCl₂ for maintaining the ionic strength). In this case, atrazine was added at an initial concentration of 0.5, 1, 5, 10, and 20 mg/l, respectively. Next, all tubes were sealed and shaken for 24 h; then, the suspension was centrifuged for 5 min at 5,000 rpm (Chefetz et al., 2004). Yue et al. (2017) claimed that the equilibrium time of atrazine adsorption was approximately 24 h as the adsorption occurred on the surface of soil organic matter. After the centrifugation, 2 ml of supernatant was filtered through a 0.45-µm pore size membrane and analyzed by High-Performance Liquid Chromatography (HPLC) with a model 490E UV detector, and a Hewlett Packard C18 column 250 mm, 5 µm particle size, 3 mm diameter. The mobile phase was a 60:40 mixture of acetonitrile and deionized water with a flow rate of 1 ml/min. Twenty µl of the sample was injected with a selected wavelength of 220 nm for each analysis with a retention time of 3.8 min; each adsorption experiment was performed in triplicate. Additionally, a blank sample (no soil) was prepared with the different initial concentrations.

The difference between the calculated initial atrazine concentration in the solution and the equilibrium concentration is the amount of absorbed atrazine in the soil. K_d was estimated as the ratio of the adsorbed atrazine concentration to the remaining concentration in solution at equilibrium. If the isotherm is linear, Kd corresponds to the isotherm slope. The defined K_d values were normalized to the fractional soil organic carbon content (f_oc) of each soil to determine the adsorption coefficient (K_oc) by the following equation: K oc = K d / f oc .

Linear, Freundlich, or Langmuir adsorption models have commonly been used to describe analytical adsorption isotherms (Schwarzenbach et al., 1993). To clarify the overall adsorption process of atrazine in all soil samples, Langmuir and Freundlich sorption models have been used. The Langmuir isotherm model assumes that the energy adsorbed on the surface of the adsorbent is uniform with no interaction between the adsorbed molecules (Langmuir, 1918; Sahu et al., 2016). On the other hand, the Freundlich isotherm model is used to describe heterogeneous surface equations such as the heterogeneity of the adsorbent surface, the adsorption energy, and the exponential distribution of the adsorption point (Freundlich, 1907). In this study, isotherms were modeled using Langmuir and Freundlich models and given their expressions since they show related results.

HYDRUS-1D is an environmental model for simulating fluid movement and contaminant transport in one-dimensional various saturated media (Šimunek et al., 2012). The Richards equation and Fickian-based advection-dispersion equations are used to govern solute transport in the liquid phase and diffusion in the gaseous phase (Ladu and Zhang, 2011). In this study, HYDRUS-1D was used to predict the hydraulic conductivity with the following equation, which is modified from Richard’s equation (Ali et al., 2021): ∂ θ ∂ t = ∂ ∂ x [ K ( ∂ h ∂ x + 1 ) ] , where h is the water pressure head (L), θ is the volumetric water content (L³L⁻³), t is leaching time (minutes), x is the spatial coordinate (positive equivalent to upwards), and K is the unsaturated hydraulic conductivity (LT⁻¹). It should be noted that Richard’s equation does not account for the effects of the vapor phase in water flow but considers only the liquid phase in the water mass balance.

The attenuation/retardation factor model (AF/RF model) is a tier-1 or simple model used to determine the leaching potential of a pesticide through the soil profile (Rao et al., 1985). This model can be used to solve large-scale problems because it does not require an excessive amount of data. The Attenuation Factor, or AF value, ranges from 0 to 1. The value indicates the possibility of pesticides leaching from the soil surface to the groundwater table. If the AF value is 1, it means that the pesticide is non-adsorbed and has a high groundwater pollution risk. In contrast, if the AF is 0, the pesticide is considered a strongly adsorbed pesticide and has a low groundwater pollution risk (De Paz and Rubio, 2006). The AF value can be defined by the following equation: AF = exp ( ( − ln 2 ⋅ d ⋅ RF ⋅ θ FC ) / ( q ⋅ t 1 / 2 ) ) , where d is the groundwater depth (m), θ_FC is the water content at field capacity, q is the water recharge through soil or water flow (m/d), and t_1/2 is the half-life of the pesticide (d). Additionally, RF is the retardation factor and can be computed by the following equation: RF = 1 + ( ( ρ b ⋅ f oc ⋅ K oc ) / θ FC ) , where ρ_b is the bulk density of soil (kg/m³), f_oc is the organic carbon content fraction, and K_oc is the adsorption coefficient (m³/kg). To evaluate the degree of leaching risk to groundwater of the pesticide under study, the AF value is classified into five classes, as shown in Table 1 (Khan and Liang, 1989). The required data were stored in a GIS, which was integrated with the AF/RF model to evaluate the potential of leaching groundwater of the pesticide.

Moreover, the Groundwater Ubiquity Score (GUS) index was used to assess the vulnerability of atrazine in the study area. This screening method uses the K_oc value and t_1/2 to rank the leaching potential of the pesticide to groundwater. A pesticide with a GUS value less than 1.8 is described as a “non-leacher,” a pesticide with a value greater than 2.8 is described as a “leacher,” and a pesticide with a value between 1.8 and 2.8 is considered “transitional.” The GUS value can be calculated by the following equation: GUS = ( log t 1 / 2 ) ( 4 − log ( K oc ) ) .

The chi-square test is the sum of the squared difference between the estimated value provided by the model and the experimental data (Ho and Ofomaja, 2006). The chi-square test was used to assess the adsorption isotherm with the following equation: X 2 =   ∑ ( ( q e − q e , model ) 2 /   q e , model )   , where q_e is the equilibrium data from the experiment (mg/g) and q_{e, model} is the amount of adsorbed atrazine per mass of soil (mg/g). If X² is a small number, it can be concluded that the observed data and calculated data are almost the same.

All 13 groundwater samples were filtered through a 0.45-µm pore size membrane to avoid clogging during injection, and stable isotopes in water (oxygen and deuterium isotopes: δ¹⁸O and δ²H) were analyzed using cavity ring-down spectroscopy (model L2130-i, Picarro Inc., Santa Clara, United States) at the Thailand Institute of Nuclear Technology (TINT). Moreover, the compositions of nitrogen isotopes (δ¹⁵N) and δ¹⁸O were analyzed by the denitrifier method (Sigman et al., 2001) at the UC Davis Stable Isotope Facility, University of California, United States. The analytical precision for δ¹⁵N and δ¹⁸O were approximately ±0.1‰ and ±0.3‰, respectively. All stable isotopic compositions (delta, δ) were reported in per mil (‰) notation using the following equation: δ ( ‰ ) = ( R sample − R standard ) R standard x 1000 , where R_sample is the isotope ratio of the sample and R_standard is the isotope ratio of the standard.

For analysis of cations (Ca²⁺, Mg²⁺, and Na⁺), all groundwater samples were acidified to pH < 2 with nitric acid (HNO₃) and loaded into an atomic adsorption spectrophotometer (AAS) (Analytik Jena model ZEEnit 700, Jena, Germany). For analysis of anions (F⁻, Cl⁻, NO₃ ⁻, NO₂ ⁻, SO₄ ²⁻, HCO₃ ⁻, and CO₃ ²⁻), the samples were loaded into an ion chromatography (IC) system (model Dionex ICS-2500, Sunnyvale, United States). Then, total alkalinity was determined from the samples to also measure HCO₃ ⁻ and CO₃ ²⁻ due to the limitation of AAS.

The results of the on-site parameter measurements indicated that the pH of the groundwater samples ranged from 7.10 to 8.40, with an average of 7.77, indicating normal to mildly alkaline conditions (Supplementary Table S2). The temperature of the groundwater samples was between 27 and 30.8°C. However, the temperature of some groundwater samples could not be measured because the samples were not collected directly from the wells. Dissolved oxygen (DO) varied in the range of 4.7–7.1 mg/l. The high amount of oxygen in the samples can be explained by the diffusion of oxygen from the atmosphere due to the shallow water level of the groundwater wells (Rose and Long, 1988). It is interesting that high DO was found although W11 was a deep well; this might imply that the groundwater in this well was oxygenated by the effect of groundwater pumping (Bonte et al., 2017). Electrical conductivity (EC) ranged from 356 to 4,850 µS, and the total dissolved solids (TDS) ranged from 268 to 2,430 mg/l. Fluctuation of the oxidation–reduction potential (ORP) was also found in this study, indicating that both oxidizing and reducing environments occur in the groundwater environments.

The result of the cation and anion groundwater analysis displayed that Cl⁻ and F⁻ were in the range of 1.00–2,105.70 mg/l and 1.97–20.58 mg/l, respectively. Interestingly, NO₃ ⁻ concentrations were relatively high in W3, W5, W6, and W10, ranging from 22.29 to 96.01 mg/l, as shown in Figure 2; however, NO₃ ⁻ concentration groundwater detected in all samples was lower than the drinking water standard (50 mg/l). High NO₃ ⁻ concentrations were mainly found in an agricultural area in the western part of Song Phi Nong District, which is one of the possible sources of NO₃ ⁻. Moreover, Ca²⁺, Mg²⁺, and Na⁺ of these groundwater samples ranged from 3.50 to 204.40 mg/l, 1.50–102.00 mg/l, and 4.50–477.00 mg/l, respectively. SO₄ ²⁻ and NO₂ ⁻ could not be detected in some of the samples due to their low concentrations.

In addition, the result of the cation and anion measurements is normally used to define the water type by the Piper diagram (Zhang et al., 2015), as indicated in Figure 3. This diagram graphically displays the relative proportions of major cations and anions to facilitate the identification of geochemical facies. Six types of water were found in this experiment, including Ca–Cl, Ca–Na–Cl, Ca–HCO₃–Cl, Na–Cl, Na–HCO₃–Cl, and Ca–Na–HCO₃, as shown in Supplementary Table S3. It can be concluded that different types of water arose from different sources, and some of the groundwater samples had Cl as a dominant constituent, indicating that groundwater in some areas was affected by human activities, especially fertilizer application and wastewater discharge (Ogrinc et al., 2019).

The isotopic ratios of the water samples indicated that the δ¹⁸O and δ²H of the samples ranged from −7.75‰ to −3.44‰ and from −53.66‰ to −25.63‰, respectively. The relationship between δ¹⁸O and δ²H is illustrated in Figure 4 as an evaporation trend (orange line) and has the linear equation of δ¹⁸O = 4.23776 δ²H—17.392. The slope of the line is different from one area to another depending on the local climate (Zhang et al., 2013). The relationship was compared to the Bangkok Local Meteoric Water Line (BKK LMWL) because Suphan Buri Province has almost the same elevation and distance from the sea as Bangkok. According to Figure 4, the evaporation trend was almost the same as the BKK LMWL, showing that the groundwater table was shallow and that most of the waters were contaminated by human activities on the ground surface, particularly fertilizer application. This occurred because the land had been used mainly for sugarcane fields for a long period of time. The result also indicated that surface water and groundwater are connected. In contrast, W12 was far from the BKK LMWL, indicating that the groundwater level was deeper than those of other wells. The ORP of W12 revealed a reducing condition, indicating that the water comes from a deep confined aquifer. As shown by the map, the NO₃ ⁻ in W12 may not be from fertilizer application due to the deep groundwater level. These results can be implied that groundwater resources in this study area can be contaminated by atrazine as the area is used for agriculture. Moreover, the δ¹⁵N_Air and δ¹⁸O_VSMOW in the study area varied from 3.51 to 28.42‰ and −3.92 to 19.79‰, respectively, as illustrated in Supplementary Table S4. The values of δ¹⁵N_Air and δ¹⁸O_VSMOW at W12 and W4 were mostly higher than those at the other sampling points.

According to the result of the nitrogen isotope analysis shown in Figure 5, groundwater in the study area was affected by three sources, including domestic wastewater, farming, and fertilizer. However, it was found that most of the areas were urban zones, indicating that the NO₃ ⁻ in groundwater was from sewage manure or domestic wastewater (Matiatos, 2016). Furthermore, the majority of the groundwater samples in the study area belong to the Ca–Na–Cl and Ca–Cl types, demonstrating that Cl⁻ is a dominant ion in the contaminated water. As mentioned, the NO₃ ⁻ and Cl⁻ detected in the area were mainly derived from anthropogenic activities.

The properties of each soil sample collected in a sugarcane field in this study area (i.e., pH, organic matter (OM), bulk density, hydraulic conductivity, and soil texture) are listed in Table 2 According to the table, the pH values of the samples ranged from 6.80 to 7.90, indicating mildly alkaline conditions. The result was nearly the same as the pH of the groundwater samples collected in this area because of sediments from the weathered limestone and the use of an alkaline pesticide in the sugarcane field. Based on soil texture, there are four different soil types (e.g., clay, clay loam, sandy clay loam, and loam). The OM of the soil samples was in the range of 1.07–2.62%. Only two samples (S5 and S6) indicated OM values higher than 2%. The bulk density of the soil samples ranged from 1.51 to 1.70 g/cm³. Hydraulic conductivity values were in the range of 0.003–0.103 m/d, corresponding to the soil textures. The hydraulic conductivity of all soil samples was not higher than 0.01 m/d, except for S2 and S3, which had values of 0.016 and 0.103 m/d, respectively. In addition, due to an error in soil sample collection, some unreasonable values were obtained from the experiment for estimating hydraulic conductivity; thus, the Neural Network Prediction (NNP) option available in HYDRUS-1D was applied by assigning the values of bulk density as well as the sand, silt, and clay percentages measured from the samples. According to the result, the hydraulic conductivity values of S4, S7, and S8 were defined by HYDRUS-1D. As a result, it can be concluded that most of this area has low water holding capacity, as inferred from soil texture.

A batch adsorption experiment was also used to assess the adsorption coefficients of soil samples collected in the study area. The results of the experiment are presented in Figure 6. According to the results, most of the samples can be fitted with the Freundlich isotherm, indicating that adsorption may occur with heterogeneous active site energy distribution on the soil surface. Therefore, the Freundlich constant (K_f) was suitable for determining K_oc. K_f was used as K_d, as in the following equation (Martins et al., 2018): K oc = K f / f oc .

Sample S6 was identified as the sample with the highest adsorption capacity (K_d = 0.301 l/kg, K_f = 0.822 m³/kg, and QM = 6.575 mg/g) because it also had the highest %OM (2.62%), which is considered to be an important factor influencing the adsorption capacity of the soil. Previous studies have also found that OM was a significant factor in the adsorption of atrazine in soil (Yue et al., 2013; Jing et al., 2020; Yu et al., 2020).

According to the result of the soil analysis, the measured properties displayed that clay content and OM affected the adsorption coefficient (K_d or K_f). A previous study also reported soil OM content played an important role in atrazine adsorption in soil and sediment (Yue et al., 2013; Yue et al., 2017). Moreover, the adsorption isotherms for atrazine of the soils in the study area were described by the Freundlich isotherm with a regression coefficient (R²) between 0.964 and 0.997 as well as a chi-square between 0.066 and 0.696. The linear, Freundlich, and Langmuir constants, representing the adsorption properties of the soils, are presented in Table 3.

Previous studies also reported the adsorption capacity of atrazine in a different isotherm, as illustrated in Table 4. The K_d from other research was in the range of 0.51–4.00 l/kg, whereas the K_d measured in this study ranged from 0.100 to 0.301 l/kg, which is lower than those of the other studies. It can be concluded that the soil samples in this study contained less clay content or lower OM. Additionally, the K_f values derived from previous studies were from 0.60 to 3.90 l/kg, while the K_f measured in this study was partly within the range of the values from other research, with values of 0.284–0.822 l/kg.

In this study, the leaching risk of atrazine was derived from the AF/RF model. The leaching potential of atrazine was classified into five classes: very unlikely (0–0.00001), unlikely (0.00001–0.01), likely (0.01–0.1), moderately likely (0.1–0.25), and very likely (0.25–1) (Khan and Liang, 1989); for large field-scale evaluation, the results are best presented in the form of a map. As illustrated in Figure 7, the leaching potential of atrazine is high in most of the agricultural areas in the western part of the study area due to the low adsorption capacity of the soil (K_oc). This study only focused on agricultural areas which are mostly located in the western part of the area, while soil samples were not densely collected from the urban eastern area. The adsorption capacity in this area was in the range of 0.017–0.121 m³/kg; it was found that different leaching potentials dependent on the soil properties. Soil with low clay content showed low water holding capacity and high hydraulic conductivity (Li et al., 2018). It should be noted that most of the soils with low water holding capacity and low OM appeared to have a high risk for atrazine leaching.

To evaluate the performance of the model, the result of the model was compared with the NO₃ ⁻ concentration detected in the groundwater. Interestingly, Department of Groundwater Resource (2009) detected only a low concentration of atrazine (lower than the detection limit) in groundwater samples collected within the study area as atrazine was strongly adsorbed by soil. It has been reported from studies in several regions that NO₃ ⁻ can be an indicator for pesticide and herbicide leaching to the groundwater table (Kross and Hallberg, 1990; Hallberg, 1997; Vonberg et al., 2014); most of the areas with high leaching potential was also contaminated with high NO₃ ⁻ concentration compared to the area with lower leaching risk. However, the area with high leaching risk might be affected by other factors which were not considered in this model, such as the aquifer horizontal hydraulic conductivity, microbial degradation rate, volatilization, crop-root uptake, and dose of pesticide usage. Furthermore, the GUS score can also be used to evaluate leaching potential. According to Table 5, the GUS scores in this study area were in the range of 3.68–4.83; as a result, the agricultural area could be described as the leachable area because GUS scores there were higher than 2.8.

Additionally, the soil with high clay content expressed a high adsorption coefficient; in contrast, loam and sandy clay loam soils were found to have low adsorption coefficients due to the low clay contents found in such soil samples. Khan (2016) also reported a positive correlation between the adsorption coefficient and clay content. It was also reported that using partial rank correlation coefficient analysis (PRCC) based on f_oc was the most effective leaching potential analysis (D'Alessio et al., 2018). The leaching evaluation map of this study provides an overview for estimating pollution potential; generally, the AF index is used to identify areas with a high potential of groundwater contamination from chemicals. From the result of this evaluation, the area of a high AF index should be monitored first when establishing a groundwater monitoring system. Moreover, the GUS index estimated in this study confirmed that soil in the area had high atrazine leaching potential, meaning that atrazine can leach to the groundwater table as NO₃ ⁻, the pesticide and herbicide indicator, was detected in the subsurface water resource.

Yue et al. (2013) also reported that atrazine in the range of 1.9–54.9 μg/kg was absorbed by soil residue in China. A previous study also reported atrazine concentration in groundwater in an area near the selected area of this study (Thapinta and Hudak, 2003); atrazine concentration was detected in 13 of 90 wells, with the highest concentration in the vicinity of 1.89 ppb. Interestingly, although soil in this area was described as having a high leaching risk as a result of this study, the Department of Groundwater Resources (2009) found atrazine with a low concentration because other factors affecting leaching potential, such as horizontal hydraulic conductivity and specific half-life of each soil sample, was not used for estimating the AF value. Furthermore, it was reported that atrazine with an average concentration of 0.133 mg/kg in topsoil and 0.183 mg/kg in subsoil was detected in soil with a pH of 5 in the Huay Kapo Watershed, Nam Nao District, Phetchabun Province, Thailand (Phewnil et al., 2010). It can be concluded that topsoil had higher hydraulic conductivity, leading to higher atrazine concentration in groundwater as NO₃ ⁻ was a presence with higher concentration. Moreover, as atrazine is slightly alkaline, the lower pH of the soil, the higher the adsorption of atrazine in the soil (Jing et al., 2020). At soil pH of below 7, known as an acidic condition, atrazine may pick up hydrogen ions from the soil solution providing a positive charge to the atrazine molecule increasing the attraction between atrazine molecule and negatively charged soil colloids.