ENDOCRINE, PHYSIOLOGICAL AND HISTOPATHOLOGICAL RESPONSES OF FISH AND THEIR LARVAE TO STRESS WITH EMPHASIS ON EXPOSURE TO CRUDE Science and Technology, Special Review (2000) 265-280 ©2000 Sultan Qaboos University 265 Modeling Flow and Transport in Unsaturated Porous Media: A Review H.H. Al-Barwani, M. Al-Lawatia, E. Balakrishnan, and A. Purnama Department of Mathematics and Statistics, College of Science, P.O. Box 36, Sultan Qaboos University, Al-Khod 123, Muscat, Sultanate of Oman مراجعة وتحليل: ة نماذج حركة المياه و النقل في مناطق التهوية في األوساط المسامي حمدى البرواني، محمد اللواتيا، اسوارن بالكرشنان ، انتون برناما ذل كل جهد لفهم الطرق و الوسائل الممكنة لالستفادة يجب أن نب المياه الجوفية مورد طبيعي و حيوي و :خالصـة التي تحد من أعلى بالسطح اليابس و من اسفل هي تلك منطقة التهوية . فـيها بكفاءة عالية مـنها و التصـرف هي المنطقة التي من خاللها ترشح المياه مع المياه الملوثة لتصل إلى المياه الجوفية، و عليه وبالطاولـة المائـية تحديد كل من كمية و نوع المياه التي فالعملـيات المختلفة التي تحدث داخل منطقة التهوية تلعب دوراً أساسيا في ير غ بتدفق المياه و نقل المياه الملوثة في الوسط المسامي تنبؤللالطرق التقليدية . ترشـح داخـل المـياه الجوفية الطبيعة و ذلك نسبة لوجود مسامات كبيرة و عناصر على حاالت مشـبع بصورة عامة غير وافية عند تطبيقها ال مخططات النقل الملوث تحتاج إلى حل متزامن لمعادالت . عند المواقع المختلفة ص التربة منظمة و تغيرات لخوا قد ازداد كثيرا استعمال الحلول العددية و المحاكات في اآلونة األخيرة مشـبع ومعادالت النقل و الالـتدفق غـير جازه و آخر ما تم التوصل و في هذا البحث سوف نستعرض ما تم إن. الحاسـوبية المبنية علي المخططات العددية إلـيه فـي مخططـات التدفق المائي و النقل الملوث في منطقة التهوية و سوف نؤكد علي ما يجب التركيز عليه .مستقبال و خاصة ما له عالقة بسلطنة عمان ABSTRACT: 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 CONTENTS 1. Introduction 266 2. Mathematical Equations of Water and Transport In Unsaturated Soils 267 3. Unsaturated Flow Models 270 4. Solute Transport Models 272 5. Concluding Remarks 275 6. Acknowledgment 275 7. References 275 1. Introduction Oman 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. 266 AL-BARWANI, AL-LAWATIA, BALAKRISHNAN, and PURNAMA 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. 2. 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, z Kq ∂ Ψ∂ −= , (1) 267 MODELING FLOW AND TRANSPORT IN UNSATURATED POROUS MEDIA 268 K zh + 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 G z q t = ∂ ∂ + ∂ ∂θ , (2) leads to Richards equation for water flow in the unsaturated zone. In the mixed form, the unsaturated flow equation is written as G z K z h K zt + ∂ ∂ +      ∂ ∂ ∂ ∂ = ∂ ∂θ , (3) 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 G z K z h K zt h C + ∂ ∂ +      ∂ ∂ ∂ ∂ = ∂ ∂ , (4) where the specific moisture capacity C is the slope of the soil water retention curve ( )hθ ; and the moisture content θ -based equation G z K z D zt + ∂ ∂ +      ∂ ∂ ∂ ∂ = ∂ ∂ θθ , (5) where it is customary to define the coefficient of moisture diffusivity as CK=D . 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 and 1989b). 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 AL-BARWANI, AL-LAWATIA, BALAKRISHNAN, and PURNAMA 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) ( )mnrsr hαθθθθ +−+= 1/)( , (6) where rθ and are residual and field-saturated volumetric water contents, respectively, and , n, and m are empirical constants. For the hydraulic conductivity functions, Mualem’s (1976) formulation has more potential than other models as it is more applicable to a wider variety of soils (van Genuchten and Nielsen, 1985). αsθ 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 ( ) ( ) Pcq z c D zt c t s +      − ∂ ∂ ∂ ∂ = ∂ ∂ + ∂ ∂ θ θρ , (7) 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). 269 MODELING FLOW AND TRANSPORT IN UNSATURATED POROUS MEDIA 270 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). 3. 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). AL-BARWANI, AL-LAWATIA, BALAKRISHNAN, and PURNAMA 271 Exact or quasi-analytical solutions have been derived for simplified unsaturated flow equations and mostly used for verification of numerical solutions (Philip, 1955 and 1957; 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 5- point 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 two- region 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 and Beven, 1985 and 1986), where the self-sharpening wave (pulse) corresponds to an isolated wetting front drained to field capacity. MODELING FLOW AND TRANSPORT IN UNSATURATED POROUS MEDIA 272 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). 4. 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 and Aral, 1992a, 1992b; Segol, 1993; Leij and Toride, 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. AL-BARWANI, AL-LAWATIA, BALAKRISHNAN, and PURNAMA 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: cks = , (8) where k is often referred to as the distribution coefficient . The source-sink term of equation (7) is often approximated by dK γµ +−= cP . Assuming steady state flow, the solute transport equation (7) simplifies to the classical linear convection-dispersion equation γµ +− ∂ ∂ − ∂ ∂ = ∂ ∂ c z c v z c D t c R 2 2 , (9) where the retardation factor R is given by θρ kR += 1 , µ 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 ( )csfts ,=∂∂ have been examined. The simplest first-order linear kinetics is commonly assumed; and ignoring the degradation terms, it leads to the coupled system ( )sck t s z c v z c D t c t s −= ∂ ∂ ∂ ∂ − ∂ ∂ = ∂ ∂ + ∂ ∂ α θ ρ ; 2 2 , (10) 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 mobile- immobile water to include adsorption-desorption processes ( )imiiimmmmmmiiimmm cct c R z c v z c D t c R t c R −= ∂ ∂ ∂ ∂ − ∂ ∂ = ∂ ∂ + ∂ ∂ αθθθθθ ; 2 2 (11) 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. 273 MODELING FLOW AND TRANSPORT IN UNSATURATED POROUS MEDIA 274 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., 1991, 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, 1993, 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). AL-BARWANI, AL-LAWATIA, BALAKRISHNAN, and PURNAMA 275 5. 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. 6. 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. 7. References AL-LAWATIA, M., SHARPLEY, R.C. and WANG, H. 1999. Second-Order Characteristic Methods for Advection- Diffusion Equations and Comparison to Other Schemes. Adv. Wat. Resour., 22: 741-768. MODELING FLOW AND TRANSPORT IN UNSATURATED POROUS MEDIA 276 ARYA, L.M. and PARIS, J.F. 1981. A Physicoempirical Model to Predict the Soil Moisture Characteristic from Particle-Size Distribution and Bulk Density Data. Soil Sci. Soc. Am. J., 45: 1023-1030. ASSOULINE, S., TESSIER, D. and BRUAND, A. 1998. A Conceptual Model of the Soil Water Retention Curve. Wat. Resour. Res., 34: 223-231. BAKR, A.A., GELHAR, L.W., GUTJAHR, A.L. and MACMILLAN, J.R. 1978. Stochastic Analysis of Spatial Variability in Subsurface Flows. 1. Comparison of One-Dimensional Flows. Wat. Resour. Res., 14: 263-271. BARRY, D.A. and SPOSITO, G. 1989. Analytical Solution of a Convection-Dispersion Model with Time-Dependent Transport Coefficients. Wat. Resour. Res., 25: 2407-2416. BEAR, J. 1972. Dynamics of Fluids in Porous Media. Elsevier, New York. BEVEN, K. and GERMANN, P. 1982. Macropores and Water Flow in Soils. Wat. Resour. Res., 18: 1311-1325. BIRD, N.R., BARTOLI, A.F. and DEXTER, A.R. 1996. Water Retention Models for Fractal Soil Structures. Eur. J. Soil Sci., 47: 1-6. BRANDT, A., BRESLER, E., DINER, N., BEN-ASHER, I., HELLER, J. and GOLDBERG, G. 1971. Infiltration from a Trickle Source. Soil Sci. Soc. Am. J., 35: 675-682. BRESLER, E. 1973. Simultaneous Transport of Solutes and Water Under Transient Unsaturated Flow Conditions. Wat. Resour. Res., 9: 975-986. BRESLER, E. 1987. Modeling of Flow, Transport, and Crop Yield in Spatially Variable Fields. Advances in Soil Science (edited by Stewart, B.A.), 7: 1-51. BRESLER, E. and DAGAN, G. 1981. Convective and Pore Scale Dispersive Solute Transport in Unsaturated Heterogeneous Fields. Wat. Resour. Res., 17: 1683-1689. BROADBRIDGE, P. and WHITE, I. 1988. Constant Rate Rainfall Infiltration: A Versatile Nonlinear Model. 1. Analytical Solution. Wat. Resour. Res., 24: 145-154. BROOKS, R.H. and COREY, A.T. 1964. Hydraulic Properties of Porous Media. Hydrology Paper 3, Colorado State Univ., Fort Collins, CO. BRUCH, J.C. and STREET, R.L. 1967. Two-Dimensional Dispersion. J. Sanit. Engng. Div., Proc. ASCE, 93: 17-39. CELIA, M.A., HERRERA, I., BOULOUTAS, E. and KINDRED, J.S. 1989. A New Numerical Approach for the Advective-Diffusive Transport Equation. Numer. Methods Partial Differential Equations, 5: 203-226. CELIA, M.A., BOULOUTAS, E. and ZARBA, R. 1990. A General Mass-Conservative Numerical Solution of the Unsaturated Flow Equation. Wat. Resour. Res., 26: 1483-1496. CHARBENEAU, R.J. 1984. Kinematic Models for Soil Moisture and Solute Transport. Wat. Resour. Res., 20: 699-706. CHILDS, E.C. and COLLIS-GEORGE, N. 1950. The Permeability of Porous Materials. Proc. Roy. Soc. Lond., A201: 392-405. CLEARY, R.W. and ADRIAN, D.D. 1973. Analytical Solution of the Convective-Dispersive Equation for Cation Adsorption in Soils. Soil Sci. Soc. Am. J., 37: 197-199. DAGAN, G. 1986. Statistical Theory of Groundwater Flow and Transport: Pore to Laboratory, Laboratory to Formation and Formation to Regional Scale. Wat. Resour. Res., 22: 120s-134s. DAGAN, G. and BRESLER, E. 1979. Solute Transport in Unsaturated Heterogeneous Soil at Field Scale. 1. Theory. Soil Sci. Soc. Am. J., 43: 461-467. DAVIDSON, M.R. 1985. A Numerical Study of Saturated-Unsaturated Infiltration in a Cracked Soil. Wat. Resour. Res., 21: 709-714. DE MARSILY, G. 1986. Quantitative Hydrogeology: Groundwater Hydrogeology for Engineers. Academic Press, Orlando, Florida. ECHING, S.O., HOPMANS, J.W. and WENDROTH, O. 1994. Unsaturated Hydraulic Conductivity from Transient Multistep Outflow and Soil Water Pressure Data. Soil Sci. Soc. Am. J., 58: 687-695. ELLSWORTH, T.R., JURY, W.A., ERNST, F.F. and SHOUSE, P.J. 1991. A Three-Dimensional Field studies of Solute Tranport through Unstructured Layered Porous Media. Methodology and, Mass Recovery and, Mean Transport. Wat, Resour. Res., 27: 951-965. ELLSWORTH, T.R. and BOAST, C.W. 1996. Spatial Structure of Solute Transport Variability in an Unsaturated Field Soil. Soil Sci. Soc. Am. J., 60: 1355-1367. ELLSWORTH, T.R., SHOUSE, P.J., SKAGGS, T.H., JOBES, J.A. and FARGERLUND, J. 1996. Solute Transport in Unsaturated Soil: Experimental Design, Parameter Estimation, and Model Discrimination. Soil Sci. Soc. Am. J., 60: 397-407. ERIKSSON, K. and JOHNSON, C. 1993. Adaptive Streamline Diffusion Finite Element Methods for Stationary Convection Diffusion Problems. Math. Comp., 60: 167-188. AL-BARWANI, AL-LAWATIA, BALAKRISHNAN, and PURNAMA 277 EWING, R.E. (Editor). 1983. The Mathematics of Reservoir Simulation Frontiers in Applied Mathematics. Vol. 1. SIAM, Philadelphia. EWING, R.E and WANG, H. 1995. Optimal-Order Convergence Rate for Eulerian-Lagrangian Localized Adjoint Method for Reactive Transport and Contamination in Groundwater. Numer. Methods Partial Differential Equations, 11: 1-31. FEYEN, J., JACQUES, D., TIMMERMAN, A. and VANDERBORGHT, J. 1998. Modelling Water Flow and Solute Transport in Heterogeneous Soils: A Review of Recent Approaches. J. Agric. Engng. Res., 70: 231-256. FLURY, M. 1996. Experimental Evidence of Transport of Pesticides Through Field Soils - A Review. J. Environ. Qual., 25: 25-45. FREEZE, R.A. 1971. Three-Dimensional, Transient, Saturated-Unsaturated flow in a Groundwater Basin. Wat. Resour. Res., 7: 347-366. FREEZE, R.A. and CHEERY, J.A. 1979. Groundwater. Prentice-Hall, Englewood Cliffs, N.J. GALEATI, G. and GAMBOLATI, G. 1989. On Boundary Conditions and Point Sources in the Finite Element Integration of the Transport Equation. Wat. Resour. Res., 25: 847-856. GARDNER, W.R. 1958. Some Steady State Solutions of Unsaturated Moisture Flow Equations with Application to Evaporation from a Water Table. Soil Sci., 85: 228-232. GAUDET, J.P., JEGAT, H., VACHAUD, G. and WIERENGA, P.J. 1977. Solute Transfer, with Exchange between Mobile and Stagnant Water, through Unsaturated Sand. Soil Sci. Soc. Am. J., 41: 665-671. GELHAR, L.W. 1986. Stochastic Subsurface Hydrology from Theory to Applications. Trends and Directions in Hydrology (edited by Burges, S.J.), Wat. Resour. Res., 22: 135S-145S. GERKE, H.H. and VAN GENUCHTEN, M.Th. 1993. A Dual-Porosity Model for Simulating the Preferential Movement of Water and Solutes in Structured Porous Media. Wat. Resour. Res., 29: 305-319. GERMANN, P. and BEVEN, K. 1985. Kinematic Wave Approximation to Infiltration into Soils with Sorbing Macropores. Wat. Resour. Res., 21: 990-996. GERMANN, P. and BEVEN, K. 1986. A Distribution Function Approach to water Flow in Soil Macropores Based on Kinematic Wave Theory. J. Hydrol., 83: 173-183. GHODRATI, M. and JURY, W.A. 1990. A Field Study Using Dyes to Characterize Preferential Flow of Water. Soil Sci. Soc. Am. J., 54: 1558-1563. GOLTZ, M.N. and ROBERT, P.V. 1986. Interpreting Organic Solute Transport data From a Field Experiment Using Physical Nonequilibrium Models. J. Contam. Hydrol., 1: 77-93. GRAY, W.G. and HASSANIZADEH, S.M. 1991. Paradoxes and Realities in Unsaturated Flow Theory. Wat. Resour. Res., 27: 1847-1854. GVIRTZMANN, H., PALDOR, N., MAGARITZ, M. and BACHMAT, Y. 1988. Mass Exchange Between Mobile Freshwater and Immobile Saline Water in the Unsaturated Zone. Wat. Resour. Res., 24: 1638-1644. HILLS, R.G., PORRO, I., HUDSON, D.B. and WIERENGA, P.J. 1989a. Modeling One-Dimensional Infiltration Into Very Dry Soils. 1. Model Development and Evaluation. Wat. Resour. Res., 25: 1259-1269. HILLS, R.G., HUDSON, D.B., PORRO, I. and WIERENGA, P.J. 1989b. Modeling One-Dimensional Infiltration Into Very Dry Soils. 2. Estimation of the Soil Water Parameters and Model Predictions. Wat. Resour. Res., 25: 1271- 1282. HILLS, R.G., WIERENGA, P.J., HUDSON, D.B. and KIRKLAND, M.R. 1991. The Second Las Cruces Trench Experiment: Experimental Results and Two-Dimensional Flow Predictions. Wat. Resour. Res., 27: 2707-2718. HUYAKORN, P.S., THOMAS, S.D. and THOMPSON, B.M. 1984. Techniques for Making Finite Elements Competitive in Modeling Flow in Variably Saturated Porous Media. Wat. Resour. Res., 20: 1099-1115. HUYAKORN, P.S., MERCER, J.W. and WARD, D.S. 1985. Finite Element Matrix and Mass Balance Computational Schemes for Transport in variably Saturated Porous Media. Wat. Resour. Res., 21: 346-358. ILLANGASEKARE, T.H. and DOLL, P. 1989. A Discrete Kernel Method of Characteristics Model of Solute Transport in Water Table Aquifers. Wat. Resour. Res., 25: 857-867. JOHNSON, C. 1987. Numerical Solution of Partial Differential Equations by the Finite Element Method. Cambridge University Press, Cambridge. JU, S.-H. and KUNG, K.-J.S. 1997. Impact of Funnel Flow on Contaminant Transport in Sandy Soils: Numerical Simulation. Soil Sci. Soc. Am. J., 61: 409-415. JURY, W.A. and GRUBER, J. 1989. A Stochastic Analysis of the Soil and Climate Variability on the Estimation of Pesticide Groundwater Pollution Potential. Wat. Resour. Res., 25: 2465-2474. JURY, W.A. and K. ROTH. 1990. Transfer Functions and Solute Movement Through Soil: Theory and Application. Birkhauser Verlag, Basel, Switzerland. MODELING FLOW AND TRANSPORT IN UNSATURATED POROUS MEDIA 278 KIRKLAND, M.R., HILLS, R.G. and WIERENGA, P.J. 1992. Algorithms for Solving Richards Equation for variably Saturated Soils. Wat. Resour. Res., 28: 2049-2058. KNIGHT, J.H. 1983. Infiltration Equations from Exact and Appropriate Solutions of Richards Equation. In Advances in Infiltration, ASAE Publ. 11-83, New York. KOOL, J.B., PARKER, J.C. and VAN GENUCHTEN, M.Th. 1985. Determining Soil Hydraulic Properties from One- Step Outflow Experiments by Parameter Estimation. I. Theory and Numerical Studies. Soil Sci. Soc. Am. J., 49: 1348-1354. KOOL, J.B., PARKER, J.C. and VAN GENUCHTEN, M.Th. 1987. Parameter Estimation for Unsaturated Flow and Transport Models-A Review. J. Hydrol., 91: 255-293. LEIJ, F.J., SKAGGS, T.H. and VAN GENUCHTEN, M.Th. 1991. Analytical Solutions for Solute Transport in Three- Dimensional Semi-Infinite Porous Media. Wat. Resour. Res., 27: 2719-2733. LEIJ, F.J. and TORIDE, N. 1995. Discrete Time- and Length-Averaged Solutions of the Advection-Dispersion Equation. Wat. Resour. Res., 31: 1713-1724. LEIJ, F.J. and TORIDE, N. 1998. Analytical Solutions for Solute Transport in Finite Soil Columns with Arbitrary Initial Distributions. Soil Sci. Soc. Am. J., 62: 855-864. LINDSTORM, F.T. and BOERSMA, L. 1989. Theory of Chemical Transport with Simultaneous Sorption in a Water Saturated Porous Medium. Soil Sci., 110: 1-9. MILLER, C.T., WILLIAMS, G.A., KELLEY, C.T. and TOCCI, M.D. 1998. Robust Solution of Richards Equation for Nonuniform Porous Media. Wat. Resour. Res., 34: 2599-2610. MORIDIS, G.J. and REDDEL, D.L. 1991. The Laplace Transform Finite Difference Method for Simulation of Flow Through Porous Media. Wat. Resour. Res., 27: 1873-1884. MUALEM, Y. 1976. A New Model for Predicting the Hydraulic Conductivity of Unsaturated Porous Media. Wat. Resour. Res., 12: 513-522. NARASIMHAN, T.N., WITHERSPOON, P.A. and EDWARDS, A.L. 1978. Numerical Model for Saturated- Unsaturated Flow in Deformable Porous Media. 2. The Algorithm. Wat. Resour. Res., 14: 255-261. NIELSEN, D.R., VAN GENUCHTEN, M.Th. and BIGGAR, J.W. 1986. Water Flow and Solute Transport Processes in the Unsaturated Zone. Trends and Directions in Hydrology (edited by Burges, S.J.), Wat. Resour. Res., 22: 89S- 108S. NKEDI-KIZZA, P., BIGGAR, J.W., VAN GENUCHTEN, M. Th., WIERENGA, P.J., SELIM, H.M., DAVIDSON, J.M. and NIELSEN, D.R. 1983. Modeling Tritium and Chloride 36 Transport Through an Aggregated Oxisol. Wat. Resour. Res., 19: 691-700. PANICONI, C., ALDAMA, A.A. and WOOD, E.F. 1991. Numerical Evaluation of Iterative and Noniterative Methods for the Solution of the Nonlinear Richards Equation. Wat. Resour. Res., 27: 1147-1163. PARLANGE, J.-Y. 1980. Water Transport in Soils. Ann. Rev. Fluid Mech., 12: 77-102. PARLANGE, J.-Y., BARRY, D.A., PARLANGE, M.B., HOGARTH, W.L., HAVERKAMP, R., ROSS, P.J., LING, L. and STEENHUIS, T.S. 1997. New Approximate Analytical Technique to Solve Richards Equation for Arbitrary Surface Boundary Conditions. Wat. Resour. Res., 33: 903-906. PHILIP, J.R. 1955. Numerical Solution of Equations of the Diffusion Type with Diffusivity Concentration-Dependent. Trans. Faraday Soc., 51: 885-892. PHILIP, J.R. 1957. The Theory of Infiltration. 1. The Infiltration Equation and Its Solution. Soil Sci., 83: 345-357. PINDER, G.F. and GRAY, W.G. 1979. Finite Element Simulation in Surface and Subsurface Hydrology. Academic Press, New York. PICKENS, J.F. and GILLHAM, R.W. 1980. Finite Element Analysis of Solute Transport Under Hysteretic Unsaturated Flow Conditions. Wat. Resour. Res., 16: 1070-1078. RICHARDS, T.L. and STEENHUIS, T.S. 1988. Tile Drain Sampling of Preferential Flow on a Field Scale. J. Contam. Hydrol., 3: 307-325. ROCKHOLD, M.L., SIMMONS, C.S. and FAYER, M.J. 1997. An Analytical Solution Technique for One- Dimensional, Steady Vertical Water Flow in Layered Soils. Wat. Resour. Res., 33: 897-902. ROTH, K., JURY, W.A., FLUHLER, H. and ATTINGER, W. 1991. Transport of Chloride Through an Unsaturated Field Soil. Wat. Resour. Res., 27: 2533-2541. RUSSO, D. 1993. Stochastic Modeling of Macrodispersion for Solute Transport in Heterogeneous Unsaturated Porous Formation. Wat. Resour. Res., 29: 383-397. RUSSO, D. 1995. On the Velocity Covariance and Transport Modeling in Heterogeneous Anisotropic Porous Formation. 2. Unsaturated Flow. Wat. Resour. Res., 31: 139-145. AL-BARWANI, AL-LAWATIA, BALAKRISHNAN, and PURNAMA 279 RUSSO, D., ZAIDEL, J. and LAUFER, A. 1998. Numerical Analysis of Flow and Transport in a Three-Dimensional Partially Saturated Heterogeneous Soil. Wat. Resour. Res., 34: 1451-1468. SEGOL, G. 1993. Classical Groundwater Simulations: Proving and Improving Numerical Models. Prentice-Hall, Englewood Cliffs, N.J. SRIVASTAVA, R. and YEH, J.T.-C. 1991. Analytical Solutions for One-Dimensional, Transient Infiltration Toward the Water Table in Homogeneous and Layered Soils. Wat. Resour. Res., 27: 753-762. STARR, J.L. and PARLANGE, J.-Y. 1976. Solute Dispersion in Saturated Soil Columns. Soil Sci., 121: 364-372. SUDICKY, E.A. 1986. A Natural Gradient Experiment on Solute Transport in a Sand Aquifer: Spatial Variability of Hydraulic Conductivity and its Role in the Dispersion Process. Wat. Resour. Res., 22: 2069-2082. SUDICKY, E.A. 1989. The Laplace Transform Galerkin Technique: A Time-Continuous Finite Element Theory and Application to Mass Transport in Groundwater. Wat. Resour. Res., 25: 1833-1846. SUDICKY, E.A. and FRIND, E.O. 1982. Contaminant Transport in Fractured Porous Media: Analytical Solutions for a System of Parallel Fractures. Wat. Resour. Res., 18: 1634-1642. TANG, D.H., FRIND, E.O. and SUDICKY, E.A. 1981. Contaminant Transport in Fractured Porous Media: Analytical Solution for a Single Fracture. Wat. Resour. Res., 17: 555-564. TANG, Y. and ARAL, M.M. 1992a. Contaminant Transport in Layered Porous Media. 1. General Solution. Wat. Resour. Res., 28: 1389-1397. TANG, Y. and ARAL, M.M. 1992b. Contaminant Transport in Layered Porous Media. 2. Applications. Wat. Resour. Res., 28: 1399-1406. TORIDE, N. and LEIJ, F.J. 1996. Convective-Dispersive Stream Tube Model for Field-scale Solute Transport: 1. Moment Analysis. Soil Sci. Soc. Am. J., 60: 342-252. VANDERBORGHT, J., GONZALEZ, C., VANCLOOSTER, M., MALLANTS, D. and FEYEN, J. 1997. Effects of Soil Type and water Flux on Solute Transport. Soil Sci. Soc. Am. J., 61: 372-389. VANDERBORGHT, J., VANCLOOSTER, M., MALLANTS, D., DIELS, J. and FEYEN, J. 1996. Determining Convective Lognormal Solute Transport Parameters from Resident Concentration Data. Soil Sci. Soc. Am. J., 60: pp. 1306-1317. VAN DER HEIJDE, P., BACHMAT, Y., BREDEHOEFT, J., ANDREWS, B., HOLT, D. and SEBASTIAN, S. (Editors). 1985. Groundwater Management: The Use of Numerical Models. Geophys. Monogr. Ser. Vol. 5. AGU, Washington, D.C. VAN GENUCHTEN, M.Th. 1980. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Sci. Soc. Am. J., 44: 892-898. VAN GENUCHTEN, M.Th. 1982. A Comparison of Numerical Solutions of the One-Dimensional Unsaturated- Saturated Flow and Mass Transport Equations. Adv. Wat. Resour., 5: 47-55. VAN GENUCHTEN, M.Th., DAVIDSON, J.M. and WIERENGA, P.J. 1974. An Evaluation of Kinetic and Equilibrium Equations for the Prediction of Pesticide Movement in Porous Media. Soil Sci. Soc. Am. J., 38: 29- 35. VAN GENUCHTEN, M.Th. and NIELSEN, D.R. 1985. On Describing and Predicting the Hydraulic Properties of unsaturated soils. Annales Geophysicae, 3: 615-627. VAN GENUCHTEN, M.Th. and WIERENGA, P.J. 1976. Mass Transfer Studies in Sorbing Porous Media. 1. Analytical Solutions. Soil Sci. Soc. Am. J., 40: 473-480. VAN GENUCHTEN, M.Th. and WIERENGA, P.J. 1986. Mass Transfer Studies in Sorbing Porous Media. Soil Sci. Soc. Am. J., 40: 473-480. WAGNER, B.J. and GORELICK, S.M. 1986. A Statistical Methodology for Estimating Transport Parameters: Theory and Applications to One-Dimensional Advective-Dispersive Systems. Wat. Resour. Res., 22: 1303-1315. WARRICK, A.W. 1974. Simultaneous Solute and Water Transfer for an Unsaturated Soil. Wat. Resour. Res., 7: 1216- 1225. WARRICK, A.W., BIGGAR, J.W. and NIELSEN, D.R. 1971. Simultaneous Solute and Water Transfer for an Unsaturated Soil. Wat. Resour. Res., 7: 1216-1225. WARRICK, A.W., ISLAS, A. and LOMEN, D.O. 1991. An Analytical Solution to Richards Equation for Time-Varying Infiltration. Wat. Resour. Res., 27: 763-766. WARRICK, A.W. and YEH, T.-C.J. 1990. One-Dimensional, Steady Vertical Flow in a Layered Soil Profile. Adv. Wat. Resour., 13: 207-210. WILDENSCHILD, D., HOGH JENSEN, K., HOLLENBECK, K.J., ILLANGASEKARE, T.H., ZNIDARCIC, D., SONNENBORG, T. and BUTTS, M.B. 1997. A Two-stage Procedure for Determining Unsaturated Hydraulic Characteristics using a Syringe Pump and Outflow Observations. Soil Sci. Soc. Am. J., 61: 347-359. MODELING FLOW AND TRANSPORT IN UNSATURATED POROUS MEDIA WILSON, J.L. and GELHAR, L.W. 1981. Analysis of Longitudinal Dispersion in Unsaturated Flow. 1. The Analytical Method. Wat. Resour. Res., 17: 122-130. WU, Y.S., KOOL, J.B., HUYAKORN, P.S. and SALEEM, Z.A. 1997. An Analytical Model for Nonlinear Adsorptive Transport Through Layered Soils. Wat. Resour. Res., 33: 21-29. YEH, T. –C. J. 1992. Stochastic Modelling of Groundwater Flow and Solute Transport in Aquifers. Hydrol. Processes, 6: 369-395. ZHANG, D., WALLSTROM, T.C. and WINTER, C.L. 1998. Stochastic Analysis of Steady-State Unsaturated Flow in Heterogeneous Media: Comparison of the Brooks-Corey and Gardner-Russo Models. Wat. Resour. Res., 34: 1437-1449. Received 30 January 2000 Accepted 24 June 2000 280 Modeling Flow and Transport in Unsaturated Porous Media: A Review H.H. Al-Barwani, M. Al-Lawatia, E. Balakrishnan, and A. Purnama Sultan Qaboos University, Al-Khod 123, Muscat, Sultanate of Oman Introduction Mathematical Equations of Water and Transport in Unsaturated Soils Unsaturated Flow Models Solute Transport Models Concluding Remarks