Modeling Flow and Transport in Unsaturated Porous Media : A Review

Underground water is a vital natural resource and every effort should be made to understand ways and means of efficiently using and managing it. The unsaturated zone, bounded at its top by the land surface and below by the water table, is the region through which water, together with pollutant carried by the water, infiltrates to reach the groundwater. Therefore, various processes occurring within the unsaturated zone play a major role in determining both the quality and quantity of water recharging into the groundwater. Classical methods of predicting water flow and contaminant transport processes in unsaturated porous media are generally inadequate when applied to natural soils under field conditions, due to the occurrence of macropores, structured elements and spatial variability of soil properties. Contaminant transport models also require the simultaneous solution of the unsaturated flow and transport equations. For applications to field conditions, numerical solutions and computer simulations based on numerical models have been increasingly used. Advances and progress in modeling water flow and contaminant transport in the unsaturated zones are reviewed, and specific research areas in need of future investigation especially relevant to Oman are outlined. MODELING FLOW AND TRANSPORT IN UNSATURATED POROUS MEDIA

O man is an arid country with limited natural water resources.The rainfall is unreliable and records indicate periods of several consecutive dry years.There are few perennial streams; surface flow is usually confined to periods of a few days or even a few hours after storms.Infiltration through streambeds of surface run off from the mountains is the main source of recharge to the groundwater.The main user of water in Oman is agriculture.The climate will support very little rainfed agriculture and irrigation is essential.Water for domestic use is traditionally obtained from aflaj, springs or shallow wells.
In most traditional communities in Oman, groundwater has long been extracted and directly utilised through the falaj (plural aflaj).It may consist of simple conduits conveying water from wadi channels to domestic use and irrigation systems, or consist of tunnels constructed at the water table to exploit the groundwater resource.The falaj is self-regulating; discharge declines as the upstream storage in the aquifer falls during droughts, and rises after the aquifer is recharged.In coastal regions, water is extracted from shallow wells for various purposes.
Desalination plants have also been established to provide domestic supplies for some coastal areas; but agriculture, which is responsible for more than 95% of freshwater consumption, is almost entirely reliant upon groundwater.Since the Sultanate's water resources are fully utilised, any pollution may have immediate and serious consequences.Therefore, every effort should be made to ensure that all valuable aquifers are not contaminated.
Increased public awareness and concern over many environmental problems is, at least partly, attributed to the fact that these problems are so visible.However, the situation is quite different for groundwater, which is hidden from sight, and renders the characterization of impacts on its quality and availability extremely difficult.Consequently, problems of over-use, pollution or saltwater intrusion may develop to a critical level before being noticed.In the absence of any reliable scientific knowledge applicable to groundwater, this valuable resource cannot be effectively protected, conserved and managed.
A quantitative study of water flow and contaminant transport in the unsaturated zone is a key factor in the improvement and protection of the quality of groundwater supplies.This is the region bounded above by the land surface and below by the ground water table.It is the region through which water derived from precipitation and irrigation infiltrates and transports contaminants to reach the groundwater.Over-pumping of groundwater for domestic, agricultural and industrial consumption will lower the water table and can accelerate the movement of pollutant-laden surface water into the groundwater.
In unsaturated soil, the water content is less than the soil porosity, and the soil water pressure head (matric potential) is negative, being less than that of free water at the same location.The unsaturated zone is inextricably involved in many aspect of hydrology: infiltration, evaporation, groundwater recharge, soil moisture storage, and soil erosion.It also contributes to the spatial and temporal distributions of plant communities under naturally occurring rainfed conditions and serves as a modifying influence on the production of cultivated crop species.Thus the unsaturated zone represents the conduit through which liquid and gaseous constituents are attenuated and transformed as they are exchanged in both directions between the ground surface and the ground water table.
Traditional deterministic modeling approaches to determine downward movement of contaminants, based upon the convection-dispersion equation for contaminant transport and the Richards equation for water flow, work relatively well for homogeneous field soils and laboratory soil columns.Experimental investigations, however, have shown that flow and transport processes in most fields are heterogeneous.Thus, to use classical approaches for describing transport processes in natural soils under field conditions would be a misapplication of the methods.Soil heterogeneity, occurring in structured, cracked or macroporous soils, and due to the inherent spatial variability of soil properties, has a distinct effect on the water flow and contaminant transport processes.Since macropores act as preferential flow paths, occurrences of macropores at the field scale may increase the risk of leaching of pollutants to the groundwater.Therefore, a quantitative prediction of contaminant transport in macroporous, heterogeneous soils at field or regional scales is an important matter in dealing with environmental problems.
This review is an attempt to summarize some aspects of the main scientific advances and progresses in modeling water flow and contaminant transport in unsaturated porous media.Although water flow and contaminant transport phenomena should be treated simultaneously, for convenience, we first discuss water flow and then deal more specifically with contaminant transport processes.We focus primarily on conceptual and mathematical, deterministic aspects of unsaturated zone flow and transport processes.
The review starts with the basic physical concepts that are classically used to describe water flow and solute transport processes in unsaturated porous media.Existing analytical and numerical solutions of the unsaturated water flow and contaminant equations will be discussed.For applications to field conditions, numerical solutions and computer simulations based on numerical models have been increasingly used.Thus recent progress in numerical methods employed for these purposes will be examined.
Oman has extensive areas of active and fossil karst.The complex geological structure of Oman consists of well-developed fracture systems in the rocks forming the northern mountain aquifers.The highest average annual rainfall for the Sultanate is in the mountainous region with more than 300mm.Perennial spring discharge from the limestones in the northern mountains can maintain surface flow over a substantial distance.The essence of fractured aquifers is the existence of an interconnected set of voids capable of conveying water and contaminants at speeds much higher than the typical Darcian flow rate.As classical representations are not intended to describe flow and transport processes in natural soils under field conditions where the occurrence of macropores, fractures, karstic conduits and spatial variability of soil properties is dominant, more research is still needed for a quantitative prediction of water flow and contaminant transport in the unsaturated porous media.

Mathematical Equations of Water and Transport in Unsaturated Soils
In this section, the general classical deterministic forms of the water flow and contaminant transport equations in the unsaturated zone are described.A comprehensive derivation is provided by Bear (1972), Pinder and Gray (1977), and Freeze and Cherry (1979).Note that there is considerable skepticism as to the validity of these equations, and Gray and Hassanizadeh (1991) proposed a new set of equations to describe the unsaturated flow processes obtained from averaging theory coupled with an interface thermodynamic analysis.However, this set of equations contains many more unknowns than that of classical equations.
The liquid in the unsaturated zone consists of a solution of water and dissolved solid and gaseous constituents.Therefore, the forces acting on the liquid cannot be restricted merely due to the earth's gravitational field (Nielsen et al., 1986).However, since direct measurement or calculation of other forces still remains undeveloped, most studies to date assumed that the net force acting on the liquid is the matrix potential and gravitational force.Thus, the rate at which water moves through the unsaturated zone is according to the Darcy's velocity,

K z h + h
where is the hydraulic conductivity, Ψ the (hydraulic head) total potential which can be written as in which is the (pressure) matric potential head, and z is the gravitational potential head, measured positively upward.In unsaturated flow, the void space is partly filled by air and partly by water.As air is infinitely mobile, it is generally assumed that the trapped air in the void space is at constant atmospheric pressure.

= Ψ
Combining the Darcy's velocity with the conservation of mass equation leads to Richards equation for water flow in the unsaturated zone.In the mixed form, the unsaturated flow equation is written as where θ is the volumetric water content, t is time, z is the (positively upward) vertical distance, and G represents sources and sinks of water in the system.Other two alternative forms of Richards equation are the pressure head h-based equation where the specific moisture capacity C is the slope of the soil water retention curve ( ) where it is customary to define the coefficient of moisture diffusivity as . Although these three forms are mathematically equivalent, some give better results when numerical approximation methods are applied.
Note that these equations show that the soil water retention curve ( ) h θ and the unsaturated hydraulic conductivity K are essential components for predicting water flow in unsaturated zones.Retention curves show how much is retained in the soil against gravity; its shape depends on the soil pore-size distribution and pore shapes of the porous medium.However, its derivation ignores soil matrix and fluid compressibilities and assumes that the fluid density ρ is a function of the pressure only, and is independent of concentration.Equation (3) also assumes that Darcy's law is valid for unsaturated flow.
The hysteretic nature of the measured soil water retention curve ( ) h θ implies that the soil water content is not a unique function of h, but depends on the previous history of the soil (Parlange, 1980).The effects are more pronounced by the presence of entrapped air, and by alternate rates of wetting and drying (Bear, 1972).This effect introduces a nonlinearity into the Richards equation (3); and hence, it is necessary to state whether a drying or wetting process is taking place along the boundary of the system (Hills et al., 1989a and1989b).Hysteresis is rarely considered under field situation, due to the difficulties in measurement.
Numerous methods have been developed to determine the soil water retention and hydraulic conductivity functions of unsaturated soils using both in situ field and laboratory procedures (van Genuchten and Nielsen, 1985;Bird et al., 1996, Assouline et al., 1998).While in situ field measurements are the most representative of actual flow conditions, current methods are likely to remain approximate in nature (Kool et al., 1987).This is due to simplifying assumptions inherent in the methods, and to problems of obtaining undisturbed samples.Kool et al. (1985) developed a method to determine water retention and hydraulic conductivity functions simultaneously from transient flow data by parameter estimation.The methodology appears promising, although problems of computational efficiency, convergence, and parameter uniqueness remain unresolved (Eching et al., 1994;Wildenschild et al., 1997).
Direct field methods to determine the unsaturated hydraulic conductivity are time consuming, expensive, and usually subject to simplifying assumptions.An alternative method is the theoretical calculation from more easily measured field or laboratory soil water retention data (Childs and Collis-George, 1950;Mualem, 1976;Arya and Paris, 1981;van Genuchten and Nielsen, 1985;Assouline et al., 1998).This method produces analytical (non-tabular) functions that are particularly useful for inclusion in simulation models and also for a rapid comparison or scaling of the hydraulic properties of different soils.Among the retention equations that have realistic shapes is the equation proposed by van Genuchten (1980)

(
) The physical processes that control the solute distribution in the unsaturated zone are advection by water flow and dispersion by mixing and molecular diffusion.Other processes have been noted by Nielsen et al. (1986), but only summarily known, if understood at all.The classical equation that is thought to describe solute transport ( ) ( ) where c and s are solute concentrations associated with the solution and solid phases of the soil, D is the dispersion coefficient, and P is the rate of solute removal or supply not specifically included in s.Loss and gain of solute mass can occur as a result of chemical reactions or radioactive decay.
The first term of the equation ( 7) describes the rate at which the solute interacts or exchanges with the solid matrix.Both equilibrium and kinetic rate laws have been used to describe adsorption-exchange processes (Nielsen et al., 1986), which are traditionally represented by the retardation coefficient, R. In practice, is used as an empirical parameter that includes all of the solute spreading mechanisms.However, in concept, the coefficient is commonly assumed to include only the effect of mechanical dispersion and molecular diffusion (Bear, 1972).

D
The flow and transport equations previously described are found not adequate to describe field-scale spreading of contaminants transported in unsaturated soils due to the occurrence of macropores and spatial variability of soil hydraulic properties.Results from field experiments (Richards and Steenhuis, 1988) have indicated that groundwater contamination may not be associated with the unsaturated flow.Pollutant could be transported very expeditiously through certain pathways in unsaturated soils into the groundwater (Ju and Kung, 1997).Nevertheless, the scope and impact of this preferential flow on contaminant transport under field conditions have not yet been fully characterized (Flury, 1996).This is mainly because it is difficult to obtain representative patterns of preferential flow with soil coring methods (Ghodrati and Jury, 1990).
The irregular and erratic character of hydraulic conductivity and dispersion coefficient distributions observed in natural soils has led to the development of a stochastic framework for transport modeling.In this stochastic-convective approach, contaminant is assumed to move at different velocities in isolated stream tubes without any mixing of contaminant between the tubes (Dagan and Bresler, 1979).
In light of the uncertainty in the preferential flow mechanism and hence inadequate deterministic predictions, Jury and Gruber (1989) developed a stochastic approach to assess the groundwater contamination.This model estimates the probability of a contaminant reaching groundwater, using the same deterministic equations that are used in the homogeneous unsaturated soils, but assuming the parameters of the equations and the variables are stochastic functions (Yeh, 1992;Feyen et al., 1998).Another approach involves extrapolation techniques to get large scale effective parameters from the small-scale measurements (Wagner and Gorelick, 1986;Vanderborght et al., 1996).
A different approach involves the dual-porosity model proposed by Gerke and van Genuchten (1993), which is gaining widespread acceptance for modeling preferential water flow and solute transport.The two domains consist of the soil matrix and macropores.Flow and solute transport equations for each domain are then coupled by means of exchange term accounting for the mass transfer of water and solutes between both regions (van Genuchten and Wierenga, 1976;Gaudet et al., 1977;Feyen et al., 1998).
Other approaches, which account for the presence of macropores that form a separate network, consider the soil as a collection of stream tubes (Toride and Leij, 1996).Within each stream tube, the water and solute transport equations are used.It is assumed that there is no exchange of water and solute between these tubes and that within each tube, the parameters determining water flow and solute transport are deterministic, but vary between tubes.It is worth noting that in this approach neither the spatial structure of the heterogeneity of soil properties nor hydraulic conductivity and dispersion coefficient are considered (Feyen et al., 1998).

Unsaturated Flow Models
Analytical solutions, when available, provide better insight into the physics behind the transport phenomena and are computationally efficient and simple to use.However, analytical approaches are for the most part limited to situations of simple geometry domains, linear governing equations, and homogeneous systems.Along with efficient numerical methods and rapidly updated computer hardware a large number of numerical models have been developed.However, the numerical technique cannot replace the analytical approaches completely, since numerical methods themselves need verification against analytical solutions because of discretization errors and convergence and stability problems that may be especially troublesome for advection-dominated and nonlinear adsorption processes.Some models are intense in computational demand when modeling a heterogeneous porous medium system with uncertainties in physical properties and distributions using stochastic approaches.
Analytical solutions to the Richards equation for unsaturated flow under various boundary and initial conditions are difficult to obtain because of the nonlinearity in soil hydraulic parameters.This difficulty is exaggerated in the case where soil is heterogeneous.Generally, one has to rely on numerical approaches for predicting moisture movement in unsaturated soils, even for homogeneous soils.However, numerical approaches often suffer from convergence and mass balance problems (Celia et al., 1990).Numerous analytical and numerical solutions of the unsaturated flow equation have been developed (Narasimhan et al., 1978;van der Heijde et al., 1985;Warrick et al., 1991;Parlange et al., 1997;Rockhold et al., 1997;Miller et al., 1998).Only a brief review of the modeling work is given here.
Exact analytical solutions typically require specialized forms of the hydraulic conductivity and diffusivity functions.The exponential hydraulic conductivity function (Gardner, 1958) has been widely used, but it is known to have a limited range of application to many real soils.Note that other functions, such as those developed by Brooks and Corey (1964) and van Genuchten (1980) are more firmly established for practical applications.Such special forms of the hydraulic functions make it possible to linearize the governing flow equations, and hence solve them analytically.Solutions to the linearized unsaturated flow equations are generally limited to steady state flow in semi-infinite, homogeneous soils (Broadbridge and White, 1988;Warrick and Yeh, 1990), and to transient flow in homogeneous and layered soils (Srivastava and Yeh, 1991).
Exact or quasi-analytical solutions have been derived for simplified unsaturated flow equations and mostly used for verification of numerical solutions (Philip, 1955 and1957;Warrick, 1974;Knight, 1983).The method of characteristics has also been applied to gravity-dominated flow in the unsaturated zone (Charbeneau, 1984).As it neglects the role of capillary pressure gradients, the hydraulic conductivity is set equal to the Darcy (volume) flux.This approach is also useful for estimating the unsaturated hydraulic conductivity.
Analytical solutions are only applicable to highly simplified systems, and are not well suited for situations normally encountered in the field.Initially, finite difference methods were mainly developed to predict unsaturated flow only (Brandt et al., 1971;Freeze, 1971); but finite element solution techniques also became available later (Huyakorn et al., 1984;Sudicky, 1989;Paniconi et al., 1991;Segol, 1993).The nonlinearity of Richards equation is usually solved using an iterative procedure such as Newton or Picard methods (Kirkland et al., 1992).Perhaps the most important advantage of finite element techniques over standard finite difference methods is the ability to describe irregular system boundaries in simulations more accurately, as well as easily including nonhomogeneous medium properties.For one-dimensional simulations, finite difference methods are just as good as finite element schemes (Moridis and Reddel, 1991).However, several authors suggest that finite element methods lead to more stable and accurate solutions, thus permitting larger timesteps and/or coarser grid systems, and hence leading to computationally more efficient numerical schemes.
One of the most successful finite difference schemes is the block centered approach, in which the spatial domain is discretized into blocks and the solutions are found at the center of these blocks.In two spatial dimensions, the scheme couples block center to only its 4 direct neighbours; therefore it is called a 5point scheme (7-point scheme in three spatial dimensions).The generated system is regularly structured (5 bands), symmetric and can be solved efficiently.This scheme was also shown to be equivalent to the lowest order mixed finite element method (Ewing, 1983).
To numerically solve coupled systems of equations, the solution process requires some manipulation at each timestep so that the dependence of one equation on the solution of the other is dealt with accurately.One way to overcome this is to use a fully implicit approach to solve the equations simultaneously.Any nonlinearity of the generated system can be handled by Newton's method.The implicit nature of this scheme allows for larger timesteps in the simulation to find stable solutions as compared to the timesteps for explicit schemes.However, this algorithm requires solving a system of equations of order equal to twice the number of grid points for each Newton iteration within each timestep.
An alternative to the fully implicit scheme is to use the mixed implicit-explicit scheme.However, the explicit part of the scheme means that this algorithm is now subject to a stability constraint which severely restricts the size of the timestep and introduces numerical artifacts.
Another approach is to decouple the two equations, and then solve both equations by implicit schemes, where different methods and space-time grids can be used on each one.The procedure is to solve the first equation using an extrapolated approximation of the other solution, and then solve the second equation.
Infiltration and flow through structured soils can be substantially different from the infiltration and flow in relatively homogenous materials (Beven and Germann, 1982) because the structured soils contain relatively large voids, such as inter-aggregate pores, drying cracks in clay soils, decayed root channels, and other type of "macropores".Attempts to describe flow in structured soils have generally centred on tworegion approaches.One region consists of the soil matrix in which the flow is described by the conventional, Darcian-based unsaturated flow equation, and the other region consists of either a single macropore or of a statistical network of macropores through which water flows primarily under the influence of gravity.The two regions are connected through some common boundary conditions, or by means of a simple source-sink term (Davidson, 1985).Other approaches include the use of kinematic wave theory (Charbeneau, 1984;Germann andBeven, 1985 and1986), where the self-sharpening wave (pulse) corresponds to an isolated wetting front drained to field capacity.
In order to apply the unsaturated flow model to a field situation, the hydraulic conductivity of the site has to be specified.However, the field hydraulic conductivity distribution commonly exhibits a high degree of spatial variability (Sudicky, 1986;Yeh, 1992).In addition, the spatial variability also varies with the scale of the problem (Dagan, 1986;Gelhar, 1986).The classical deterministic approach implies that hydraulic conductivity values in a model are known and specified at all points in the solution domain.To acquire such a detailed hydraulic conductivity distribution at the scale of tens of kilometers is generally considered impractical and unfeasible.The alternative is to utilize a small number of samples to estimate the variability of hydraulic conductivity in a statistical framework.That is, the spatial variability of the hydraulic conductivity is characterized by its probability distribution estimated from samples.However, analyses of field data show that although the conductivity values vary significantly in space, the variation is not entirely random, but correlated in space (Bakr et al., 1978;Yeh, 1992).Such a correlated nature implies that the values are not statistically independent in space and they must be treated as a stochastic process, instead of a single random variable (Zhang et al., 1998).Note that the stochastic approach does not provide information about the values of hydraulic conductivity at a specific location, but it provides a way to quantify its spatial variability.That is, it determines where the mean value lies and how widely the values spread around the mean value.However, in a numerical model it is required to estimate the values at certain points; in this case, the ordinary kriging is used to incorporate spatial correlation structure obtained from site measurements (de Marsily, 1986).
Many stochastic approaches are available including the geostatistics method, the effective parameter approach, the Monte Carlo simulation, and the conditional simulation (Yeh, 1992).Although each of these methods has limitations, one should bear in mind, as already mentioned above, that the uncertainty estimates afforded by the stochastic models are themselves uncertain (de Marsily, 1986).
For heterogeneous soils that contain macropores, a different modeling approach is needed, as the presence of macropores in these soils may form a separate network for water flow (Beven and Germann, 1982).The common approach is to introduce two-region flow domain, one for macropores and the other for the soil matrix.In each flow domain, hydraulic conductivity is prescribed separately (Feyen et al., 1998).
Another approach is to consider the heterogeneous soil as a collection of stream tubes (Toride and Leij, 1996).It is assumed that there is no exchange of water between these tubes, and that within each tube, the hydraulic conductivity is defined, but varies between tubes (Feyen et al., 1998).

Solute Transport Models
Transport of dissolved solutes in soils is commonly described by the advection-dispersion equation.Analytical solutions have been derived for a variety of boundary and initial conditions (Starr and Parlange, 1976;Barry and Sposito, 1989;Tang andAral, 1992a, 1992b;Segol, 1993;Leij andToride, 1995, 1998;Wu et al., 1997).Although these solutions are obtained for narrowly defined conditions, they have many applications such as the verification of computer codes, prediction of solute movement for large times or distances where the use of numerical models become impractical, and the determination of transport parameters from soil column studies.The majority of analytical solutions pertain to semi-infinite and infinite media.Solutions are obtained by a variety of mathematical techniques including Green's functions (Lindstorm and Boersma, 1989), separation of variables (Bruch and Street, 1967), characteristics method (Wilson and Gelhar, 1981;Charbeneau, 1984;Illangasekare and Doll, 1989), Laplace transforms (Leij et al., 1991) and Fourier transforms (Cleary and Adrian, 1973).
Prediction of solute migration under field conditions requires the simultaneous solution of the unsaturated flow and solute transport equations (Warrick et al., 1971;Bresler, 1973).Numerical solutions of the combined unsaturated flow and solute transport equations have been studied by Pickens and Gillham (1980);van Genuchten (1982); and Huyakorn et al. (1985) among others.First approximations involve or assume steady flow and constant water contents.
The common approach for modeling the sorption is to assume instantaneous adsorption or exchange reactions as well as simple linearity between s and c, the solute concentrations associated with the solid and the solution phases of the soil, respectively: where k is often referred to as the distribution coefficient .The source-sink term of equation ( 7) is often approximated by . Assuming steady state flow, the solute transport equation ( 7) simplifies to the classical linear convection-dispersion equation where the retardation factor R is given by , µ is first-order solute decay rate, and γ is the zero-order solute production rate.Application of this equation to transport through disturbed soil columns in the laboratory and in relatively uniform field soils involving nonreactive or only weakly reactive solutes has been fairly successful.
Chemically controlled kinetic rate reactions of the form ( ) have been examined.The simplest first-order linear kinetics is commonly assumed; and ignoring the degradation terms, it leads to the coupled system ( ) where α is a first-order rate coefficient.Success of this model has generally been limited (van Genuchten et al., 1974), thus suggesting that other than linear first-order kinetic processes dominate the sorption process during water movement.An alternative approach is to partition soil water into mobile (flowing) and stagnant (immobile or nonmoving) phases (Gaudet et al., 1977).That is, convective-dispersive transport is confined to only a fraction of the liquid-filled pores, and the remainder of the pores contain stagnant water.This stagnant water has been visualized as thin films around soil particles, as dead-end pores, or as relatively isolated regions associated with unsaturated flow.Van Genuchten and Wierenga (1986) extended the concepts of mobileimmobile water to include adsorption-desorption processes ( ) where the subscripts m and i refer to mobile and immobile regions, respectively.The success of this model for a variety of solutes in laboratory scale transport processes has been clearly demonstrated (Nkedi-Kizza et al., 1983).Adsorption-exchange inside aggregate is an important and often overlooked phenomenon in transport studies involving reactive tracers (Gvirtzmann et al., 1988).Sorption and exchange sites are often not located along the main fluid flow lines, but along dead-end pores, or inside clusters of particles that have larger internal surface areas.Even with a small amount of stagnant water in the system, an uneven distribution of sorption sites can lead to significant preferential transport.
Since the two-region modeling approach is conceptually pleasing and has resulted in improved prediction capabilities, it has then been extended to simulate transport in soils that contain well-defined macropores or that are made up of relatively large, uniformly sized and shaped aggregates.It is commonly assumed that the solute is transported through a single, well-defined pore or crack of known geometry, or through the voids between well-defined, uniformly sized aggregates.In addition, diffusion-type equations are used to describe the transport of solute from the larger pores into the micropores of the soil matrix (Tang et al., 1981;Sudicky and Frind, 1982;Goltz and Robert, 1986;Gerke and van Genuchten, 1993).
Because of the natural complexity of unsaturated flow, methods of predicting solute transport have relied largely on finite difference or finite element approximations of the governing equations (Galeati and Gambolati, 1989).Since advection-diffusion transport equations normally have no closed form analytical solutions, it is important to have accurate numerical approximations.When diffusion dominates the physical process, standard finite difference and element methods work well in solving these equations.However, when advection is the dominant process, these equations may present many numerical difficulties.In fact, any standard finite difference or element method will produce a solution, which exhibits non-physical oscillations (Johnson, 1987).A class of Eulerian methods has been developed to overcome these difficulties.The scheme uses Eulerian fixed grid and improved techniques, such as upstream weighting, to generate more accurate approximations (Celia et al., 1989).Although it is simple to formulate and to implement the scheme, these solutions suffer from excessive time truncation errors, as the scheme requires small timesteps to generate stable solutions (Eriksson and Johnson, 1993;Al-Lawatia et al., 1999).
Another class of characteristic methods can also overcome the effect of advection dominated problems.The method uses a combination of Eulerian fixed grids to treat the diffusive component and Lagrangian coordinates by tracking particles along the characteristics to treat the advective component.It has the advantage of greatly reducing the time truncation errors, thus allowing for large timestep, but the scheme has difficulty in conserving mass and treating general boundary conditions.
The Eulerian-Lagrangian Localized Adjoint Method (ELLAM), which was first proposed by Celia et al. (1990), provides a consistent framework for solving advection-dominated equations with general boundary conditions in a mass conservative manner.The ELLAM formulation was based on definition test functions that satisfy the homogeneous space-time adjoint equation locally, on elements that have boundaries defined by space-time characteristic curves of the equation.A second-order Runge-Kutta approximation of the characteristics within the framework of ELLAM has been developed by Al-Lawatia et al. (1999) for advection-diffusion equations and compared to other widely used and well known methods.ELLAM techniques have been extended to simulate contaminant transport in groundwater (Ewing and Wang, 1995).
One of the distinctive features of the porous media on the field scale is the spatial heterogeneity of transport properties (Vanderborght et al., 1997).These features have a distinct effect on the spatial distribution of contaminant concentration, as has been observed in field experiments (Roth et al., 1991;Ellsworth et al., 1991Ellsworth et al., , 1996) ) and demonstrated by simulation of contaminant transport in unsaturated, heterogeneous soil (Russo et al., 1998).Description of the mixing process due to spatial variability of the unsaturated hydraulic conductivity has been advanced with the development of numerical solutions, which assume spatially variable soil properties (Hills et al., 1991;Ellsworth and Boast, 1996); stochastic models (Bresler and Dagan, 1981;Bresler, 1987;Russo, 1993Russo, , 1995)); and stochastic stream tube models, which decompose the field into a set of independent vertical soil columns (Jury and Roth, 1990;Toride and Leij, 1996).Note that the predictions of these stochastic models on contaminant spreading are based on ensemble mean statistics, and hence, they may not be appropriate for particular site-specific applications.The stream tube models are distinguished by the degree of lateral mixing, which do not allow horizontal mixing; and the concentration for each tube represents a discrete value in the horizontal plane.This model is not popular due to the limited discussion of its theoretical foundation and a lack of procedures for its application (Toride and Leij, 1996).

Concluding Remarks
Groundwater is of crucial importance because it provides the source for most of the Sultanate's irrigation and domestic supplies.Every effort should be made to understand ways and means of efficiently using and managing it.The unsaturated zone is the region through which water, together with pollutant carried by the water, must pass to reach the groundwater.Therefore various processes occurring within the unsaturated zone play a major role in determining both the quality and quantity of water recharging into the groundwater.
In this review we have attempted to illustrate the advances made towards the progress in the understanding of the flow and contaminant transport processes in the unsaturated zone as our soil and groundwater resources are increasingly subjected to the dangers of long-term pollution.We concentrated initially on conceptual aspects of various deterministic-mathematical approaches for modeling flow and contaminant transport in the unsaturated zone, and also on various stochastic and statistical approaches geared toward field scale variability and heterogeneity.Unsaturated flow in the porous media is classically modeled using Richards equation, which is often solved using numerical approximation methods, such as finite difference or finite element methods.
Predicting water flow and contaminant transport on a field-scale based on the current monitoring and modeling techniques is a difficult task.There are large uncertainties in predictions mainly due to our inability to depict detailed spatial distributions of soil hydraulic properties on the field-scale.Due to the high costs of data acquisition, few field measurements are usually available for the characterization of flow and contaminant transport, even though the spatial distribution of a contaminant plume may be highly irregular.Thus the quantitative methods used for predicting the fate of a contaminant and its potential threat to groundwater quality must provide descriptions of transport based on statistical information, which can realistically be extracted from the data.
Classical, deterministic representations are generally not adequate to describe flow and transport in natural soils under field conditions, due to the occurrence of macropores and structured elements.Stochastic methods are simply mathematical tools used to overcome uncertainties in acquiring detailed spatial distributions of flow parameters in field scale aquifers.That is, instead of dealing with physical flow parameters, models dealt with the mean value, variance and correlation scales of the flow parameters.Quantitative descriptions of water flow and contaminant transport rely mainly upon model simulations in which, due to difficulty in obtaining its input parameters, hypothetical parameters are often used.Therefore, extensive validation of the models is required before the stochastic models can be widely accepted.
Oman has extensive areas of active and fossil karst.As aquifers are inherently heterogeneous at various observation scales, the complex geological structure of Oman reflects well-developed fracture systems in the rocks forming the northern mountain aquifers.The essence of fractured aquifers is the existence of an interconnected set of voids capable of conveying water and contaminants at speeds much higher than the typical Darcian flow rate.Therefore, more research associated with water flow and contaminant transport in the unsaturated zone of aquifers containing fractures and karstic conduits is needed for future investigations.

Acknowledgement
The authors wish to thank the referees for their comments which the authors found to be very valuable in revising the paper.Specially, to one referee for his detailed and thoughtful comments.