MODELLING OF SOLUTE TRANSPORT THROUGH AN UNSATURATED ZONE EXTENDING FROM GROUND SURFACE TO WATER TABLE

Abstract

Waste disposal on land and application of fertilizers and pesticides to crop lands has become a common practice universally. Water infiltrating at the ground, dissolves such matter and carries it downward through the unsaturated zone. Many types of waste material (e.g., heavy metals, radioactive material) do not decompose easily. Such pollutants travelling through the unsaturated zone join the water table and may affect the water quality adversely. Further, fertilizers which are not utilized by crops are transported below the root zone by percolating water and pose a potential threat to the groundwater quality. In a reverse situation, evapotranspiration may lead to an accumulation of pollutants in the root zone. This may lead to a fall in crop yield and deterioration of top soil conditions. To avoid such problems and to design safe disposal systems, the time variant rate of pollutant transfer to the water table as well as depth and time variant concentrations need to be estimated. In the present study an attempt has been made to develop a numerical model for simulating one dimensional (vertical) solute transport from ground to the water table. The mechanisms of solute transport accounted for are convection, hydrodynamic dispersion, lateral diffusion into/out of an immobile phase (in case of two phase solute transport) and linear adsorption - desorption in either or both the phases. The model is developed in the following four stages. Stage I Single phase non-reactive solute transport. Stage II Single phase reactive solute transport, accounting for first order linear kinetic adsorption-desorption. Stage III Two phase non-reactive solute transport. Stage IV Two phase reactive solute transport accounting for linear equilibrium adsorption-desorption. The numerical methods employed for solving the solute transport equations are the method of characteristics (MOC) and the finite differences. The convective component is solved using MOC to overcome the problem of numerical dispersion encountered in solving convection dominated flow problems. Change in concentration due to hydrodynamic dispersion and adsorption-desorption is accounted for subsequently, using an implicit finite difference scheme. The soil moisture and flux distribution required for solving the solute transport equations is obtained by solving (Mohan Rao, 1986) the head form of Richards equation using a Crank-Nicolson finite difference scheme. The problem of non linearity arising due to dependence of scpecific moisture capacity and capillary conductivity on soil moisture (or capillary head) was taken care of by using Picard's Iteration method. To account for solute transport due to convection by MOC the domain under consideration is discretized by a finite number of moving packets of a pre-assigned strip thickness. Each moving packet is defined by two co-ordinates (representing its upper and lower bounds) and the solute and water volumes contained in it. During simulation, the movement of these packets is traced. During each time step the new positions of moving packets are obtained by ensuring a compatibility between the cumulative water profiles obtained by the considerat of flow and transport. Solute volumes (per unit plan area) of these moving packets were further redistributed amongst themselves to account for solute transport due to hydrodynamic dispersion. This is done by solving the ii governingdifferential equation using an implicit finite difference scheme. To compute further change in concentration due to adsorption-desorption of solute by the soil matrix or lateral diffusion of solute into/out of the immobile phase and subsequent adsorption-desorption a fixed grid system is superposed on the moving co-ordinate system. Concentration distribution of this grid is computed by identifying moving packets lying wholly or partially in the area of influence of any node. The governing differential equations are then solved using an implicit finite difference scheme. Further, change in solute volume (per unit plan area) at the nodes is attributed to the moving packets. Thus, the model is capable of simulating spatial and temporal distribution of solute concentration and quantifying the volume of solute (per unit plan area) joining the water table. The model was validated by comparing model simulated solute transport with the results of two analytical solutions of van Genuchten and Alves, 1982 (cited in Parker and van Genuchten, 1984) and Parker and van Genuchten, 1984. The analytical solution of van Genuchten and Alves', 1982 (cited in Parker and van Genuchten, 1984) pertains to flow conditions, accounting for linear equilibrium adsorption-desorption. Neglecting solute matrix interaction this solution was used to validate the stage I model. Results obtained by the two methods showed an excellent agreement. The analytical solution of Parker and van Genuchten (1984) pertains to single phase and two phase solute transport under steady state flow conditions. For single phase solute transport the sorption sites present in the soil matrix are assumed to comprise of two fractions i.e., equilibrium adsorption ('type-1' sites) and kinetic iii equilibrium adsorption ('type-2' sites). For two phase solute transport the interaction between solute and soil matrix in both phases, is described by a linear equilibrium adsorption-desorption isotherm. By assigning appropriate values to the parameters, this solution was used to validate the stage II, III and IV model. An excellent agreement was obtained for most of the simulations. The proposed model was also used to simulate reported experimental data of two field experiments (Warick et al, 1971; Bottcher and Strebel, 1989). The model (Stage I) was used to simulate Chloride concentration profiles in depth under conditions identical to the experiments of Warrick et al. (transport of CaCl and water in Panoche clay loam). The simulated and measured concentration profiles compared reasonably well, except for a lag between the simulated and measured depth of solute travel. The simulation was repeated considering the presence of an Immobile phase (Stage III model) and neglecting solute transfer into/out of the immobile phase (equivalent to a case of anion exclusion, assuming the effect of osmotic potential on the fluid flow to be negligible). For 6 = 0.06, the lag was almost eliminated. A bromide leaching experiment was conducted by Bottcher and Strebel, 1989 (unpublished data) and breakthrough curves at 51 locations at depths of 120 cm were measured. Profiles of bromide amounts in depth at 26 boring sites were also made on two dates. The measured experimental data exhibited a considerable lateral variation in solute transport. Although, the proposed model does not account for horizontal transport, a reasonable agreement was observed between the model simulated and measured mean concentration distribution. IV Model application to real life problems was demonstrated by simulating two problems of solute transport through the unsaturated zone. The model was used to simulate salt accumulation in the root zone of two crops (wheat and rice) assumed to grow over a period of one year and irrigated by considerably saline water (1.5 mmho/cm). Concentration profiles (ground to water table) at different discrete times, covering the entire period were also simulated. Evapotranspiration by the crops was accounted for (Doorenbos et al., 1979). Salt accumulation in the root zone was also estimated using the salt storage equation (Van der Molen, 1973). A considerable deviation was observed between the salt accumulation as computed by the salt storage equation and by the model. This deviation was possibly caused by the assumptions on which the salt storage equation is based. The model was used to simulate solute travel of a conservative pollutant (assumed to be abundantly available on the ground) in two types of soils (loam and clay) under conditions of heavy monsoon rainfall. At the end of the simulation period (140 days) the solute joining the water table in case of clay was negligible - caused only by dispersion. However, in case of loam, the convective front reached the water