Geostatistics as a basis to the CMLS pesticide simulation model with validation in soil columns

Geoestatística como suporte ao modelo de simulação de agrotóxico CMLS com validação em colunas de solos

SOILS AND PLANT NUTRITION

Geostatistics as a basis to the CMLS pesticide simulation model with validation in soil columns

Geoestatística como suporte ao modelo de simulação de agrotóxico CMLS com validação em colunas de solos

Gilberto Nicolella^{I, *}; Archimedes Perez Filho^II; Manoel Dornelas de Souza^I; Vera Lúcia Ferracini^I

^IEmbrapa Meio Ambiente, C.P. 69 - 13820-000 - Jaguariúna, SP - Brasil

^IIUNICAMP - Instituto de Geociências, Cidade Universitária "Zeferino Vaz" - 13083-979 - Distrito de Barão Geraldo - Campinas, SP - Brasil

ABSTRACT

The use of simulation models is probably the most efficient means for predicting the behavior of pesticides in the soil-plant-water system. The CMLS (Chemical Movement in Layered Soils) simulation model for predicting the fate of pesticides was used for studying the behavior of tebuthiuron, a herbicide used in sugar cane crops, from a sampling grid with 111 sampling points 200 m apart from one another and encompassing three types of soil: Ustic Quartzipsamment, Rhodic Hapludox and Typic Hapludox, all with medium and clay textures. The 373 points assessed by the simulator, generated from samples coming from the original grid and through the geostatistical methods of variography and ordinary kriging, returned the depth values reached by the herbicide after six years of simulation (1989-1995). For the Ustic Quartzipsamment, tebuthiuron, in four simulated points, returned depth values above 43 m and a maximum 50 m, with a certain amount of the product still remaining in the soil that was close to 10% of the original 1.1 kg ha^-1 applied. Results from the column assay used for validating the model showed that the model overestimated the depth reached by the herbicide in 6.6% as compared to the column value for the Ustic Quartzipsamment. The depth was underestimated in 4.5% and 20% for the Typic Hapludox and the Rhodic Hapludox, respectively. These data support the adequacy of the model for assessing the fate of tebuthiuron in both Ustic Quartzipsamment and Typic Hapludox.

Key words: microbasin, sugar cane, tebuthiuron, variography, kriging

RESUMO

O uso de modelos de simulação é provavelmente a maneira mais eficiente para predizer o comportamento de agrotóxicos no sistema solo/água/planta. O modelo de simulação de destino de agrotóxicos CMLS (Chemical Movement in Layered Soils), foi usado para estudar o comportamento do herbicida tebuthiuron, utilizado na cultura de cana-de-açúcar, a partir de uma grade de amostragem composta de 111 pontos amostrais, equi-espaçados de 200 m e englobando três tipos de solo: Neossolo Quartzarênico (RQ), Latossolo Vermelho distroférrico (LVdf) e Latossolo Vermelho distrófico (LVd), texturas média e argilosa. Os 373 pontos avaliados pelo simulador, gerados das amostras da grade original e através dos métodos geoestatísticos da variografia e krigagem, produziram valores de profundidade atingida pelo herbicida, após seis anos de simulação (1989-1995). Para Neossolo Quartzarênico (RQ), o tebuthiuron, em quatro pontos simulados, apresentou valores de profundidade acima de 43 m e com máximo de 50 m, com uma quantidade de produto permanecendo ainda no solo, próximo a 10% do valor aplicado de 1,1 kg ha^-1. Resultados provenientes do ensaio em coluna, usados para a validação do modelo, mostraram que o mesmo superestimou em 6,66% a profundidade atingida pelo herbicida quando comparado ao valor produzido pela coluna, para o Neossolo Quartzarênico (RQ), e subestimou com valores de 4,5% e 20,0%, para o Latossolo Vermelho distrófico (LVd) e o Latossolo Vermelho-distroférrico (LVdf), respectivamente. Esses dados confirmam, para o Neossolo Quartzarênico(RQ) e o Latossolo Vermelho-distrófico (LVd) a adequação do modelo, na avaliação do destino do herbicida.

Palavras-chave: microbacia, cana-de-açúcar, tebuthiuron, variografia, krigagem

INTRODUCTION

Mathematical simulation models have been largely used to evaluate fate of pesticides through soil compartments mainly because of physical, chemical and biological processes involved in the soil-plant-water system. The use of simulation models can provide guidance in determining which pesticides should be used for a particular combination of soil, crop and climate and also in estimating the pesticide level to apply and the most appropriate application time to maximize crop protection and minimize adverse impacts. Likewise, pesticides that should not be approved due to their high potential of contamination of one or more compartments of the environment, can also be identified "a priori" by using these models.

The CMLS (Chemical Movement in Layered Soils) (Nofziger & Hornsby, 1987) model for simulating the fate of pesticides, which was used in this work for assessing the behavior of the herbicide tebuthiuron (N-[5-(1,1-dymethilethyl)-1,3,4-thiadiazol-2-yl]-N-N -dimethylurea), although classified in the literature (Addiscott & Wagenet, 1985) as an instructional model, has been used for various purposes. In association with the DRASTIC model, the CMLS has been used for identifying which combinations of local settings and pesticides can cause the contamination of groundwater (Ehteshami et al., 1991). The CMLS model can also be used as a tool in a system for managing and evaluating the behavior of pesticides in agriculture termed AGCHEMS (Agricultural Chemical Evaluation and Management System). This system is used for assessing the impact of various scenarios of pesticide management (Haan et al., 1993).

The CMLS model was also used in association with a climate simulation model named WGEN (Weather Generation) for evaluating the impact of possible and equally probable sequences of climate events on the movement time and the relative amount of pesticides that penetrate the soil to a certain depth (Haan et al., 1994). Furthermore, the ability to predict experimental results obtained in lysimeters (Nofziger et al., 1994), to simulate the concentrations of herbicides (which are measured in samples of drainage water), and to evaluate the impact of the CMLS model itself on the feeding of data to resolution maps, are the other relevant aspects addressed by the CMLS. Temperature effects on the atrazine degradation process in a loamy sand under nonirrigated summer crops, were introduced in the CMLS model (Wu & Nofziger, 1999).

The main purpose of this study was to characterize and predict the migration of tebuthiuron, a herbicide widely used in Brazilian agriculture, mainly in sugar cane crops, in order to evaluate possible contamination of ground water.

MATERIAL AND METHODS

The experimental area selected for this study is located in the microbasin of the Espraiado pond, which is situated at the border of the municipalities of Ribeirão Preto, Cravinhos and Serrana, SP, Brazil, with a total area of 4,140 ha. The area is situated between 21°05' and 21°20' S and 47°40' and 47°50' W. The average altitude is 600 m above sea level, and the dominant relief is of the light wavy type. Soils are predominantly Typic Hapludox, Rhodic Hapludox next to hillside and Ustic Quartzipsamment in the lowlands and downstream of the microbasin. The original vegetation is composed of tropical to low deciduous forests. The geology of the Espraiado microbasin is composed of two lythotypes: the basaltic rocks from the Serra Geral Formation both midstream and upstream, and the arenites from the Botucatu Formation downstream to the mouth of the river. From the hydrological point of view, both the regional climatology and hydrogeology are worth mentioning. The climate of Ribeirão Preto is tropical with a dry winter, which is typical of the savanna (AW) according to Köppen's classification. The average annual temperature oscillates between 21 and 22^oC and the annual precipitation ranges from 1,300 to 1,500 mm. In relation to the regional hydrogeology, the Guarani aquifer with an area of 16,000 km² is the most relevant. The aquifer recharge occurs only in the outcropping areas of the Botucatu and Pirambóia formations (DAEE, 1974; Gomes, 1995).

To meet the requirements of the CMLS solute transport model, both deformed and undeformed soil samples were extracted from October 1995 to July 1996, from a sampling grid composed of 111 equispaced sampling points (200 to 200 m), at depths of 0-20, 20-40, 40-60, 60-80, 80-100 and 100-120 cm. The physical and chemical parameters measured were the bulk density, total porosity, field capacity, wilting point and organic carbon.

The pesticide adsorption study was conducted on the three most frequently occurring soils of the Espraiado microbasin in Ribeirão Preto, SP. Samples of 2 mm sieved air-dried soil of Typic Hapludox, Rhodic Hapludox and Ustic Quartzipsamment, of depths of 0-10 cm and 10-20 cm were used for the experiment. Results are shown in Table 1.

To determine the adsorption isotherms, 5 g triplicate samples of each soil were put into 50-mL polypropylene capped tubes containing 25 mL of 0.01 mol L^-1 CaCl₂ solution with the following concentrations of the active ingredient (tebuthiuron, commercial formulation, 500 g i.a.L^-1 ): 0, 1,2, 4, 8 and 14 mg L^-1. The tubes were shaken at a speed of 170 oscillations per minute during 24 hours at 24^oC. Following agitation, samples were centrifugated at 3,000 rpm for 10 minutes. The overflow was then filtered through a 47 mm diam x 0.45 mm porosity membrane (ME 25, Scheleicher & Schüll). The filtrate was analyzed using an LC 10 AD Shimadzu liquid cromatograph with a SPD 10AV to 254 nm ultraviolet detector. The column used was a C18 Bondesil (4.6 mm × 25 cm × 5 mm) with a 0.8 mL min^-1 flow, MeOH:H₂O (63:37; v/v) moving phase, 20 mL injection volume. Data obtained were adjusted to the Freundlich equation. The equation's K_f was used as an estimator of K_d.

The leaching tests were carried out in 15 cm diam × 80 cm tall PVC columns, one being used for each soil type. Columns were filled with soil in the same sequence of layers as those found in the field. Air-dried fine soil with densities similar to those found on the field was used (1.27 g cm^-3 average soil density for the Typic Hapludox, 1.24 g cm^-3 for the Rhodic Hapludox and 1.56 g cm^-3 for the Ustic Quartzipsamment). After filling the columns, soil solution extractors were placed every 10 cm of depth of the column. Extractors were made of 12.52 mm porous capsules glued to a 20 cm long × 12.52 mm diam PVC tube. The PVC tube was sealed with a rubber cork and a 3.17 mm diam hole was punched into the center of the cork through which a nylon spaguetti tube was introduced into the capsule until reaching the bottom. The other end of the nylon tube was introduced into a two-holed cork of a collection tube (200 mL amber tube). A 20 cm long piece of the spaguetti tube was introduced into the other hole through which vacuum was created with the help of an electric pump, to extract soil solution.

The columns were previously saturated through a hydraulic load and an ascending flux. After saturation, the conductivity of the saturated soil was estimated by Klute's method (1965). The excess water was gravitationally drained for three days and 1,000 ml of tebuthiuron solution containing 10.36 mg L^-1 (equivalent to two times the in-field dosage of 2.5 kg ha^-1) of the active ingredient were then applied at the top of each column. The volume of solution used corresponds to a 57 mm water depth.

Based on the values of hydraulic conductivity obtained for the saturated soils, the time for each solution to reach the base of the column was calculated and the vacuum pump was then connected, starting the collection. After the extraction was completed, the solution was filtered and the filtrate was analyzed by the process described above in order to measure the tebuthiuron concentration.

For each parameter of soil attributes and for the six depths mentioned above, the level of spatial dependency was studied so that one could become familiar with the covariance structure and conduct the kriging, i.e., increase the number of points in the sampling space and, from an amplified set of points, use the CMLS simulation model for determining the fate of pesticides in order to create the agrotoxic load map for the soil. The geostatistical analyses of variography, cross validation and kriging (Cressie, 1993; Vieira, 1982; Vieira et al., 1983; Isaaks & Srivastava, 1989) were made using a geostatistical software (Englund & Sparks, 1991) while the simulation model was calculated for each of the 111 sampling points to which the interpolated points were added via ordinary kriging, with the interpolated points being placed 100 m away from the original points. Plateau-based theoretical methods, namely the spherical, exponential, Gaussian and nested models, were used for adjusting the data.

The CMLS Model

The CMLS (Chemical Movement in Layered Soils) model used to simulate the movement of a chemical through the soil profile, estimates the depth of the center of mass of a non-polar chemical, as a function of time after application. It assumes that chemicals move only in the liquid phase in response to soil-water movement. Since the chemical moves with soil water, the change in depth of the solute depends upon the amount of water passing through the soil. If D_c(t) represents the chemical depth at time t, q the amount of water passing that depth during time dt, the depth of the chemical at time t + dt is given by:

where q_fc is the volumetric water content of the soil at "field capacity" and R is the retardation factor for the chemical in the current layer of the soil.

Assuming a linear and reversible equilibrium adsorption model, the retardation factor is written as:

where r is the soil bulk density and K_d is the partition coefficient for the chemical in the soil. K_d was calculated as the product of the organic carbon partition coefficient, K_oc, and the organic carbon content of the soil, according to the equation:

Some other assumptions in chemical movement assumed by CMLS deals with aspects related to soil, and that chemical properties are uniform within a soil layer, all water in the profile takes part in the flow process, water already in the soil profile being pushed ahead of infiltrating water in a piston-like manner and dispersion of the chemical being ignored.

Evaluating depth chemical implies in estimating q at each day and calculating the chemical depth at that day, considering the depth of the preceding day and the change in depth for that day. The process is then repeated every day during the period of time contemplated by the simulation (Haan et al., 1993; Nofziger & Hornsby, 1987).

RESULTS AND DISCUSSION

An exploratory analysis of the data, as a first indicator of the hypotheses of stationarity (which is fundamental for the interpolation of data, i.e., for the kriging procedure) could or not be accepted, showed that in a certain direction in which extreme values were predominant, and for a given number of lag, the semivariogram tended not to stabilize, that is, it did not reach a plateau and, consequently, the intrinsic hypothesis was not proven true. The solution outlined for this problem was based on a work by Burgess and Webster (Burgess & Webster, 1980) according to which the tendency of non-stationarity large lag numbers does not always significantly affect interpolation by kriging as long as a local stationarity exists along short distances for which the semivariogram is used. Therefore, as short distances it was decided to take those for which |h|< (L/2), where L is the maximum distance within the incomplete grid. Information generated from the directional semivariograms showed that in general the maximum direction of continuity was located around 135^o and the minimum direction of continuity, at 45^o. Thus from the assumed simplification, the modeling approach consisted of considering the two possible variability models, that is, the isotropic or omnidirectional model, and the anisotropic or geometrical model, in the main directions of the sampling grid, that is 0°,45°,90 ° and 135°.

Both the Gaussian (GAUSS) and spherical (SPH) theoretical models have proven to be the most appropriate for adjusting estimated semivariances (Table 2). The variation of the amplitude of reach was situated between 3.0 (600 m) and 7.60 (1520 m) for organic carbon at 20-40 cm and 0-20 cm, respectively. The predominant variability structure for the thirty variables was the isotropic one, except for five cases.

The cross validation technique, which consists of estimating values at each sampled position in the area through ordinary kriging with neighboring sampling values excluding the value of the point being estimated, was used for comparing the omnidirectional semivariogram models to the anisotropic models at 0^o and 90^o and 45^o and 135^o, which involved different kriging strategies with respect to search radius, maximum number of points per sector and minimum number of points to be used. In view of the peculiarities of the sampling grid and the simplifications assumed for the variography, a search radius smaller or equal to three units (600 mm) was adopted, with the maximum number of points always being less than 10 and the minimum equal to 3. The statistics generated by cross validation is shown in Table 3. It refers only to bulk density (x10²) for depths of 0-20 and 20-40 cm.

An observation of Table 3 and comparison of the first two columns for both depths, shows that the estimated values have less variability than the observed values, an expected result since the estimated values in the cross validation are the average of the neighboring values and their oscillations are hardly attenuated. Both the average and the variance values for the standard differences or errors are close to 0 and 1 respectively for both depths, which meets two of the criteria used for considering the validity of the adjustment model.

Once the parameters for the most appropriate semivariogram model were selected through cross validation, the ordinary punctual kriging followed to generate a set of points resulting from the interpolation. The maximum number of neighboring points considered was ten, which were situated inside a circle (isotropic model) with a search radius equivalent to half the reach around the point to be estimated. Due to the fact that the sampling grid was equispaced, albeit incomplete, and a reasonable number of points was enough for establishing the scenarios for contamination risk by tebuthiuron, the interpolated mesh was composed of points placed at an equal distance of 100 m from one another, totaling 263 estimated values.

The procedures adopted for the selection of the best adjustment model were the criteria for the standardized residue of Y as a function of X, which must be situated in the 2 to + 2 interval, and of the normal probability of residue according to which the normal expected and actually observed residues are compared.

Both the standardized residue and the normal residue probability criteria were used for all data and except for the case of the Ustic Quartzipsamment at the depth of 0-10 cm, Freundlich's model showed the best adjustment.

According to Freundlich's model, when the values of (1/n) are numerically equal to the coefficient b of the linear equation and both are close to or equal to 1, the value of K_f is equivalent to a solute partition coefficient between the solution and the solid surface, that is, K_f = K_d (adsorption coefficient). With the exception of the Ustic Quartzipsamment at 0-10 cm of depth, which presented a value far from 1 as an estimate for (1/n), the other values found in Table 5 meet this assumption.

The highest value of K_d for the Rhodic Hapludox points to a higher degree of herbicide adsorption in relation to the Ustic Quartzipsamment, which has a lower value for this coefficient, presenting, therefore, a higher degree of leaching of the chemical.

The values of K_d at 10-20 cm depth would naturally be chosen to represent this coefficient in the 0-20 cm interval for the three considered soils. However, the values actually used in the CMLS simulator are described and justified below.

The Column Experiment

During the column assay, after estimating the remaining concentration of tebuthiuron, all the K_d values were tested to validate the model, and only for the Rhodic Hapludox a correspondence was found between the K_d that best adjusted to the adsorption data (Freundlich's model) with closer values generated by both, the simulator and the column. For the other two types of soil, namely the Ustic Quartzipsamment and the Typic Hapludox, the average values of K_d obtained for the linear model at the depths of 0-10 and 10-20 cm were the closest to the validation results. The supposition that the CMLS model, with reference to the adsorption process, is represented by a linear model that is in equilibrium and reversible, can explain this fact, as long as the models of Freundlich and Langmuir are not placed in this category, although there is a possibility that they can be linearized. The results of data obtained in both the columns and the simulator are shown in Table 6.

With reference to depth (cm) and the Ustic Quartzipsamment, the model overestimated in 6.66%. With regard to the Typic Hapludox and the Rhodic Hapludox, the model was overestimated returning values corresponding to 95.5 and 80.0% respectively of the values generated by the column (Table 6). With respect to the remaining concentration (mg L^-1) and to the Typic Hapludox, the model was overestimated in 1.98%, and underestimated for both the Ustic Quartzipsamment and the Rhodic Hapludox, which presented values corresponding to 96.36% and 95.33%, respectively, of the values generated by the column.

Simulation of the CMLS Model

Basic statistics resulting from the application of the CMLS model for 373 points encompassing the original observations of the sampling grid and those produced by interpolation via kriging are summarized in Table 7 for the Rhodic Hapludox, Typic Hapludox and Ustic Quartzipsamment.

The Ustic Quartzipsamment is the soil that poses the highest risk of contamination of the water table by the herbicide tebuthiuron since the maximum depth value of 50 m reached by the product was obtained only 4 years after the simulation took place (1989-1993). In addition, four more points (16.5; 5), (16.5;1), (18;0) and (9;3) presented depth values above 39 m, which is not regarded as worrisome since these points are located next to the drainage water, thereby posing a high risk of contamination of the water table. The other two soils, namely the medium and clay textured Typic Hapludox and the Rhodic Hapludox can be regarded as having low risk of contamination, particularly the Rhodic Hapludox which, due to its intrinsic characteristics, opposes resistance to leaching of the product. The depth value of 29 m reached at point (12;0) corresponding to the Typic Hapludox can also be attributed to the proximity of the water table of the microbasin of the Espraiado Pond and consequently the point of the Ustic Quartzipsamment.

The simulation results must not be regarded as definitive but rather as an indication of the high contamination risk of rivers receiving this water. The CMLS simulator has proven to be an adequate management tool due to its low demand for parameters and easiness to handle when compared to the more sophisticated simulators used in research models. However, in order to corroborate to the simulation results, hydrogeological and geological aspects, and other must be assessed, serving as a reference for the final conclusions about the contamination risk of a certain compartment of the environment by any product that is potentially a pollutant.

Received October 31, 2003

Accepted November 01, 2004

• ADDISCOTT, T.M.; WAGENET, R.J. Concepts of solute leaching in soils: A review of modeling approaches. Journal of Soil Science, v.36, p.411-424, 1985.
• BURGESS, T.M.; WEBSTER, R. Optimal interpolation and isarithmic mapping of soil properties. I. The semivariogram and punctual kriging. Journal of Soil Science, v.31, p.315-331, 1980.
• CRESSIE, N. Statistics for spatial data. 1.ed. New York: John Wiley, 1993. 900p.
• DEPARTAMENTO DE ÁGUAS E ENERGIA ELÉTRICA. Estudo de águas subterrâneas: Região Administrativa 6 - Ribeirão Preto. São Paulo, 1974. v.1. Resumo.
• EHTESHAMI, M.; PERALTA, R.C.; EISELE, H.; TINDALL, T. Assessing pesticide contamination to groundwater: a rapid approach. Ground Water, v.29, p.862-868, 1991.
• ENGLUND, E.; SPARKS, A. GEO-EAS. 1.2.1 Assessment Software. User's guide. Las Vegas: EPA, 1991. (US-EPA, 600018).
• GOMES, M.A.F. Impacto da agricultura intensiva de cana-de-açúcar sobre as propriedades do solo e da água subterrânea. In: ENCONTRO ANUAL DA SECÇÃO BRASILEIRA DA INTERNATIONAL ASSOCIATION IMPACT ASSESSMENT, 4., Belo Horizonte, 1995. Anais Belo Horizonte: IAIA, 1995. p.383-387.
• HAAN, C.T.; NOFZIGER, D.L.; GREGORY, M. An agricultural chemical evaluation, and management system. Journal of Agricultural Engineering Research, v.56 p.3301-312, 1993.
• HAAN, C.T.; NOFZIGER, D.L.; AHMED, F.K. Characterizing chemical transport variability due to natural weather sequences. Journal of Environmental Quality, v.23, p.349-354, 1994.
• ISAAKS, E.H.; SRIVASTAVA, R.M. An introduction to applied geostatistics 1.ed. New York: Oxford University Press, 1989. 560p.
• KLUTE, A. Laboratory measurement of hydraulic conductivity of saturated soil. In: BLACK, C.A. (Ed.) Methods of soil analysis Madison: ASA, 1965. p.210-221.
• NOFZIGER, D.L.; HORNSBY, A.G. Chemical movement in layered soils: user's manual. Gainesville: University of Florida, Institute of Food and Agricultural Sciences, Cooperative Extension Service, 1987. 44p. (Circular, 780).
• NOFZIGER, D.L.; CHEN, J.S.; HAAN, C.T. Evaluating the chemical movement in layered soil model as a tool for assessing risk of pesticide leaching to groundwater. Journal of Environmental Science and Health, v.29, p.1133-1155, 1994.
• VIEIRA, S.R. Geostatistical analysis of some agronomical observations. Davis: University of California, 1982. 261p. (Ph.D. Thesis).
• VIEIRA, S.R.; HATFIELD, J.L.; NIELSEN, D.R.; BIGGAR, J.W. Geostatistical theory and application to variability of some agronomical properties. Hilgardia, v.51, p.1-75, 1983.
• WU, J.; NOFZIGER, D.L. Incorporating temperature effects on pesticide degradation into a management model. Journal of Environmental Quality, v.28, p.92-100, 1999.

Publication Dates

• Publication in this collection
  22 Feb 2005
• Date of issue
  Jan 2005

History

• Accepted
  01 Nov 2004
• Received
  31 Oct 2003