HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Schoups, G. Articles by Hopmans, J. W. Search for Related Content PubMed Articles by Schoups, G. Articles by Hopmans, J. W. GeoRef GeoRef Citation Agricola Articles by Schoups, G. Articles by Hopmans, J. W. Related Collections Vadose Zone Processes and Chemical Transport Ground Water Quality Irrigation

Analytical Model for Vadose Zone Solute Transport with Root Water and Solute Uptake

Hydrologic Sciences, Department of Land, Air, and Water Resources, University of California, One Shields Avenue, Davis, CA 95616
* Corresponding author (jwhopmans{at}ucdavis.edu)

Received 3 December 2001.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSION REFERENCES

 
An efficient method is presented for calculating one-dimensional solute transport through the vadose zone in the presence of vertically distributed root water and solute uptake. The method is based on an analytical solution of the convective transport equation, which is solved by the method of characteristics. The solution requires a time-invariant leaching fraction, and assumes that transport occurs by convection alone; that is, hydrodynamic dispersion and molecular diffusion are neglected. From this solution, two variables of practical importance are calculated, (i) the average solute concentration in the root zone, weighted by the root distribution, and (ii) the cumulative solute flux to the groundwater. Model parameters are related to water management, crop type, and soil type, through the following parameters: leaching fraction, root distribution, and soil hydraulic properties. Simulation results are presented for both downward and upward flow scenarios. Simulations with a nonuniform moisture profile indicate that nearly identical results are obtained by replacing the profile by a uniform moisture content equal to the arithmetic average or the uptake-weighted average moisture content, so that closed-form analytical expressions for travel time and travel depth can be obtained. A sensitivity analysis of the relevant model parameters shows that solute concentrations are largely determined by the magnitude of the effective water application ratio, defined by the fraction of infiltrated water that is removed by evapotranspiration. We present suggestions of how dispersive mixing and temporal changes in leaching fraction are readily incorporated into the current model. Although potentially simplistic, the presented methodology can be extremely useful for longer-term and field-to-regional scale characterization of contaminant movement.

Abbreviations: BVP, boundary value problem • IVP, initial value problem

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSION REFERENCES

 
THE MOVEMENT OF POLLUTANTS from diffuse sources applied at the surface is recognized as one of the major contributors to soil and water contamination worldwide (Duda, 1993). One important example is soil and groundwater salinization in irrigated agriculture. Sustainable management of soil and water resources depends largely on the ability to predict the impacts of current and future management practices on soil and groundwater quality. Due to the large time and spatial scales involved, assessment of non-point source pollution is a challenging task. An integrated system of advanced information technologies, including GIS, remote sensing, and solute transport modeling is needed to cope with this complex environmental problem (Corwin et al., 1999). Solute transport models are very useful, because they represent causal relationships between land use decisions and soil and water contamination levels, which may not be immediately apparent in the field, for example, due to large solute travel times. Furthermore, solute transport models provide a means of predicting impacts of future management decisions. Jury and Tseng (1999) stressed the need for a compromise between deterministic–mechanistic and statistical–empirical approaches in choosing a suitable transport model, especially for large-scale hydrologic systems. For example, a one-dimensional model that includes the most relevant processes and parameters at the space and time scale of interest can be applied to the local scale, and distributed over the larger region, either stochastically (e.g., Bresler and Dagan, 1983; Toride and Leij, 1996) or deterministically by its coupling to a GIS (e.g., Corwin and Wagenet, 1996; Loague et al., 1998). This approach may be referred to as streamtube modeling.

Solute transport through the vadose zone is affected by a large number of factors, such as soil physical and chemical properties, plant water and solute uptake rates, and water and solute application rates. In areas with a shallow water table, capillary rise of water and solutes from the water table into the root zone is also an important factor. Water and solute uptake by plants affects root zone concentrations and also the flux and travel time of water and solutes to the water table. For example, evapoconcentration of the soil solution by root water uptake is the main mechanism of salt buildup in irrigated lands (Tanji, 1990). The effect of root water uptake on vadose zone transport commonly is described by the leaching fraction. This dimensionless parameter is defined as the fraction of surface-infiltrated water that drains below the root zone. Root zone salt concentrations have been shown to be strongly dependent on the leaching fraction, with low leaching fractions resulting in high root zone concentrations (Hoffman and van Genuchten, 1983).

Root water uptake has been modeled from both a microscopic and a macroscopic perspective. The microscopic approach considers radial flow to a single root of specified geometry (Gardner, 1960). Alternatively, the macroscopic approach (e.g., Nimah and Hanks, 1973; Molz, 1981; Clausnitzer and Hopmans, 1994) avoids these geometric complications and incorporates root water uptake into continuum water flow models through a distributed sink term. This sink term, r_w, can be represented in one dimension as (Simunek et al., 1998),

[1]

where z denotes spatial location within the root zone, t is time, is the stress response function as a function of matric head (h) and osmotic head (), T_pot represents the potential transpiration rate by the plant, and r(z) is a root distribution function. Several functional forms have been proposed for r: uniform (Feddes et al., 1976), linear (Molz and Remson, 1970; Prasad, 1988), trapezoidal (Hoffman and van Genuchten, 1983), exponential (Raats, 1974b) as well as variations thereof (Dwyer et al., 1988; Li et al., 1999; Vrugt et al., 2001), and a power law function (Ojha and Rai, 1996). Reductions in root water uptake due to water and salinity stress are represented by empirical functions, (h,), which may be linear (Feddes et al., 1978) or nonlinear (van Genuchten, 1987). Homaee (1999) and Homaee and Feddes (1999) provided a detailed discussion of reductions in water uptake as a result of combined water and salinity stress.

Root solute uptake has also been modeled from a microscopic (Nye and Marriott, 1969) and a macroscopic viewpoint (Somma et al., 1998). The macroscopic approach incorporates a sink term in the transport equation. For root solute uptake, a distinction is made between passive and active uptake mechanisms (Epstein, 1960). Passive solute uptake occurs by both advection and diffusion from high to low concentrations. Active uptake, on the other hand, is driven by specific energy-driven ion carriers or occurs through ion channels embedded in the root cell membrane (Marschner, 1995). Active uptake requires energy from plant respiration and allows solute uptake against a concentration gradient. In their study of coupled uptake of water and solutes to a single root for steady flow conditions, Dalton et al. (1975) showed that the rate of solute uptake is proportional to the water uptake rate, even when active uptake is dominant.

Several numerical models of vadose zone solute transport are available that account for distributed root water and solute uptake using the macroscopic approach. These models solve the Richards equation with a sink term representing root water uptake, and the convection–dispersion equation for transport, in either one dimension (Simunek et al., 1998), two dimensions (Simunek et al., 1999), or three dimensions (Somma et al., 1998). Although these models are very attractive at the laboratory or plot scale, their use is more challenging at the field to regional scale, with time scales spanning years or decades. Therefore, it is appropriate to develop and use simplified models that focus on the main processes operating at the pertinent temporal and spatial scales. The error introduced by using a simpler model is assumed to be small compared with the large uncertainties associated with estimating the time-varying boundary conditions (e.g., evapotranspiration) and spatially variable soil parameter values (Chen et al., 1994).

This paper presents a simple, yet efficient and robust, method for calculating one-dimensional solute transport through the vadose zone in the presence of vertically distributed root water and solute uptake. The method is based on an analytical solution of the convective transport equation, which is derived by the method of characteristics. The main assumptions are that the leaching fraction is time-invariant and that transport occurs by convection alone (i.e., hydrodynamic dispersion and molecular diffusion are neglected). These theoretical results generalize and extend the collective results of Raats (1975) and Ginn and Murphy (1997), who both studied convective transport of a tracer through the root zone. The proposed model differs from other analytical solutions (van Genuchten and Alves, 1982; Leij et al., 1991; Leij and Toride, 1998) in that these do not consider distributed root water or solute uptake. The model's intended application is to represent local transport in a streamtube model for long-term, regional-scale predictions of solute transport, specifically as related to soil salinity in both shallow and deep vadose zones. As such, it may provide an attractive alternative to more involved methods based on a numerical solution of the Richards and convection–dispersion equations (e.g., Simunek et al., 1998). The model may also be useful for estimating fluxes of water and solute from solute concentration measurements. For example, Rhoades et al. (1999) indicated the need for a model that allows for temporal variation of applied water and salts to estimate leaching fraction and salt loadings from soil salinity measurements.

The paper is organized as follows. In the first part, we introduce the governing equations and assumptions, and an analytical solution is derived using the method of characteristics. In the second part, we present simulation results for downward and upward flows and explore the sensitivity of the model results to different parameters. Finally, we discuss how the model can be extended to account for time-varying leaching fraction and dispersive mixing.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSION REFERENCES

 
Governing Equations and Assumptions
The continuity equation for one-dimensional vertical flow through the root zone is

[2]

where (z,t) is volumetric water content, q(z,t) is soil water flux (L/T, positive downwards), r_w(z,t) is a sink term representing root water uptake and evaporation (1/T), t is time (T), and z is vertical coordinate (L, positive downwards). The continuity equation for one-dimensional transport of a solute undergoing root uptake and equilibrium sorption, while neglecting dispersion and diffusion, is given by

[3]

where c is the dissolved solute concentration (M/L³), expressed as mass of solute per volume of soil solution; c_a is adsorbed solute concentration (M/M), expressed as mass of sorbant per mass of dry soil; _b is dry bulk density (M/L³); and r_s is a sink term representing root solute uptake [M/(L³T)]. Combining Eq. [2] and [3] leads to (Appendix A),

[4]

where v_R = v/R = q/ is the retarded solute velocity (L/T), v(z,t) is the pore water velocity (L/T), R = 1 + is the retardation coefficient, where F' is the slope of the sorption isotherm at c (L³/M). Note that Eq. [4] is generally valid for transient flow and transport. Solution of Eq. [4] requires specification of initial and boundary conditions for c, as well as the space–time variation of the following variables: q, , r_w, R, and r_s. In the remainder of this section, each of these will be discussed.

Let us first consider the root water and solute uptake terms, r_w(z,t) and r_s(z,t). The expression for r_w(z,t) is based on Eq. [1],

[5]

where moisture and salt stress are ignored ( = 1), ET(t) (L/T) is the evapotranspiration rate (i.e., the sum of evaporation and transpiration rates), and r(z) (1/L) is the root distribution function (including evaporation from the top layers). The distribution of r(z) can take any functional form with the restriction that zr0rdz = 1, where zr is the root zone depth (L). Note that evaporation and transpiration are both represented as a distributed sink, which is convenient though not completely physically correct (e.g., Raats, 1981). Root solute uptake is considered to be coupled to root water uptake in the following way,

[6]

where a (a 0) is a parameter that may depend on the type of solute and crop species. For a = 0, Eq. [6] defines that root solute uptake is zero, whereas a = 1 corresponds to passive root solute uptake at a rate determined by the concentration in the soil solution and the rate of root water uptake. Values for a may be estimated from an expression derived by Dalton et al. (1975), which accounts for diffusive, convective, and active uptake mechanisms. A value of a > 1 would correspond to active uptake.

Solving Eq. [4] also requires calculation of the flow velocity field, v(z,t) or q(z,t), and (z,t). The flow velocity can be obtained by numerically solving the Richards equation with a sink term for root water uptake (e.g., Simunek et al., 1998). However, at large time scales (months to years to decades), one is not interested in the short-term flow dynamics due to rain or irrigation events, and a simplified description of flow may suffice. It is assumed that short-term (hourly, daily, weekly) fluctuations in flow and their effects on solute transport are not important at the observation scale t, if t is equal to or larger than a month. Kavvas (1999) refers to this phenomenon as the coarse-graining of hydrologic processes with increasing scale.

Let us first assume that flow is at a steady state at the time scale t; that is, temporal changes in water flux and moisture storage at smaller time steps are ignored. At this time scale t, flow is characterized by the effective water application ratio, p, which is defined as the fraction of surface-infiltrated water that is removed from the soil by evapotranspiration, or

[7]

where ₀ is the average infiltration rate from irrigation and rainfall during t (L/T), and is the average evapotranspiration rate during t (L/T). Since both ₀ > 0 and > 0, it follows that 0 < p < . Parameter p is defined for both downward and upward flows. In the case of upward flow, > ₀ and p becomes greater than 1. On the other hand, for downward flow, we have < ₀ and p is smaller than 1. Moreover, for downward flow, parameter p is related to the leaching fraction, LF, by

[8]

since the leaching fraction is defined as the fraction of surface-infiltrated water that drains below the root zone. An alternative way of handling upward flows is to introduce a dimensionless net upward flow rate as in Raats (1974b).

We may extend the steady-state assumption to allow for temporal changes in the infiltration rate q₀(t) and the evapotranspiration rate ET(t), while still ignoring soil moisture storage changes, and assuming that p (or LF) remains constant, so that the soil water flux, q(z,t), can be separated in space and time, or

[9]

where (z) is defined in Eq. [10] and (z) is the time-invariant soil moisture content profile. These assumptions result in a transient flux model with a "stiff" flow geometry, where changes in the boundary flux are immediately translated to velocity changes along the column (Ginn and Murphy, 1997). The physical argument for assuming a constant value for p is that temporal variations in ET rates (be they intra-annual or inter-annual) cause likewise variations in applied water rates q₀, thereby minimizing variations in p. In other words, farmers will apply more water when crop water demand increases, thereby keeping p relatively constant. It is important to note that the assumed proportionality of ET and q₀ is for average quantities over t, which doesn't require proportionality of ET and q₀ during specific irrigation and leaching events. The hypothesis of a time-invariant moisture content has been tested using flow simulations of intermittent irrigation with the Richards equation, both with and without root water uptake (Wierenga, 1977; Beese and Wierenga, 1980; Destouni, 1991). Wierenga (1977) showed that the prediction of solute transport based on a steady-state flow model is comparable to that based on a transient flow model if the breakthrough curve is presented as a function of cumulative drainage instead of time. Beese and Wierenga (1980) found even better agreement when accounting for root water uptake and equilibrium sorption processes. Similar conclusions were reached by Destouni (1991) when comparing solute travel times predicted by a transient flow model using daily boundary fluxes and a steady-state flow model using an annual average boundary flux. The time-invariant moisture profile, , may be obtained by solving the one-dimensional Richards equation with the storage term set to zero and with constant flow boundary conditions given by ₀, , and zero pressure head at the water table. Both numerical (Whisler et al., 1968) and analytical (Raats, 1974b) solutions are available. Soil water flux, on the other hand, can be directly calculated by inserting Eq. [5] and [9] into Eq. [2] and integrating. The resulting expression for (z) is,

[10]

where (z) is the normalized soil water flux at depth z (L), defined as the soil water flux divided by the soil surface infiltration rate. If p < 1, or ET(t) < q0(t), flow is downward throughout the profile, and (z) > 0 for all depths. On the other hand, if p > 1, or ET(t) < q0(t), flow is downward in the upper part of the root zone [(z) > 0], and upward in the lower part [(z) < 0]. In this case, the depth at which flow converges (zc) is given by = 0, or pzc0rd = 1. Note that since rw accounts for evaporation in addition to root water uptake, q = q0 when ET = 0.

The retardation coefficient is also needed in Eq. [4]. It is specified by assuming sorption to be linear; hence,

[11]

where K_d is the distribution coefficient (Jury et al., 1991).

Finally, initial (IC) and boundary (BC) conditions for solute concentration are specified,

[12]

[13]

where c_i(z) is the depth-dependent initial solute concentration at initial time t₀, and c₀(t) is the time-dependent solute concentration at the boundary z₀. In the case of downward flow, the boundary z₀ corresponds to the soil surface , and c₀(t) represents the concentration of the infiltrating water. In the case of capillary rise from a shallow water table, the boundary z₀ corresponds to the water table depth, and c₀(z) is the concentration of the shallow groundwater.

General Solution by Method of Characteristics
Given the assumptions above, we now proceed by solving Eq. [4] by the method of characteristics. Equation [4] is equivalent to the following system of ordinary differential equations (Garabedian, 1998)

[14]

The first equation in [14] describes the space–time trajectory, or characteristics path, of a solute parcel in differential form. The second equation in [14] describes the change of dissolved concentration, c, along the trajectory, also in differential form. The latter can be integrated analytically along the trajectory, from a starting point (z_s,t_s) in the space–time domain to a point (z,t), resulting in (Appendix B),

[15]

Equation [15] reduces to Eq. [B14] of Ginn and Murphy (1997), and Eq. [21] and [23] of Raats (1975) for the special case of zero root solute uptake and no sorption. Note that if a = 1, it follows that c = c, or the concentration of a solute parcel does not change along its trajectory.

In Eq. [15], the calculation proceeds backward in time along the trajectories until either an initial condition, c_i(z_s), or a boundary condition, c₀(t_s), is reached. This procedure is illustrated in Fig. 1, which shows solute parcel trajectories in the space–time domain for steady downward flow in the presence of root water uptake without solute uptake (p = 0.8, a = 0, exponential root distribution). In this example, the solute concentration at depth z₁ and time t₁ depends on the initial concentration at depth z_s, which is the starting depth of the trajectory through (z₁,t₁). On the other hand, the solute concentration at the same depth z₁, but at a later time t₂ depends on the boundary condition at t_s, where t_s is the starting time of the trajectory through (z₁,t₂). The trajectory through the origin (z₀,t₀) of the space–time domain is defined as the limiting characteristic (thick line in Fig. 1), and divides the space–time domain into two regions: the region where solute concentrations depend on the initial conditions (initial value problem, IVP), and the region where solute concentrations depend on the boundary conditions (boundary value problem, BVP).

View larger version (70K): [in this window] [in a new window]   Fig. 1. Solute parcel trajectories for downward flow in the presence of root water uptake: p = 0.8, a = 0, and root distribution is exponential. BVP, boundary value problem; IVP, initial value problem.

[16]

with,

[17]

[18]

where T₀(t) (T) is the time required to infiltrate the cumulative surface infiltration depth between times t₀ and t when infiltration rate is equal to the annual average rate, ₀. In other words, T₀(t) represents a rescaled time that accounts for deviations of the infiltration rate, q₀(t), from its annual average (₀). In the case that the infiltration rate remains steady, T₀ = t - t₀. Note that according to Eq. [9], T₀(t) may also be defined in terms of the evapotranspiration rate ET(t). The variable (z) (T) is defined as the travel time of a solute parcel moving from z₀ to z for time-averaged boundary conditions, ₀ and .

Equation [16] is solved first to find an expression for the limiting characteristic curve, separating the solution domain into a subdomain BVP and a subdomain IVP. The limiting characteristic is the trajectory with starting depth z_s = z₀ and starting time t_s = t₀, so that both (z_s) in Eq. [18] and T₀(t_s) in Eq. [17] are zero. It then follows from Eq. [16] that the limiting characteristic is described by,

[19]

Given a location (z,t) in the space–time domain, travel time (z) and rescaled time T₀(t_s) are calculated.

If (z) < T₀(t), then (z,t) lies on a trajectory which originated at a model domain boundary, hence its starting depth z_s = z₀. This is a boundary value problem (BVP). To solve for the solute concentration at (z,t), the starting time (t_s) must be computed (see Fig. 1). This is the time at which the solute particle that currently is at depth z was introduced at the boundary. The value of t_s is obtained from Eq. [16] by letting z_s = z₀ so that = 0. It then follows that,

[20]

where T^-1₀ is the inverse function of T₀(t).

If, however, for location (z,t), (z) > T₀(t), then (z,t) lies on a trajectory that originated at a depth within the model domain. For this initial value problem (IVP), the starting time of the trajectory through (z,t) is known, that is, t_s = t₀, but the starting depth z_s needs to be calculated (see Fig. 1). This is the depth from which the solute particle that currently is at depth z originated. The value of z_s is obtained from Eq. [16] by letting t_s = t₀, so that T₀ = 0. It then follows that,

[21]

where ^-1 is the inverse function of (z). Note that both (z) and T₀(t) are monotone increasing functions, so that their inverse exist.

Combining Eq. [15] and [19] results in the following explicit expression for dissolved solute concentration,

[22]

where H is the Heaviside function; that is, H[x] = 1 for x > 0, and H[x] = 0 otherwise, and c_i and c₀ are the initial and boundary conditions defined in Eq. [12] and [13], respectively. The first term on the right-hand side of Eq. [22] solves the IVP, whereas the second term applies to the BVP. The separation of these two terms follows from the assumption of piston displacement.

For the special case of steady-state flow, Eq. [22] reduces to (Appendix D),

[23]

where z_s = ^-1 and t_s = t - .

Solution Procedure
Summarizing, Eq. [20] through [22] represent the solution to the convective transport Eq. [4] for initial and boundary conditions given by Eq. [12] and [13], and using the assumptions presented above. The dissolved solute concentration, c(z,t), given by Eq. [22], is obtained as follows:

1. Specify boundary conditions for flow. This consists of a time series of soil surface infiltration rates, q₀(t), or evapotranspiration rates, ET(t). The average infiltration rate, ₀, or evapotranspiration rate, , is computed. Specify the average depth to the water table z_wt.
   
2. Specify initial and boundary conditions for transport.
   
3. Parameterize
   • Leaching fraction (LF), or efficiency (p),
     
   • Distribution coefficient (K_d),
     
   • Root solute uptake parameter (a),
     
   • Root distribution function, r(z),
     
   • Soil hydraulic properties.
     
   
   
4. Calculate the normalized soil water flux, (z), from Eq. [10]; if (z) > 0, flow is downward, and if (z) < 0, flow is upward.
   
5. Specify or calculate the steady-state moisture profile, as a function of soil hydraulic properties and average flow boundary conditions.
   
6. Calculate solute travel time, (z), from z₀ to z using Eq. [18], z₀ for downward flow, and z₀ = z_wt for upward flow from a shallow water table.
   
7. Calculate rescaled time T₀(t) from Eq. [17], or in case of steady-state flow, T₀ = t - t₀.
   
8. Determine whether we have a BVP or an IVP:
   • If (z) < T₀(t), we have a BVP. Calculate starting time t_s from Eq. [20], or in the case of steady-state flow, t_s = t - .
     
   • If (z) > T₀(t), we have an IVP. Calculate starting depth z_s from Eq. [21].
     
   
   
9. Calculate c(z,t) from Eq. [22], or in the case of steady-state flow, Eq. [23].
   

The model outlined above provides a robust and efficient method for determining changes in vadose zone solute concentration (e.g., soil salinity), as influenced by irrigation and drainage management (leaching fraction, irrigation water quality, water table depth), crop type (crop water use, root distribution, root solute uptake), climate, and soil type. Two useful quantities can be derived from the model. First, from an agricultural production standpoint, the solute concentration in the root zone weighted by the root distribution function, C_wm (M/L³), is of interest, since it is related to crop yield (Raats, 1974a; Hoffman and van Genuchten, 1983),

[24]

where ₀rdz = 1. Second, environmental issues mainly concern solute loadings below the root zone into the groundwater, as influenced by soil properties and soil water management decisions. Cumulative solute flux, S_cum (M/L²), is given by,

[25]

where s = qc is the convective solute flux [M/(L²T)].

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSION REFERENCES

 
Whereas the presented solutions are valid for solute transport in general, the results of the examples will be viewed in the context of soil salinity and salt transport in soils and through the deeper vadose zone. Only simulation results for the case of an exponential root distribution, r(z), will be presented. Similar calculations could easily be done for other root distributions as well, following the procedure described above. However, in the case of using the exponential root distribution function, steady-state flow can be solved using the analytical solution of Raats (1974b), which is valid for a uniform soil profile and an exponential conductivity function,

[26]

[27]

where K(z) is the hydraulic conductivity at depth z, K_s is saturated hydraulic conductivity, is a soil parameter describing the change in hydraulic conductivity as a function of pressure head h, is a parameter of the exponential root distribution function, r = 1/exp (Raats, 1974b), and z_wt is water table depth. The pressure head profile given by Eq. [27] is used to calculate the steady-state moisture content profile by specifying the soil retention curve. Here, we use the retention curve of Russo (1988),

[28]

where _r is residual moisture content, _s is moisture content at saturation, and µ is a parameter accounting for pore tortuosity and connectivity. From the moisture content profile, solute travel time is calculated using Eq. [18] by numerical quadrature (Press et al., 1990). The solute starting times (Eq. [20]) and starting depths (Eq. [21]) are calculated using a bisection algorithm (Press et al., 1990).

Further simplifications of the analytical expressions for the solute travel time and the solute starting depth may be obtained by replacing the moisture profile with a uniform moisture content, , that may be estimated independently, or calculated from the moisture profile as,

[29]

where w is a weighting function. For example, w = 1/z gives the arithmetic depth-averaged moisture content. For an exponential root distribution and a uniform moisture content , the resulting expressions for solute travel time and starting depth are given in Table 1. The expression for solute travel time when p 1 reduces to Eq. [8] of Raats (1975) for the case of no sorption and steady-state flow. Additional solutions for linear and uniform root distributions (Hoffman and van Genuchten, 1983) may be computed as well. Alternatively, depth and time distributions of water content can be computed from numerical solution of the steady-state form of Eq. [2], (Whisler et al., 1968).

View this table: [in this window] [in a new window]   Table 1. Analytical expressions for normalized soil water flux, (z), solute travel time, (z), and solute starting depth, zs, for the case of an exponential root distribution (Raats, 1974) and a uniform moisture content equal to .

View larger version (15K): [in this window] [in a new window]   Fig. 4. Influence of parameters p, a, Kd, and soil type on weighted average root zone solute concentration, Cwm. The full line corresponds to the reference case (Table 2), whereas the dashed lines represent variations from the reference values, as indicated.

View larger version (14K): [in this window] [in a new window]   Fig. 5. Influence of parameters p, a, Kd, and soil type on cumulative solute flux, Scum, at the 2-m depth. The full line corresponds to the reference case (Table 2), whereas the dashed lines represent variations from the reference values, as indicated.

View larger version (30K): [in this window] [in a new window]   Fig. 2. Solute transport for downward flow to a deep water table: triangles (steady flux) and diamonds (transient flux) represent a nonuniform moisture distribution, whereas full lines are solutions for uniform moisture profiles.

View larger version (25K): [in this window] [in a new window]   Fig. 3. Solute transport for upward flow from a shallow water table: triangles (steady flux) and diamonds (transient flux) represent a nonuniform moisture distribution, whereas full lines are solutions for uniform moisture profiles.

Upward Flow from a Shallow Water Table
Figure 3 presents similar simulation results for the case of upward flow from a shallow water table at the 2-m soil depth (z_wt = 2m). The root distribution extends to 2 m, and the moisture content increases with depth (Fig. 3a). The value of the effective water application ratio, p, is assumed larger than 1 (Table 2), representing underirrigation with part of the crop water demand supplied by capillary rise from a shallow water table. In this case, flow is downward in the upper part of the root zone [(z) > 0], and upward in the lower part of the root zone [(z) < 0]. At a depth of around 0.50 m the flow converges, which results in the solute concentration profiles of Fig. 3b. Since no leaching occurs, there is a net accumulation of salts within the root zone, leading to high solute concentrations, especially where downward and upward flows converge. The sharp peak at a depth of 0.50 m results from neglecting diffusive and dispersive forces and the assumed absence of mineral precipitation. This peak in solute concentration continues to grow as time goes beyond the 1-year simulation time. Figure 3c shows the change in time of C_wm for steady and transient infiltration rates, for both uniform (full lines) and nonuniform moisture content values (triangles or diamonds). As for downward flow, a transient infiltration rate introduces additional variation in C_wm. However, whereas C_wm attained a constant value in the downward flow case, for upward flow C_wm continues to rise after the first year due to a continued accumulation of salts in the root zone. In other words, in the upward flow case, no steady-state situation will be attained. Figure 3c also shows that the nonuniform moisture profile can be replaced by its weighted average value, _wm, to produce nearly identical results for C_m. Finally, Fig. 3d shows the change of S_cum at 2-m groundwater table depth during a 1-yr period. Negative values indicate upward salt flux by capillary rise from the water table into the root zone.

Model Parameter Sensitivity
Figures 4 and 5 show the influence of different model parameters on C_wm and S_cum at the 2-m depth. The thick lines in these figures represent reference conditions, corresponding to the parameter values listed in Table 2, while the thin lines represent solutions using the indicated parameter values while keeping the other parameters at their reference values. These simulations were done for steady infiltration and a uniform moisture content, equal to either _wm (for C_wm) or _m (for S_cum). The analytical solutions of Table 1 were used since moisture content is uniform. Figure 4 shows the sensitivity results for C_wm, with respect to the parameters p, a, and soil type. In Fig. 4a, the influence of the effective water application ratio, p, on C_wm is presented, with values ranging from 40 to 120%. As the results show, C_wm is very sensitive to p. As expected, increasing values of p result in higher root zone concentrations. For p = 1.2 ( > ₀), C_wm increases most rapidly due to the accumulation of salts in the root zone by capillary rise. In Fig. 4b, the root solute uptake parameter, a, is increased from 0 to 1. As the value of a increases, more solute is taken up by the roots and therefore root zone solute concentration decreases. For a = 1, solute is taken up at the same rate as water, and no evapoconcentration occurs. Figure 4c shows the influence of soil texture and solute sorption on C_wm. The coarse-textured soil provides for a faster breakthrough than the fine-textured soil, resulting in a faster increase in C_wm. However, the final steady-state value of C_wm is identical for both soils. The faster breakthrough in the coarse-textured soil is due to its lower moisture content, and therefore higher flow velocity, as compared with the fine-textured soil. The influence of solute sorption (K_d = 0.2 cm³ g^-1) is to slow down the system response without any effect on the steady-state value of C_wm. Finally, Fig. 5 shows the influence of various parameters on the cumulative solute flux S_cum at the 2-m depth of the water table. As was the case for C_wm, S_cum is also very sensitive to the value of p (Fig. 5a). However, the effect of increasing p is now to decrease S_cum. A larger value for p results in higher concentrations but lower water fluxes, and therefore lower solute fluxes. From Fig. 5b it can be seen that S_cum is much less sensitive to the solute uptake parameter a than is C_wm, at least for the first 3 yr. However, since breakthrough at the 2-m depth occurs in the third year, the two curves shown in Fig. 5b will increasingly diverge for t > 3 yr. Figure 5c shows the influence of soil type and solute sorption on S_cum. The results are similar to those for C_wm. Hence, for the coarse soil solute, breakthrough occurs sooner than for the fine-textured soil, and sorption introduces some delay in solute breakthrough compared with a nonsorbing solute.

Discussion
The model discussed so far assumes that the leaching fraction or effective water application ratio remains constant with time. However, when simulating over multiple years or decades, changes in land use may occur, thereby changing the leaching fraction with time. In this case, the model can be applied sequentially for each year n (major time stepping increment) with a constant LF or p value, so that the initial conditions for year n are obtained from the end concentrations of year n - 1. This was already suggested by Raats (1981) in his equations [20] and [21]. An example of this is illustrated in Fig. 6, which shows yearly increasing solute concentration profiles for a linear increase in the value of p from 0.1 in Year 1 to 0.8 in Year 8. The result is a sequence of solute waves that travel downwards through the root zone. As the value of p increases, the solute waves introduced at the surface become bigger, while displacing the earlier, smaller waves down the soil profile.

View larger version (17K): [in this window] [in a new window]   Fig. 6. Yearly solute concentration profiles for a linear increase in p from 0.1 in Year 1 to 0.8 in Year 8.

As discussed in the introduction, plant water and solute uptake may be reduced due to moisture and salt stress conditions in the root zone. This is represented conceptually by the (h,) function in Eq. [1], which takes on a functional value between 0 and 1 as determined by the corresponding matric head (h) and osmotic head () values. Model development in this paper neglected moisture and salt stress. However, the general solution of Eq. [22] remains valid under conditions of moisture stress, as long as the root water uptake distribution in Eq. [10] is modified to account for these stress conditions. Salt stress is more difficult to incorporate since it involves a two-way feedback between flow and transport; that is, the rate of water uptake in Eq. [1] now also depends on the soil solution concentration and its variations with depth and time. Salt stress can still be accounted for though by iterative calculation of the stress factor (h,) and water uptake function, using the initial salt concentrations at the start of the major timestepping increment (e.g., year). Solute concentrations are solved using Eq. [22] for the current year, which results in updated solute concentrations and values for (h,) and r_w. Subsequently, using these updated values for r_w, transport can be recalculated for that year. The solution procedure is repeated until the solution converges.

Finally, it should be noted that in addition to salt transport, the model developed in this paper may also be applied to other chemicals, such as for nitrate movement in shallow and deep vadose zones. However, the presented solution does not account for solute transformations and multiple species or precipitation–dissolution reactions.

	CONCLUSION

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSION REFERENCES

 
An efficient method is presented for calculating one-dimensional solute transport through the vadose zone in the presence of vertically distributed root water and solute uptake. An explicit analytical solution to the convective transport equation is derived using the method of characteristics, so no discretizations in space or time are required. Model parameters are related to root water and solute uptake and soil type. Two variables of practical interest are calculated from this solution, namely the weighted average root zone solute concentration and the cumulative solute flux for a specific depth below the root zone. The method can be applied for both nonuniform and uniform moisture profiles, and for both steady and transient infiltration and evapotranspiration rates, without accounting for changes in moisture storage over time. Both downward and upward flows are considered. Simulations with a nonuniform moisture profile indicate that nearly identical results are obtained by replacing the profile by a uniform moisture content equal to either the average or weighted average moisture content. The advantage of using a uniform moisture content is that closed-form analytical solutions can be derived for solute travel time. Results from the sensitivity analysis indicate that the effective application parameter p is the most sensitive to presented solutions.

The main limitations of the presented model are that the leaching fraction is required to be time invariant and that transport occurs by convection alone (i.e., hydrodynamic dispersion and molecular diffusion are neglected). However, as discussed, temporal variations in leaching fraction, for example due to changes in land use, are incorporated by discretizing the simulation period into smaller periods of constant leaching fraction. Furthermore, hydrodynamic dispersion at any spatial scale can be accounted for by spatial distribution of model parameters and/or boundary conditions, either deterministically, as in distributed modeling, or stochastically, as discussed by Simmons (1982). In the absence of dispersion, the model represents an estimate of average transport behavior.

The presented solutions may be very useful for predicting vadose zone solute transport in the presence of distributed root water and solute uptake for the long term (years to decades), by its application to local transport in both deterministic and stochastic stream tube models, specifically as related to soil salinity management for the vadose zone and groundwater.

APPENDIX A

Derivation of the Convective Transport Equation Without Dispersion or Diffusion
In this appendix, Eq. [4] is derived from the continuity equations for water (Eq. [2]) and solute (Eq. [3]). First, the derivatives in (Eq. [3]) are expanded,

[A1]

Assuming the bulk density (_b) to be constant in time, and inserting Eq. [2] into [A1] yields,

[A2]

Letting,

[A3]

where F(c) is the sorption isotherm with slope F'(c), yields, after combining Eq. [A2] and [A3],

[A4]

Rearranging leads to,

[A5]

Defining the retardation factor R as,

[A6]

we obtain the convective transport equation Eq. [4],

[A7]

where v_R = v/R = q/.

APPENDIX B

Calculation of Solute Concentration along a Trajectory
We analytically integrate the second equation of [14] to obtain the change of solute concentration along a trajectory, starting from a point (z_s,t_s) in the space–time domain to (z,t). First, we rewrite the equation,

[B1]

Inserting Eq. [5], [6], and [9] we obtain

[B2]

Then, assuming a time-invariant leaching fraction gives

[B3]

Finally, since it follows from Eq. [10] that d = prdz, Eq. [B3] reduces to

[B4]

Equation. [B4] is integrated along a trajectory from (z_s,t_s) to (z,t) in the space–time domain, yielding Eq. [15], or

[B5]

APPENDIX C

Integration of Solute Parcel Trajectories
We rewrite the left-hand side of Eq. [14] as,

[C1]

Using Eq. [9] and [11], we obtain,

[C2]

Dividing both sides of Eq. [C2] by the annual average soil surface infiltration rate, ₀, and defining the following variables,

[C3]

where T₀ (T) is the time to infiltrate the cumulative surface infiltration depth from initial time t₀ to time t at a rate equal to ₀, and (z) (T) is the solute travel time from z₀ to z for time-averaged boundary conditions, ₀ and . Equation [C2] then reduces to,

[C4]

which after integration from (z_s,t_s) to (z,t) in the space–time domain results in the integral form of a solute parcel trajectory, Eq. [16], or

[C5]

APPENDIX D

Derivation of Transport Solution for the Special Case of Steady-State Flow
Taking q₀ = q₀ = ₀ as the steady soil surface flux, and therefore T₀ = t - t₀ as defined in Eq. [17], Eq. [16] reduces to,

[D1]

The same procedure as used for transient boundary fluxes can be used to obtain expressions for the limiting characteristic, starting time t_s, and starting depth z_s. The limiting characteristic is the trajectory starting at (z₀,t₀), so its starting depth z_s = z₀ and its starting time t_s = t₀. The equation describing the limiting characteristic is then obtained from Eq. [D1] or,

[D2]

Using Eq. [D1], while substituting (z_s), yields for the BVP,

[D3]

whereas the IVP results in

[D4]

when substituting t_s = t₀ in Eq. [D1]. Combining Eq. [15], [D2], and [D3], we obtain the following transport solution for steady-state flow,

[D5]

	ACKNOWLEDGMENTS

 
The authors acknowledge financial support by the USDA Fund for Rural America, project no. 97-36200-5263, and the scientific contributions by the principal investigators Drs. W.W. Wallender, T.C. Hsiao, S.L. Ustin, R.E. Howitt, T. Harter, and G.E. Fogg. In addition, discussions with Dr. T.R. Ginn were helpful in the development of the analytical solutions of the study.

	REFERENCES

TOP ABSTRACT INTRODUCTION THEORY RESULTS AND DISCUSSION CONCLUSION REFERENCES

 

• Beese, F., and P.J. Wierenga. 1980. Solute transport through soil with adsorption and root water uptake computed with a transient and a constant-flux model. Soil Sci. Soc. Am. J. 129:245–252.
• Bresler, E. 1978. Analysis of trickle irrigation with application to design problems. Irrig. Sci. 1:3–17.
• Bresler, E., and G. Dagan. 1983. Unsaturated flow in spatially variable fields: 3. Solute transport models and their application to two fields. Water Resour. Res. 19:429–435.
• Chen, Z.Q., R.S. Govindaraju, and M.L. Kavvas. 1994. Spatial averaging of unsaturated flow equations under infiltration conditions over areally heterogeneous fields, 1. Development of models. Water Resour. Res. 30:523–533.
• Clausnitzer, V., and J.W. Hopmans. 1994. Simultaneous modeling of transient three-dimensional root growth and soil water flow. Plant Soil 164:299–314.
• Corwin, D.L., K.M. Loague, and T.R. Ellsworth. 1999. Introduction: Assessing non-point source pollution in the vadose zone with advanced information technologies. In D.L. Corwin et al. (ed.) Assessment of non-point source pollution in the vadose zone. Geophysical Monogr. American Geophysical Union, Washington, DC.
• Corwin, D.L., and R.J. Wagenet. 1996. Applications of GIS to the modeling of nonpoint source pollutants in the vadose zone, a conference overview. J. Environ. Qual. 25:403–411.[Abstract/Free Full Text]