Published online 27 May 2008
Published in Vadose Zone J 7:757-768 (2008)
DOI: 10.2136/vzj2007.0082
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: VADOSE ZONE MODELING
Evaluating Interactions between Groundwater and Vadose Zone Using the HYDRUS-Based Flow Package for MODFLOW
Navin Kumar C. Twarakavia,*,
Jirka 
im
neka and
Sophia Seob
a Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92521
b Colorado School of Mines, 1500 Illinois St., Golden, CO 80401
* Corresponding author (navink{at}ucr.edu).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 25 April 2007.
 |
ABSTRACT
|
|---|
In the past, vadose zone processes in groundwater flow models were dramatically simplified (or even neglected) due to constraints on computational resources. The one-dimensional unsaturated flow package HYDRUS, recently developed for the groundwater model MODFLOW, was evaluated and compared with other contemporary modeling approaches used to characterize vadose zone effects in groundwater models. Being fully incorporated into the MODFLOW program, the HYDRUS package provides MODFLOW with recharge fluxes at the water table, while MODFLOW provides HYDRUS with the position of the groundwater table that is used as the bottom boundary condition in the package. The performance of the HYDRUS package was analyzed for three case studies of increasing complexity: (i) a one-dimensional infiltration experiment; (ii) a two-dimensional water table recharge experiment; and (iii) a hypothetical regional-scale groundwater flow problem. The computational need and modeling efficiency of the HYDRUS package was compared with other relevant MODFLOW packages (VSF, UZF1, and REC-ET). For smaller scale problems (up to two dimensions), the VSF process and the HYDRUS and UZF1 packages performed comparably well in terms of modeling efficiency and simulation times. Because of the high computational demand, it was not feasible to use the VSF process on a typical personal computer for the hypothetical large-scale groundwater problem. The HYDRUS package provided a much more efficient alternative to VSF for this large-scale groundwater problem and could better account for vadose zone processes than the UZF1 and REC-ET packages. For large-scale groundwater problems, the HYDRUS package provides an optimal tradeoff between computational effort and accuracy of model simulations for coupled vadose zone–groundwater problems.
Abbreviations: ET, evapotranspiration • REC-ET, Recharge–Evapotranspiration package • UZF1, Unsaturated Zone Flow package • VSF, Variably Saturated Flow package
|
INTRODUCTION
|
|---|
WATER FLOW through the variably saturated (vadose) zone is an important part of the hydrologic cycle because it influences partitioning of water among various flow components. Depending on hydrologic, geologic, and soil characteristics, rain and snowmelt are partitioned at the land surface into runoff, infiltration, evapotranspiration (ET), groundwater recharge, and vadose zone storage (Fig. 1
). Water flow in the vadose zone especially affects the transfer rates between the land surface and the groundwater table, which are two key hydrologic boundaries. Evaluation of almost any hydrologic process, therefore, requires that water flow through the vadose zone is appropriately taken into account. Modeling of vadose zone flow processes, however, is a complex and computationally demanding task that is often handicapped by the lack of data necessary to characterize the hydraulic properties of the subsurface environment. Consequently, vadose zone flow processes have rarely been properly represented in hydrologic models (Ward, 2002; Scanlon, 2002; Keese et al., 2005). For example, models that simulate surface and near-surface hydrology usually oversimplify the impact of vadose zone flow processes and rarely consider three-dimensional regional groundwater flow. Similarly, regional-scale groundwater models often simplify vadose zone flow processes by calculating groundwater recharge externally without proper consideration of changes in groundwater levels (e.g., Lorenz and Delin, 2007; Shah et al., 2007; Uddameri and Kuchanur, 2007). To overcome this frequent simplification, there is an urgent need for methods that can effectively simulate water flow through the vadose zone in large-scale hydrologic models (Winter et al., 1998). This issue is especially important for groundwater models.
View larger version (30K):
[in this window]
[in a new window]
|
FIG. 1. A schematic showing the processes (including the key vadose zone processes) affecting subsurface hydrology.
|
|
Most traditional attempts at characterizing vadose zone flow processes in groundwater models follow the water budget or "residual" approach (Scanlon, 2002). Water budget methods are based on the water budget equation. The water budget equation is a mathematical representation of the fact that all water arriving at the water table leaves the system as groundwater flow, is discharged through sinks such as surface water bodies, is evapotranspirated, or is retained as storage (Scanlon, 2002). One may indirectly estimate the water table recharge from the water budget equation by measuring or estimating all other components in the water budget. An example of the water budget approach is the use of the Recharge and Evapotranspiration (REC-ET) packages in MODFLOW (Harbaugh et al., 2000), a modular three-dimensional finite-difference groundwater flow model. MODFLOW was developed by the U.S. Geological Survey and is one of the most widely used groundwater flow models. It is obvious that the REC-ET package tends to oversimplify and underestimate the effects of vadose zone flow on groundwater. In spite of its widespread use, the REC-ET approach suffers from the following major limitations: (i) the method is not reliable for deeper water tables; and (ii) the applicability of this approach is questionable for arid and semiarid regions where soil capillary pressures play a dominant role in vadose zone flow (Gee and Hillel, 1988; Lerner et al., 1990; Hendrickx and Walker, 1997).
A more promising approach to properly represent vadose zone flow processes in groundwater models involves coupling groundwater and vadose zone models. A coupled model simulates the effects of near-surface hydrologic processes on groundwater flow by linking a groundwater model with a selected vadose zone model in space and time. The majority of currently available vadose zone models are based on either the Richards equation (Richards, 1931) or the kinematic wave equation (Colbeck, 1972; Smith, 1983; Smith and Hebbert, 1983). While the Richards equation considers flow due to both capillary and gravity forces, the kinematic wave equation neglects capillarity and considers only gravity. In coupled models, the groundwater recharge is calculated internally in the model based on existing surface hydrologic conditions and water table levels. By simultaneously considering surface meteorological conditions, water table levels, and the hydraulic properties of the vadose zone, coupled models represent reality better than traditional approaches such as the REC-ET package; however, evaluation of interactions between the near-surface and groundwater flow processes using coupled models has been a desirable but difficult goal.
It is desirable to develop MODFLOW packages other than REC-ET that would better account for processes in the vadose zone. By combining these packages with MODFLOW, not only can the vast groundwater modeling capabilities of MODFLOW be harnessed, but these new numerical packages can be quickly distributed among the large number of MODFLOW users. The virtues of a coupled vadose zone–groundwater model should be evaluated based on the following criteria: (i) accuracy in representation of the physical processes that drive vadose zone flow, (ii) usability for different groundwater modeling scenarios, (iii) applicability to different spatial and temporal scales, i.e., from lab or field to regional spatial scales and from hourly to decadal temporal scales, and (iv) applicability to different meteorological and climactic conditions, such as humid, arid, and semiarid regions. Three MODFLOW packages accounting for processes in the vadose zone have been recently developed: the Variably Saturated Flow (VSF) process (Thoms et al., 2006), the Unsaturated Zone Flow (UZF1) package (Niswonger et al., 2006), and the HYDRUS package (Seo et al., 2007). Table 1
lists the strengths and weaknesses of selected currently available approaches that incorporate vadose zone flow into MODFLOW. It may be noted that all currently available coupled modeling techniques have some weaknesses and some strengths.
View this table:
[in this window]
[in a new window]
|
TABLE 1. A comparison of the Recharge–Evapotranspiration (REC-ET), Unsaturated Zone Flow (UZF1), and Variably Saturated Flow (VSF) packages for MODFLOW that incorporate the effects of water flow in the vadose zone.
|
|
Among the available packages for MODFLOW, the VSF process most robustly represents the vadose zone processes, as it can consider all major variably saturated flow processes as well as their three dimensionality. The thorough consideration of the vadose zone flow processes in the VSF process, however, makes it computationally very demanding. On the other hand, the UZF1 and REC-ET packages radically simplify vadose zone processes. As a result, these packages are computationally efficient but may not necessarily provide an accurate characterization of vadose zone flow processes.
The objective of this study was to briefly review the aforementioned approaches (i.e., the VSF process and the UZF1 and REC-ET packages) used to account for vadose zone flow in MODFLOW and then to compare them, using case studies of increasing complexity, with the HYDRUS package. The HYDRUS package for MODFLOW was developed to provide a balance between computational efficiency and accuracy. Being one dimensional, the HYDRUS package significantly simplifies the calculations but cannot consider water flow in the vadose zone in multiple dimensions.
It is important to note that when calibrated against collected field data, one may expect the REC-ET and UZF1 packages to perform relatively well for many practical applications. This study looked only at how accurately these approaches (i.e., VSF, HYDRUS, UZF1, and REC-ET) perform relative to each other when field-estimated soil hydraulic parameters or their literature values are given. Calibration, which is not considered here, may lead to estimates of soil hydraulic parameters that do not always correspond to field-estimated values.
|
MODFLOW Packages Accounting for Vadose Zone Processes
|
|---|
The Recharge–Evapotranspiration Package
The Recharge (REC) and Evapotranspiration (ET) packages can be used together to provide a simplistic characterization of vadose zone processes in MODFLOW. The REC package is used to simulate a specified downward recharge flux across the top of the model domain. The recharge flux can be varied spatially and with time. To estimate the volumetric flow rates at the water table, these fluxes are simply multiplied by the horizontal area of the cells. The REC package is a simplified representation of vadose zone flow and does not consider vadose zone flow processes such as storage and runoff. The ET package in MODFLOW is used to simulate the discharge of water to evaporation and transpiration. In the ET package, a maximum evapotranspiration rate is supplied to the model as a function of space and time. To consider the influence of the groundwater depth on evapotranspiration rates, the ET package uses a user-defined extinction depth. While the ET package may be simple to use and understand, it may tend to oversimplify the impact of vadose zone processes. Also, the necessity to supply several rather empirical parameters (e.g., the extinction depth and maximum evapotranspiration rates) increases the uncertainty of modeling results.
The Variably Saturated Flow Process
Incorporating the numerical solution of the three-dimensional Richards equation into the groundwater flow model is the most accurate way to represent the complex nature of physical processes in the unsaturated part of the subsurface. An example of such an approach is the VSF process (Thoms et al., 2006) for MODFLOW. The VSF process solves the three-dimensional form of the Richards equation for the entire MODFLOW domain. In the VSF process, the finite-difference MODFLOW domain is expanded to include the variably saturated zone, and the "mixed form" of the Richards equation is used as the governing equation. The VSF process thus offers more rigorous but much more computationally demanding treatment of water flow in both the unsaturated and saturated zones. The large computational demand stems from the fact that the numerical solution of the Richards equation requires much finer discretization of three spatial dimensions and smaller time steps than traditional groundwater models. This seriously limits the applicability of the VSF process for regional-scale groundwater flow problems (domains >100 km2) (Thoms et al., 2006).
The Unsaturated Zone Flow Package
A number of researchers (e.g., Pikul et al., 1974; Refsgaard and Storm, 1995; Niswonger et al., 2006) have proposed a simpler methodology that significantly decreases the computational demand without greatly compromising the efficiency of the coupled modeling approach. The proposed approach involves coupling a one-dimensional vadose zone flow model with a three-dimensional groundwater flow model (such as MODFLOW). Pikul et al. (1974) and Niswonger et al. (2006) noted that this approach probably provides the most efficient solution for groundwater flow models, especially for large-scale applications. This approach, i.e., consideration of only one-dimensional vertical flow in the unsaturated zone and fully three-dimensional groundwater flow, has been used, for example, in the MIKE SHE model (Refsgaard and Storm, 1995) and the UZF1 package (Niswonger et al., 2006) for MODFLOW.
The UZF1 package couples a vadose zone flow model based on the numerical solution of the one-dimensional kinematic wave equation with MODFLOW. Unlike the Richards equation, which considers both gravity and capillarity as driving forces for flow in the vadose zone, the kinematic wave equation considers only gravity-driven flow. The model, based on the kinematic wave equation, relates water fluxes directly to the degree of saturation. The Brooks and Corey model (Brooks and Corey, 1964), one of the commonly used models relating the moisture content, 
, to the unsaturated hydraulic conductivity, K(), (and the flux, q), is used in the UZF1 package:
 | [1] |
where q is the water flux [L T–1], K() is the unsaturated hydraulic conductivity [L T–1] expressed as a function of the water content (dimensionless), Ks is the saturated hydraulic conductivity (L T–1), 
is the Brooks–Corey exponent (dimensionless), and s and r are the saturated (porosity) and residual water contents (dimensionless), respectively. The UZF1 package considers evaporation and root water uptake (transpiration) by assuming that the water loss occurs instantaneously in the soil profile between the soil surface and a user-specified depth called the ET extinction depth (see Fig. 1).
Application of the kinematic wave equation for vadose zone modeling has its own advantages and disadvantages. Many researchers have debated whether or not variably saturated flow can be treated using the kinematic wave approach (e.g., Singh, 2002). The applicability of the kinematic wave equation (such as in the UZF1 package) to simulate vadose zone flow depends in large part on the soil hydraulic properties, climatic conditions, and the depth to the groundwater table. The UZF1 package requires that the unsaturated zone is homogeneous, which can significantly limit applicability of the package. While in coarse-textured soils, deep vadose zone profiles, or humid climates gravity usually dominates flow in the unsaturated zone and thus the kinematic wave approach is applicable, in fine-textured soils, profiles with shallow groundwater levels, or arid climates, neglecting capillary forces may lead to significant errors. Under such conditions, the kinematic wave equation may fail to describe the dominant flow processes. The UZF1 package can thus be applied mainly in situations where gravity-dominated water flow occurs.
The HYDRUS Package
The HYDRUS package (Seo et al., 2007) was developed for the MODFLOW-2000 (Harbaugh et al., 2000) environment to combine extensive modeling capabilities of both HYDRUS and MODFLOW. The HYDRUS package incorporates into the MODFLOW suite a vadose zone flow model based on the one-dimensional Richards equation. The package was developed to consider the effects of precipitation, infiltration, evaporation, plant water uptake, soil moisture storage, and water accumulation at the ground surface and in the vadose zone. It is based on the HYDRUS-1D program (imnek et al., 2005, 2008), which simulates one-dimensional water movement in the variably saturated zone.
In the coupled HYDRUS–MODFLOW system, vadose zone and groundwater flows are modeled using two separate governing equations. Similarly to the UZF1 package, groundwater flow in MODFLOW is modeled by solving the following mass-conservation equation using a finite-difference approximation:
 | [2] |
where Kx, Ky, and Kz are hydraulic conductivities [L T–1] in the direction of the x, y, and z coordinates, respectively; h is the pressure head [L], W is the volumetric flux per unit volume through sources or sinks [T–1], Ss is the specific storage of the porous material [L–1], and t is time [T].
In the HYDRUS package, vadose zone water flow is described mathematically using the modified Richards equation:
 | [3] |
where is the volumetric water content (dimensionless), h is the soil water pressure head [L], t is time [T], z is the vertical coordinate [L], S is the sink term usually accounting for root water uptake [T–1], and K(h) is the unsaturated hydraulic conductivity [L T–1] as a function of h or . Note that the Richards equation is highly nonlinear due to the dependence of the unsaturated hydraulic conductivity, K(h), and the water content, (h), on the capillary pressure head, h. The two most widely used approaches representing these nonlinear relationships are the Brooks and Corey (Brooks and Corey, 1964) and van Genuchten–Mualem (van Genuchten, 1980) models. Both models are available in the HYDRUS package.
The computer program HYDRUS-1D (imnek et al., 2005) was adapted and simplified for the HYDRUS package. The simplification involved removal of subroutines simulating solute and heat transport, hysteresis in the soil hydraulic functions, and boundary conditions that were irrelevant for the coupled model. The final HYDRUS package thus simulates only one-dimensional water movement in variably saturated porous media. The Galerkin-type linear finite-element scheme is used in HYDRUS to numerically solve the Richards equation.
Water uptake by plant roots has a great effect on water in the root zone. Root water uptake is represented in HYDRUS as an extraction or sink term, S(h), that distributes the potential transpiration across the root zone. Feddes et al. (1978) described the sink term as
 | [4] |
suggesting that the root extraction, S, depends on the pressure head, h, the potential root water uptake rate, Sp [T–1], and the stress response function 
(h), which characterizes plant response to water stresses.
The HYDRUS package for MODFLOW, similarly to the UZF1 and REC-ET packages, does not take into account subsurface runoff because of the one-dimensional nature of the package. The impact of subsurface runoff needs to be considered independently when these packages for MODFLOW are used. On the other hand, the HYDRUS package can consider surface runoff. The HYDRUS package for MODFLOW has an option wherein any excess water on the soil surface can either accumulate there or be immediately removed by surface runoff (imnek et al., 2005).
Spatial Discretization
The efficiency of a coupled vadose zone–groundwater model depends to a large extent on how these two submodels interact with each other in space and time. The MODFLOW model uses the finite-difference approximation of the mass conservation equation to simulate groundwater flow. The groundwater modeling domain for MODFLOW is discretized into grids or blocks as described in Harbaugh et al. (2000) and the number of vadose zone profiles may be as large as the number of rows and columns of this finite-difference grid. Based on similarities in soil hydrology, topographical characteristics, and depth to groundwater, the discretized MODFLOW domain can be divided into zones, which comprise one or more cells of the MODFLOW model (Fig. 2
). One HYDRUS soil profile is then assigned to each of these zones (Fig. 2). It is assumed that the HYDRUS soil profile adequately represents vadose zone flow for the entire zone.
View larger version (47K):
[in this window]
[in a new window]
|
FIG. 2. A discretized aquifer system in MODFLOW and two associated HYDRUS soil profiles. One HYDRUS soil profile is assigned to each MODFLOW zone. Note that the discretization of HYDRUS soil profiles is much finer than that of the MODFLOW domain.
|
|
Created vertical soil profiles are then discretized vertically into finite elements. The finite-element mesh is constructed by splitting the soil profile into one-dimensional elements that are connected to each other at nodal points. Once the finite elements are constructed, they may not be changed during the simulation. Care should therefore be taken to ensure that the depth of the soil profile extends from the soil surface to the deepest possible water table that may be expected during the simulation. To ensure convergence of the numerical solution, finite-element dimensions should be relatively small at locations where sharp pressure head gradients are expected. Such smaller elements are usually needed close to the soil surface where meteorological factors can cause rapid changes in water content and pressure head gradients, and at interfaces between different soil horizons. Soil texture also needs to be considered during the discretization process. For example, coarse-textured soils generally require finer discretization than fine-textured soils due to the higher nonlinearity of their soil hydraulic properties. Once the spatial discretization of the soil profile is performed, the distribution of different soil materials in the profile needs to be described (Fig. 2) (Seo et al., 2007).
Time Discretization
The computational efficiency of the coupled HYDRUS–MODFLOW system is enhanced by simulating vadose zone and groundwater flows at their own, often different, time steps. This is needed because a proper treatment of the Richards equation requires smaller time steps than those usually used in MODFLOW simulations. Figure 3
describes the coupling procedure used in the HYDRUS package. The two models (HYDRUS and MODFLOW) interact, i.e., exchange information about the groundwater recharge and the groundwater level, only at the end of each MODFLOW time step, during which HYDRUS may perform multiple time steps to simulate unsaturated zone flow. MODFLOW receives the recharge flux from HYDRUS and calculates a new water table depth for the next time step. A new water table depth is calculated and assigned as the pressure head bottom boundary condition in the HYDRUS package for the next MODFLOW time step. The iteration procedure in the HYDRUS package is similar to that described in the HYDRUS-1D manual (imnek et al., 2005). See this reference and Seo et al. (2007) for more details.
View larger version (48K):
[in this window]
[in a new window]
|
FIG. 3. Flowchart describing the coupled modeling approach used in the HYDRUS package for MODFLOW: (a) steps shown in gray correspond to the treatment of variably saturated water flow (i, stress period; j, time step; k, soil profile number; nt, number of time steps; ns; number of stress periods; np, number of HYDRUS profiles); (b) calculations carried out by the HYDRUS package during one MODFLOW time step.
|
|
|
Case Studies
|
|---|
The performance of the HYDRUS package was analyzed and compared with other vadose zone flow packages in three case studies: (i) the Las Cruces one-dimensional infiltration experiment (Wierenga et al., 1991), (ii) the two-dimensional water table recharge experiment (Vauclin et al., 1979), and (iii) a hypothetical regional-scale groundwater flow problem. While the HYDRUS, UZF1, and VSF packages were used in the first case study, only the HYDRUS and UZF1 packages were applied in the second case study since results for the VSF process for this application can be found in the literature. Finally, the REC-ET, UZF1, and HYDRUS packages were used in the third case study.
The one-dimensional Las Cruces infiltration experiment of Wierenga et al. (1991) was used first to evaluate the effectiveness of the HYDRUS, VSF, and UZF1 packages to simulate flow in the vadose zone without considering groundwater flow. The water table recharge experiment of Vauclin et al. (1979) was then used to evaluate whether a combination of a one-dimensional vadose zone module with a groundwater model can approximate this obviously two-dimensional problem. Finally, a complex regional-scale groundwater flow problem was used to evaluate the effectiveness of different vadose zone packages in accounting for various processes in the vadose zone.
Case Study 1: One-Dimensional Infiltration Experiment
The first case study involved the one-dimensional infiltration experiment at the Las Cruces trench site (Wierenga et al., 1991). The experiment involved a comprehensive field study, conducted in southern New Mexico, the primary purpose of which was to develop a data set for validating and testing numerical models. For this purpose, the study site was heavily instrumented with neutron probes, tensiometers, and solute samplers for measuring water contents, pressure heads, and solute concentrations (Wierenga et al., 1991), respectively. More than 500 soil samples (undisturbed and disturbed) were taken at the experimental site and analyzed in the laboratory for bulk density and to find the saturated hydraulic conductivity and the soil water retention curve. The infiltration study involved application of water to a 4-m-wide area using closely spaced drips with an average surface flux of 1.82 cm d–1 for 86 d of the experiment. To reduce the disruption of the experimental conditions by rain and evaporation, the irrigated area and its surroundings were covered by a pond liner.
Wierenga et al. (1991) performed a one-dimensional simulation of the infiltration experiment using a numerical model based on the finite-difference approximation of the Richards equation. They considered a uniform soil profile with an equivalent saturated hydraulic conductivity Ks of 270.1 cm d–1. The RETC code (van Genuchten et al., 1991) was used to analyze the retention curve data for undisturbed and disturbed soils (>500 soil samples), resulting in the following retention curve parameter values (van Genuchten, 1980): s = 0.321, r = 0.083, = 0.055 cm–1, and n = 1.51. These values were then used by the HYDRUS package and VSF process. A zero extinction depth was used in the UZF1 package, as no evaporation losses were considered. Morel-Seytoux et al. (1996) developed equations describing the parameter equivalence between the Brooks–Corey exponent and van Genuchten parameters. From these equations, the Brooks–Corey exponent for the UZF1 package was estimated to be 6.92.
The one-dimensional simulation of the infiltration experiment was performed using MODFLOW with the HYDRUS, VSF, and UZF1 packages. The finite difference mesh for MODFLOW consisted of a one-cell grid. Initial pressure heads (hi = –100 cm) in the soil profile were the same as those used by Wierenga et al. (1991). While a constant water flux was used as the upper boundary condition (q0 = 1.82 cm d–1), free drainage was considered at the lower boundary. The implicit assumption in this boundary condition is that the groundwater table is deep enough so that it does not affect flow in the soil profile. The initial and boundary conditions in terms of the water content, (z, t), are described as follows:
 | [5a] |
 | [5b] |
 | [5c] |
where b is the depth of the soil profile, which must be large enough so that the wetting front does not affect the water content at the bottom of the soil profile during the simulation, init(z) and hinit(z) are initial water contents and pressure heads at depth z, respectively. A soil profile depth of 600 cm was used and the simulation was run for 35 d.
Experimental results of the Las Cruces trench infiltration experiment (Wierenga et al., 1991) are compared with results simulated using the HYDRUS, VSF, and UZF1 packages in Fig. 4. Figure 4
shows the soil water content profiles for different days of the experiment and compares model predictions of the VSF, UZF1, and HYDRUS packages with the experimental data. As expected, the HYDRUS package and the VSF process performed similarly as they both solve the same Richards equation for one-dimensional problems. The UZF1 package only slightly overpredicted water contents behind the wetting front. A comparison of the time needed for the simulations by the VSF process and the UZF1 and HYDRUS packages was done. It was noted that the computational demand of the VSF, UZF1, and HYDRUS packages was similar for the one-dimensional case study.
View larger version (19K):
[in this window]
[in a new window]
|
FIG. 4. A comparison of measured water contents and corresponding estimates calculated using the HYDRUS, Variably Saturated Flow (VSF), and Unsaturated Zone Flow (UZF1) packages for (a) Day 19 and (b) Day 35 for the Las Cruces trench experiment (data from Wierenga et al., 1991).
|
|
The UZF1, HYDRUS, and VSF packages provided similar results for this one-dimensional infiltration experiment and needed comparable computational times.
Case Study 2: Two-Dimensional Water Table Recharge Experiment
The HYDRUS and UZF1 packages were used to model the two-dimensional transient water table experiment of Vauclin et al. (1979). The same data set was used previously by Thoms et al. (2006) to evaluate the VSF process. See these references for their results.
The experimental setup consisted of a 6.0- by 2.0-m box containing sandy soil. The initial water table elevation was 0.65 m from the bottom. A constant flux of q = 3.55 m d–1 was applied across the center 1.0 m of the soil surface while the rest of the surface was covered to prevent evaporation. Due to the symmetry of the experiment, only one half of the experiment was modeled and the model domain was thus 3.0 by 2.0 m. The initial total head was set equal to 0.65 and the right boundary cells were constrained to the initial water table position throughout the 8-h simulation. The grid was discretized into uniform cells of 0.1-m width and 0.05-m depth. For comparison purposes, simulations were performed using the exact setup described in the documentation of the VSF process. Only two soil profiles representing the soil directly below the recharge zone and the rest of the transport domain were used in calculations with the HYDRUS and UZF1 packages. The saturated hydraulic conductivity of 840 cm d–1 was used. The initial and boundary conditions are described as follows:
Domain
Initial condition
 | [6] |
Boundary condition
where q(x,z,t) is the flux at spatial coordinates x and z at time t. The Brooks–Corey exponent was set equal to 6.37 in the UZF1 model, based on the estimate from Carsel and Parrish (1988). A zero ET extinction depth was assumed, as the experiment was designed to minimize all evaporative losses. The HYDRUS package used retention curve parameters similar to those used for the VSF process (Thoms et al., 2006; Vauclin et al., 1979), i.e., s = 0.30, r = 0.01, = 0.033 cm–1, and n = 4.1.
Figure 5
compares water tables simulated using the HYDRUS and UZF1 packages with the experimental data. One may also compare the performance of the VSF process by referring to Thoms et al. (2006). Water tables calculated using the HYDRUS package are similar to those simulated using the VSF process even though the numerical solution of the Richards equation in the HYDRUS package is limited to only the vertical direction. It was observed that the one-dimensional nature of the vadose zone modeling used in the HYDRUS package did not significantly affect the correspondence of simulated results with experimental data. Note that while only vertical flow was allowed in the vadose zone, horizontal flow below the water table redistributed recharged water and resulted in smooth water tables; however, a comparison of results calculated with the UZF1 package (Fig. 5b) with those obtained using the VSF and HYDRUS packages shows that the UZF1 package marginally underestimated the depth of the water table. This may be attributed to the kinematic wave approximation used in the UZF1. Also, the uncertainty in soil hydraulic parameters may be responsible for the differences. The underprediction of the water table depth by the UZF1 package is relatively larger at later times (i.e., 8 h). The calibrated UZF1 package would probably provide similar results to those by HYDRUS and VSF.
A comparison of computational times needed for the VSF, HYDRUS, and UZF1 showed that the UZF1 package required the least computational effort, followed by HYDRUS and VSF. Simulation times required by the HYDRUS and UZF packages were, however, significantly smaller than required by the VSF process.
It can be concluded that for small-scale groundwater problems (up to two dimensions, such as in Case Studies 1 and 2) with downward flow in the vadose zone, the UZF1 package has accuracy similar to the HYDRUS package and the VSF process, at least for certain cases such as those considered here. The third case study was designed to test the performance of the HYDRUS package for a regional-scale groundwater flow problem.
Case Study 3: Hypothetical Regional-Scale Groundwater Problem
The third case study involved a hypothetical large-scale groundwater flow problem in a semiarid to arid region. The geometry of the modeling domain (Fig. 6
) was based on the test example described in Prudic et al. (2004) and Niswonger et al. (2006). In this case study, we compared the effectiveness of the REC-ET, UZF1, and HYDRUS packages in characterizing vadose zone processes at a regional scale. The VSF process was not used here because of its extraordinary computational demand (Thoms et al., 2006) for such large-scale applications.
View larger version (38K):
[in this window]
[in a new window]
|
FIG. 6. Model domain, spatial distribution of hydraulic conductivities and specific yields, wells (red circles), and general head boundaries for the hypothetical regional-scale groundwater flow problem.
|
|
The model domain was designed to represent an alluvial basin with loamy soils. Figure 6 shows the model domain and other key characteristics for this hypothetical regional-scale groundwater flow problem. The flow domain was divided into uniform grids of 1524- by 1524-m size. Two cells were assigned a general head boundary condition to simulate head-dependent flux boundaries to allow flow in and out of the system. At the head-dependent flux boundary, water enters the model domain if the head in the cell is less than a certain user-defined reference head and leaves the model domain otherwise. The alluvium valley aquifer was assumed to have greater hydraulic conductivity than the upland areas. Figure 7
shows surface elevations, bedrock depth, and initial water table depths in the study area. The geologic settings were varied spatially to present a complex three-dimensional case.
View larger version (56K):
[in this window]
[in a new window]
|
FIG. 7. (a) Land surface elevation, (b) depth to bedrock, and (c) water table depth at the beginning of the simulation for the hypothetical regional-scale groundwater flow problem.
|
|
The modeling time was divided into 12 equal stress periods, each of which lasted 30.42 d. Except for the first stress period, they were modeled in MODFLOW in the transient mode with 15 time steps for each stress period. The first stress period was modeled as steady state. The meteorological conditions were assigned to represent a semiarid climate where potential evaporation rates are substantially higher than precipitation rates. Such meteorological settings require consideration of both downward and upward water fluxes in the soil profile and provide, therefore, a good case study to compare the HYDRUS package to the UZF1 package. While Table 2
provides the base precipitation, potential evaporation, and well pumping rates for the 12 stress periods, Fig. 8
shows precipitation rate factors that were used to vary precipitation rates throughout the flow domain. While different precipitation rates were assigned for each stress period, the spatial distribution of precipitation rates was considered to be the same for all stress periods.
View this table:
[in this window]
[in a new window]
|
TABLE 2. Precipitation, potential evaporation, and pumping rates for a hypothetical regional-scale groundwater flow problem.
|
|
View larger version (20K):
[in this window]
[in a new window]
|
FIG. 8. Zonation showing the spatial distribution of precipitation for the study area of the hypothetical groundwater flow problem. For any stress period, the actual precipitation rate in the zone is obtained by multiplying the precipitation rates given in Table 2 by the zone precipitation rate factors.
|
|
We assumed that the vadose zone consisted of loamy soils throughout the model domain. The following Brooks–Corey (Brooks and Corey, 1964) model parameters were used in the UZF1 package, as suggested for loam by Carsel and Parrish (1988): s = 0.30, r = 0.00, Ks = 3 x 10–5 cm s–1, inverse of the air entry pressure = 0.0897 cm–1, the pore size distribution index n = 0.22, and the pore-connectivity parameter l = 1. An ET extinction depth of 2.65 m (Shah et al., 2007) was used and the Brooks–Corey exponent, , was calculated to be 12.09.
The HYDRUS package offers a variety of models for characterizing the water retention and unsaturated hydraulic conductivity of the soil. The van Genuchten–Mualem model (van Genuchten, 1980) was used in the HYDRUS package to represent the unsaturated hydraulic conductivity and water content dependency on the capillary pressure. As suggested by Carsel and Parrish (1988), the following van Genuchten parameters for loam were used: s = 0.30, r = 0.00, Ks = 3 x 10–5 cm s–1, = 0.036 cm–1, n = 1.56, and l = 0.5.
For both HYDRUS and UZF1 packages, it was essential to first create zones that represented relatively homogeneous units with similar soil and hydrogeologic properties so that one soil profile could be assigned to each zone. To create the zones, the fuzzy c-means clustering algorithm of Dunn (1973) was used. The fuzzy c-means algorithm is a method of clustering that allows one piece of data to be split into a user-defined number of clusters such that the data points in each cluster are as similar as possible. The fuzzy c-means method was used independently of the MODFLOW–HYDRUS environment as it is not a part of it. The zones were created based on surface elevations, hydraulic conductivities, initial water table heads, and locations of cells. The number of clusters was chosen to be 20 since additional clusters did not significantly improve uniformity within each cluster (MODFLOW zone). Figure 9
shows the zones used in the HYDRUS package for the hypothetical vadose zone–aquifer interaction problem. The same zones were used for the UZF1 package.
View larger version (39K):
[in this window]
[in a new window]
|
FIG. 9. MODFLOW zones used to define HYDRUS soil profiles in the hypothetical groundwater flow problem.
|
|
The HYDRUS soil profiles corresponding to each of these zones were then created and discretized vertically into finite elements. Figure 10
shows initial soil water contents for selected soil profiles. One may note a large variation of initial soil water contents between different soil profiles even though the same soil texture was used throughout the transport domain (i.e., loamy soils). This was due to differences in other variables, such as the water table depth and the elevation of the land surface. In other packages, such as REC-ET and UZF1, it is difficult to take into account the variability of soil profiles.
View larger version (16K):
[in this window]
[in a new window]
|
FIG. 10. Initial water contents as a function of depth in HYDRUS soil profiles representing selected MODFLOW zones (colors correspond to zones in Fig. 9) in the hypothetical groundwater flow problem.
|
|
One of the benefits of the HYDRUS package is its ability to estimate transient fluxes in the groundwater table based on the hydrologic conditions and water residence time in the soil column without significantly compromising the computational requirements. Changes in the magnitude and direction of the flux can be caused by changes in water contents in the HYDRUS soil profiles as a result of time-variable surface meteorological conditions and the position of the water table. Figure 11
shows the groundwater flux zones estimated for different stress periods. Note that, depending on various hydrologic and topological conditions, the HYDRUS package predicts both positive (downward) recharge and negative (upward, capillary rise) discharge fluxes. Water table fluxes at any cell are directly influenced by surface infiltration, evaporation, and transpiration, as well as pumping rates in and around a particular cell. During the initial stress periods, the HYDRUS package predicted considerable upward fluxes, especially in cells where the precipitation rates were lower (compare Fig. 11 with Fig. 8). Cells with deeper initial water tables were less affected by evaporation than those with shallow initial water tables. As the simulation time proceeded, the number of cells with upward fluxes decreased because of the infiltration front movement toward the water table.
View larger version (25K):
[in this window]
[in a new window]
|
FIG. 11. Groundwater table fluxes (recharge vs. discharge) as predicted by the HYDRUS package at the end of Stress Periods (a) 3, (b) 6, and (c) 12 for the hypothetical groundwater flow problem.
|
|
Figure 12
compares the final water table depths estimated by the three packages. Water table depths estimated by the UZF1 and HYDRUS packages were comparable except that the HYDRUS package consistently predicted marginally deeper water tables. Water table depths predicted using HYDRUS at the end of the simulation were deeper by between 0 and 1.31 m than those calculated with the UZF1 package. This is probably due to a more accurate characterization of capillary pressures and fluxes in the vadose zone by HYDRUS. The HYDRUS package can consider upward pressure gradients that cannot be simulated by either the UZF1 or the REC-ET packages, both of which tend to predict greater downward water table fluxes.
View larger version (29K):
[in this window]
[in a new window]
|
FIG. 12. Depths to the water table for the groundwater flow problem at the end of Stress Period 12 as estimated by the (a) Recharge–Evapotranspiration (REC-ET), (b) Unsaturated Zone Flow (UZF1), and (c) HYDRUS packages.
|
|
Figure 13
shows water table depths at the end of different stress periods as a function of initial water table depths. This figure indicates the impact of vadose zone flow on model predictions of groundwater tables. Low capillary pressures in soil profiles of arid zones often lead to slow upward water movement. Final water table depths predicted using the HYDRUS package deviate from those predicted by the other packages, especially for cells with deep initial water tables. One may attribute these differences to the influence of capillary forces on vadose zone flow, which are fully considered only in the HYDRUS package. Groundwater table depths calculated using the UZF1 package also deviate from those calculated using the REC-ET package, albeit in some cells to a lesser extent than those predicted using the HYDRUS package. This is due to the delay in groundwater table fluxes resulting from the gravity-driven flow considered by the UZF1 package and neglected by the REC-ET package. One can infer that vadose zone processes resulting from different soil hydraulic characteristics can be better represented using the HYDRUS package than the REC-ET and UZF1 packages, especially for soils, such as medium- and fine-textured soils, where capillary forces play an important role in determining water fluxes.
View larger version (24K):
[in this window]
[in a new window]
|
FIG. 13. Depth to the water table estimated using the Recharge–Evapotranspiration (REC-ET), Unsaturated Zone Flow (UZF1), and HYDRUS packages at the end of Stress Periods (a) 3, (b) 6, (c) 9, and (d) 12 as a function of the initial water table.
|
|
The computational time needed to perform calculations with these vadose zone packages is of paramount importance because it seriously affects their applicability to regional-scale problems. The UZF1, REC-ET, and HYDRUS packages were run on a 1 GB RAM, 3.40 GHz Intel Pentium based personal computer. The simulation of the hypothetical regional-scale groundwater model (Case Study 3) using the REC-ET package needed, as expected, the least amount of time (5 s). The MODFLOW simulations with the UZF1 and HYDRUS packages took approximately 20 and 26 s, respectively. While the computational demand of the UZF1 and HYDRUS packages is comparable, the HYDRUS package provides a more comprehensive characterization of vadose zone flow processes.
|
Summary and Conclusions
|
|---|
We evaluated the recently developed HYDRUS package for MODFLOW and compared its performance with other MODFLOW packages (REC-ET, UZF1, and VSF) that account for processes in the vadose zone. Based on the HYDRUS-1D software, the HYDRUS package considers the effects of infiltration, soil moisture storage, evaporation, plant water uptake, precipitation, runoff, and water accumulation at the ground surface. The HYDRUS package was compared with other currently available packages (VSF and UZF1) for case studies of varying complexities. For smaller scale problems (up to two dimensions), the VSF process and the UZF1 and HYDRUS packages perform similarly. The VSF process provides the best accuracy due to its thorough consideration of vadose zone processes and is most suitable for smaller scale problems. For three-dimensional models, the HYDRUS package demonstrated a significant improvement in modeling accuracy compared with the UZF1 package and a significant decrease in computational demand compared with the VSF process. This new package represents a promising tool accounting for the effects of vadose zone processes on groundwater levels and fluxes.
Even though the UZF1 package provides comparable results to the HYDRUS package, a number of difficulties arise during its application. Although previous research (e.g., Shah et al., 2007) provides ET extinction depths for different textural classes, determination of the ET extinction depth for practical cases can be cumbersome. The presence of vegetation further complicates the determination of the ET extinction depth. In contrast, the HYDRUS package and the VSF process use soil hydraulic and plant parameters, such as the van Genuchten and Feddes parameters, respectively, that are readily available for a wide variety of soil textures and plants (Leij et al., 1996; Lilly, 1997). The HYDRUS package also offers a variety of models for describing soil hydraulic properties. The necessary parameter values required for these models in HYDRUS may be estimated experimentally or using pedotransfer functions (e.g., Schaap et al., 2001).
Another weakness of the current UZF1 package is that it simulates unsaturated flow only for homogeneous vadose zones and it cannot consider multiple soil horizons with varying hydraulic properties. User-specified layering can be easily accommodated when using the HYDRUS package. The HYDRUS package thus offers a good alternative to the UZF1 package when these factors, e.g.., vegetation or multiple horizons, are significant for a particular application. The HYDRUS package thus may be expected to perform better for regional-scale groundwater problems with complex layering in the vadose zone and with alternating recharge–discharge fluxes.
The HYDRUS package currently simulates only water flow and is distributed as an open-source code. We intend to expand the HYDRUS package to also simulate solute transport so that the MODFLOW–HYDRUS code will produce concentrations as a function of time that can be incorporated into the source function for MT3D. This is expected to be especially useful for regional-scale studies involving nonpoint-source pollution.
|
REFERENCES
|
|---|
- Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrol Pap. 3. Colorado State Univ., Fort Collins.
- Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24:755–769.[CrossRef]
- Colbeck, S.C. 1972. A theory of water percolation in snow. J. Glaciol. 2:369–385.
- Dunn, J.C. 1973. A fuzzy relative of the ISODATA process and its use in detecting compact well-separated clusters. J. Cybernetics 3:32–57.[CrossRef]
- Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of field water use and crop yield. Simulation Monogr. PUDOC, Wageningen, the Netherlands.
- Gee, G.W., and D. Hillel. 1988. Groundwater recharge in arid regions: Review and critique of estimation methods. Hydrol. Processes 2:255–266.[CrossRef]
- Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald. 2000. MODFLOW-2000, the U.S. Geological Survey modular ground-water model user guide to modularization concepts and the ground-water flow process. USGS, Denver, CO.
- Hendrickx, J., and G. Walker. 1997. Recharge from precipitation. p. 19–38. In I. Simmers (ed.) Recharge of phreatic aquifers in (semi-) arid regions. A.A. Balkema, Rotterdam, the Netherlands.
- Keese, K.E., B.R. Scanlon, and R.C. Reedy. 2005. Assessing controls on diffuse groundwater recharge using unsaturated flow modeling. Water Resour. Res. 41:W06010, doi:10.1029/2004WR003841.[CrossRef]
- Leij, F.J., W.J. Alves, M.Th. van Genuchten, and J.R. Williams. 1996. The UNSODA—Unsaturated Soil Hydraulic Database: User's manual version 1.0. Rep. EPA/600/R-96/095. USEPA Natl. Risk Manage. Res. Lab., Cincinnati, OH.
- Lerner, D.N., A.S. Issar, and I. Simmers. 1990. Groundwater recharge: A guide to understanding and estimating natural recharge. Rep. 8. Int. Assoc. Hydrogeol., Kenilworth, UK.
- Lilly, A. 1997. A description of the HYPRES database (Hydraulic Properties of European Soils). p. 161–184. In A. Bruand et al. (ed.) The use of pedotransfer functions in soil hydrology research. Proc. Worksh. of the Project "Using Existing Soil Data to Derive Hydraulic Parameters for Simulation Modelling in Environmental Studies and in Land Use Planning," 2nd, Orleans, France. 10–12 Oct. 1996. INRA, Orleans.
- Lorenz, D.L., and G.N. Delin. 2007. A regression model to estimate regional ground water recharge. Ground Water 45:196–208.[CrossRef][Web of Science][Medline]
- Morel-Seytoux, H.J., P.D. Meyer, M. Nachabe, J. Touma, M.Th. van Genuchten, and R.J. Lenhard. 1996. Parameter equivalence for the Brooks–Corey and van Genuchten soil characteristics: Preserving the effective capillary drive. Water Resour. Res. 32:1251–1258.[CrossRef]
- Niswonger, R.G., D.E. Prudic, and R.S. Regan. 2006. Documentation of the Unsaturated-Zone Flow (UZF1) package for modeling unsaturated flow between the land surface and the water table with MODFLOW-2005. Techniques and Methods 6-A19. Available at http://pubs.usgs.gov/tm/2006/tm6a19/ (verified 14 Apr. 2008). USGS, Denver, CO.
- Pikul, M.F., R.L. Street, and I. Remson. 1974. A numerical model based on coupled one-dimensional Richards and Boussinesq equations. Water Resour. Res. 10:295–302.[CrossRef]
- Prudic, D.E., L.F. Konikow, and E.R. Banta. 2004. A new Stream-Flow Routing (SFR1) package to simulate stream–aquifer interaction with MODFLOW-2000. Open-File Rep. 2004-1042. USGS, Carson City, NV.
- Refsgaard, J.C., and B. Storm. 1995. MIKE SHE. p. 809–846. In V.P. Singh (ed.) Models of watershed hydrology. Water Resour. Publ., Highlands Ranch, CO.
- Richards, L.A. 1931. Capillary conduction of liquids through porous mediums. Physics 1:318–333.[CrossRef]
- Scanlon, B.R. 2002. Choosing appropriate techniques for quantifying groundwater recharge. Hydrogeol. J. 10:18–39.[CrossRef]
- Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251:163–176.[CrossRef]
- Seo, H.S., J. imnek, and E.P. Poeter. 2007. Documentation of the HYDRUS package for MODFLOW-2000, the U.S. Geological Survey modular ground-water model. GWMI 2007-01. Int. Ground Water Modeling Ctr., Colorado School of Mines, Golden.
- Shah, N., M. Nachabe, and M. Ross. 2007. Extinction depth and evapotranspiration from ground water under selected land covers. Ground Water 45:329–338.[CrossRef][Web of Science][Medline]
- imnek, J., M.Th. van Genuchten, and M. ejna. 2005. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media, Version 3.0. HYDRUS Softw. Ser. 1. Dep. of Environ. Sci., Univ. of California, Riverside.
- imnek, J., M.Th. van Genuchten, and M. ejna. 2008. Development and applications of the HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7:587–600 (this issue).[Abstract/Free Full Text]
- Singh, V.P. 2002. Is hydrology kinematic? Hydrol. Processes 16:667–716.[CrossRef]
- Smith, R.E. 1983. Approximate sediment water movement by kinematic characteristics. Soil Sci. Soc. Am. J. 47:3–8.[Abstract/Free Full Text]
- Smith, R.E., and R.H.B. Hebbert. 1983. Mathematical simulation of interdependent surface and subsurface hydrologic processes. Water Resour. Res. 19:987–1001.[CrossRef]
- Thoms, R.B., R.L. Johnson, and R.W. Healy. 2006. User's guide to the Variably Saturated Flow (VSF) process for MODFLOW. Techniques and Methods 6-A18. Available at http://pubs.usgs.gov/tm/2006/tm6a18/ (verified 15 Apr. 2008). USGS, Reston, VA.
- Uddameri, V., and M. Kuchanur. 2007. Estimating aquifer recharge in Mission River watershed, Texas: Model development and calibration using genetic algorithms. Environ. Geol. 51:897–910.[CrossRef]
- 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.[Abstract/Free Full Text]
- van Genuchten, M.Th., S.M. Lesch, and S.R. Yates. 1991. The RETC code for quantifying the hydraulic functions of unsaturated soils. Version 1.0. U.S. Salinity Lab., Riverside, CA.
- Vauclin, M., J. Khanji, and G. Vachaud. 1979. Experimental and numerical study of a transient, two-dimensional unsaturated–saturated water table recharge problem. Water Resour. Res. 15:1089–1101.[CrossRef]
- Ward, S. 2002. Recharge and groundwater models: An overview. Hydrogeol. J. 10:110–120.[CrossRef]
- Wierenga, P.J., R.G. Hills, and D.B. Hudson. 1991. The Las Cruces trench site: Characterization, experimental results, and one-dimensional flow predictions. Water Resour. Res. 27:2695–2705.[CrossRef]
- Winter, T.C., J.W. Harvey, O.L. Franke, and W.M. Alley. 1998. Ground water and surface water a single resource. Circ. 1139. USGS, Denver, CO.
This article has been cited by other articles:
|
|
|
|
 
J. Simunek and S. A. Bradford
Vadose Zone Modeling: Introduction and Importance
Vadose Zone J.,
May 27, 2008;
7(2):
581 - 586.
[Full Text]
[PDF]
|
|
|
|
|
|
|
J. Simunek, M. Th. van Genuchten, and M. Sejna
Development and Applications of the HYDRUS and STANMOD Software Packages and Related Codes
Vadose Zone J.,
May 27, 2008;
7(2):
587 - 600.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
A. Furman
Modeling Coupled Surface-Subsurface Flow Processes: A Review
Vadose Zone J.,
May 27, 2008;
7(2):
741 - 756.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
