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 Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (3) Citing Articles via Google Scholar Google Scholar Articles by Russo, D. Articles by Laufer, A. Search for Related Content PubMed Articles by Russo, D. Articles by Laufer, A. GeoRef GeoRef Citation Agricola Articles by Russo, D. Articles by Laufer, A. Related Collections Soil Salinity Spatial Variability

ORIGINAL RESEARCH

Numerical Analysis of Transport of Interacting Solutes in a Three-Dimensional Unsaturated Heterogeneous Soil

^a Dep. of Environmental Physics and Irrigation, Institute of Soils, Water and Environmental Sciences, Agricultural Research Organization, The Volcani Center, Bet Dagan 50250, Israel
^b Currently, AMEC Earth & Environmental Ltd., 160 Traders Blvd. East, Suite 110, Mississauga, ON, Canada L4Z 3K7
^* Corresponding author (vwrosd{at}volcani.agri.gov.il)

Received 4 December 2003.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY IMPLEMENTATION RESULTS DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
The present study focuses on field-scale flow and transport in three-dimensional, heterogeneous, variably saturated formations, for the case in which the flow is coupled to the transport through the dependence of the hydraulic conductivity and water retentivity on solute concentrations. Numerical simulations of flow and transport of both tracer and mixed Na–Ca solutes were employed to analyze long-term effects of the interactions between the soil solution and the soil matrix on water and solute movement, under transient, nonmonotonic flows. The simulated flows were derived from actual irrigation and weather records, along with water uptake by plant roots. Results of this study suggest that enhanced soil solution–soil matrix interactions, induced by soil alkalinity and dilution of the soil solution, may reduce both the mean and the spatial variability of the hydraulic conductivity, and, concurrently, of the velocity. Consequently, enhanced soil solution–soil matrix interactions may slow down the tracer solute movement, decrease the spreading of the tracer solute about its center of mass (particularly in the vertical direction), diminish the skewing of the tracer solute breakthrough, and decrease both the magnitude of the effective retardation factor and the rate at which it approaches its asymptotic value. In addition, our results suggest that under realistic conditions, the three-dimensionality of the flow domain, the periodicity of the rain or irrigation events, along with the substantial redistribution periods between successive events, and the spatial heterogeneity of the hydraulic properties of the variably saturated formation may compensate in part for the adverse effects of soil alkalinity on flow and transport on the field scale. This last finding has practical implications regarding the use of sewage water for irrigation.

Abbreviations: BTC, breakthrough curve • CDE, convection–dispersion equation • CEC, cation-exchange capacity • CP, control plane • DDL, diffuse double layer • PDF, probability density function • PSD, pore-size distribution • RSF, random space function • SAR, Na adsorption ratio

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY IMPLEMENTATION RESULTS DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
INCREASING AMOUNTS OF potentially hazardous chemicals produced by various industries and by agricultural operations are entering the ecosystem. Quantitative descriptions of chemical transport in the vadose (unsaturated) zone on the field scale are essential for the management of these chemicals in the ecosystem to minimize the potential of the drainage water to pollute underlying groundwater supplies. A distinct feature of the porous formation on the field scale is the spatial heterogeneity of those of its properties that affect transport. This spatial heterogeneity is generally irregular, and it occurs on a scale not captured by laboratory samples. These features have a distinct effect on the spatial distribution of solute concentration that results from transport through heterogeneous porous formations, as has been observed in field experiments (e.g., Schulin et al., 1987; Butters et al., 1989; Ellsworth et al., 1991; Roth et al., 1991; Flury et al., 1994; Forrer et al., 1999) and demonstrated by simulation (e.g., Russo, 1991; Russo et al., 1994, 1998; Tseng and Jury, 1994) of solute transport in unsaturated, heterogeneous soils.

The aforementioned studies focused on transport of a single noninteracting solute. In reality, however, mixed-salt solutions, with ions that interact with the soil matrix, play an important role in many field problems, such as those involving the simultaneous transport of water and mixed Na–Ca salts. In such a case, in addition to convection and dispersion, the transport process is governed also by physicochemical interactions between the soil solution and the soil matrix, including cation exchange, anion exclusion, and changes in the soil pore-size distribution (PSD). The latter could affect the formation conductivity and retentivity. Because the magnitude of soil solution–soil matrix interactions depends on water content and on the concentration and composition of the soil solution, the flow is coupled to the transport.

Russo (1988a)(1988b) analyzed transport of mixed Na–Ca solutions in a homogeneous, saturated–unsaturated soil during one-dimensional vertical infiltration, taking into account soil matrix–soil solution interactions. In this approach, these interactions were quantified using theoretical considerations based on the mixed-ion diffuse double layer (DDL) theory, the structure of the clay particles, the soil's PSD, and hydrodynamic principles (Bresler, 1972, 1978; Russo and Bresler, 1977). Results suggested that in a homogeneous soil, the effect of soil solution–soil matrix interactions on the transport of water and solutes might be significant and that it would generally increase as the soil water content and the soil solution Na adsorption ratio (SAR) increased, as the soil solution concentration decreased, as the clay fraction of the soil increased, and as the soil texture became finer. In a layered soil (Russo, 1988b) the extent of the effect of these interactions on the transport process depends on the sequence of the soil layering and on the shape of the initial soil solution's composition and concentration profile.

The approach of Russo (1988a) was extended to the field scale (Russo, 1989a, 1989b) by treating the heterogeneous soil as though it were composed of a series of isolated stream tubes (vertically homogeneous, independent soil column) characterized by differing velocities. The analysis concentrated on transport of mixed Na–Ca salts under conditions of one-dimensional, time-invariant, vertical infiltration. Results of the analyses suggested that the field-average movement of both the water and the solutes might be retarded because of the soil solution–soil matrix interactions and that their spatial variability might be increased relative to the case in which these interactions had not been considered.

It should be emphasized, however, that the conceptualization of the soil spatial heterogeneity adopted by Russo (1989a)(1989b) is quite restricting and may not describe real field soils (Russo, 1995, 1997). A more realistic approach to the analysis of field-scale transport of mixed Na–Ca solutions, therefore, must take into account the three-dimensionality of the spatial variations of the soil properties and, concurrently, of the flow-controlled attributes.

The purpose of our study was twofold. First, we analyzed the transport of mixed Na–Ca salts in a realistic, three-dimensional, spatially heterogeneous, variably saturated soil, considering the effect of the soil solution–soil matrix interactions on the simultaneous movement of water and solutes. Second, we assessed the long-term effect of the soil solution–soil matrix interactions on the soil's capability to transfer water and solutes.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY IMPLEMENTATION RESULTS DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
Simulation of Transport in Heterogeneous, Variably Saturated Soil
Our purpose here is to simulate, on the field scale, the transport of waterborne mixed Na–Ca salts in a variably saturated, heterogeneous soil. In this analysis, ion precipitation or dissolution is not taken into account; that is, the simple but practical soil–water system containing only Na, Ca and Cl ions is considered here. We consider a three-dimensional flow domain that spans several meters in both horizontal and vertical directions and analyze the transport problem by means of physically based flow and transport models and a stochastic presentation of the soil properties that affect water flow and solute transport. The approach undertaken here is a "single realization" approach (e.g., Ababou, 1988; Russo, 1991; Polmann et al., 1991; Russo et al., 1994, 1998; Tseng and Jury, 1994) that may be understood as quantifying the macroscopic spreading of a single plume for particular site-specific applications. This situation should be distinguished from that of an ensemble macrodispersive flux that accounts for differences among the trajectories of random concentration replicates over the ensemble. When the solute inlet zone is sufficiently large to ensure ergodicity, the discrepancy between single-realization statistics and ensemble mean statistics might be reduced significantly (Dagan, 1989).

Governing Partial Differential Equations
We assume here that the water flow is described locally by the Richards equation, the physical parameters of which are visualized as realizations of stationary random space functions (RSFs). It is further assumed that the transport of each of the ions of the mixed Na–Ca salt is described locally by a convection–dispersion equation (CDE). When the ions in the soil solution interact with the soil matrix, the properties of the heterogeneous soil which affect the water and solute transport can be expressed as functions of the spatial coordinate vector, x, the pressure head, , the sum of the molar concentration of the cations, C₀ = c_Na + c_Ca (mol L^–1), and the ratio between the molar concentration of mono- and divalent cations, R = c_Na/c_Ca, in the soil solution.

In a Cartesian coordinate system (x₁,x₂,x₃), where x₁ (L) is directed vertically downwards, assuming local isotropy, considering water uptake by the plant roots and neglecting the effect of osmotic gradients on the flow, the "mixed" form of the Richards equation, governing saturated–unsaturated flow in a nonrigid, three-dimensional, heterogeneous soil is

[1]

where t (T) is time, = (x,t) (L) is the pressure head, = (,x) and the scalar K = K(,x) (L T^–1) are the volumetric water content and the hydraulic conductivity, respectively, of the soil in a rigid "stable" reference state in which the changes in the soil PSD due to soil solution–soil matrix interactions are negligibly small (i.e., when R 0, and C₀ ); ^r = ^r (C₀,R,,x) and K^r = K^r(C₀,R,,x) are relative water content and relative conductivity, respectively, expressing the effect of the soil solution–soil matrix interactions on the soil PSD, and S_w = S_w(x,t) (T^–1) is a sink term representing water uptake by plant roots, given (Nimah and Hanks, 1973a, 1973b) by

[2]

where R_e (L^–2) is the root effectiveness function, _r (L) is the total pressure head at the root–soil interface, and (L) is the osmotic pressure head of the soil solution.

Similarly, neglecting ion precipitation or dissolution, the equation governing the transport of mixed Na–Ca salts is

[3]

where m = 1 to 3 indicates the Cl, Na, and Ca ions, respectively; c_m (mol L^–1) is the resident solute concentration of the mth ion, expressed as mass per unit volume of soil solution; ^*_m is an "effective" volumetric water content, given (Russo, 1988a) by

[4a]

for the Cl (m = 1), and by

[4b]

for the Na (m = 2) or for the Ca (m = 3) cations, respectively; _ex and _adm are the exclusion and the adsorption volume fractions, respectively, given (Russo, 1988a) by

[5a]

and

[5b]

where A_ex (L² m^–1), A_ad (L² m^–1), and _b (M L^–3) are the specific exclusion surface area, the specific adsorption surface area, and the bulk density of the soil, respectively, and ^–_ex (mol L^–2) and ⁺_ad (mol L^–2) are the quantity of monovalent anionic charges repelled by, and the quantity of monovalent cationic charges adsorbed to a unit area of the solid surface, respectively. u_i (i = 1,2,3) (L T^–1) is the pore water velocity vector, and D_ij (i,j = 1,2,3) (L T^–2) is the pore-scale dispersion tensor. When molecular diffusion is small enough to be excluded from the analysis, D_ij is given (Bear, 1972) as

[6]

where _L and _T (L) are the longitudinal and the transverse pore-scale dispersivities; _ij is the Kronecker delta, and |u| = ^1/2.

Characterization of Flow and Transport Parameters
The local-scale soil properties that are here hypothesized to control the field-scale spreading of the mixed Na–Ca salts are the soil hydraulic functions, K() and () at a "stable" reference state in which the changes in the soil PSD caused by soil solution–soil matrix interactions are negligibly small (i.e., when R 0, and C₀ ) and the local soil coupling-interaction parameters K^r, ^r, _ex, _adm (m = 2,3), and the longitudinal, _L, and the transverse, _T, components of the pore-scale dispersivity tensor.

For a given soil material whose clay fraction is mostly montmorillonite, a procedure to estimate the local soil coupling-interaction parameters, all as functions of , R, and C₀, was outlined by Russo (1988a). The inherent soil data required are the soil hydraulic functions, K() and (), at the "stable" reference state, the soil cation-exchange capacity (CEC), the soil specific surface area, A_s, and the soil bulk density, _b. These soil properties vary in space in an irregular fashion, and the length scales of these variations are much greater than the pore scale. To quantify the spatial variations in these properties, we regard the soil as a continuum, the properties of which are continuous functions of the spatial coordinates. The elementary volume representing the soil properties of the field is characterized by a macroscopic scale that is small in comparison with the scale of the field heterogeneities, but contains many pores. Our main concern is the effect of the spatial variations in the hydraulic properties and in the soil coupling-interaction parameters on the transport of mixed Na–Ca salts. Therefore, deterministic, constant values were adopted for pore-scale dispersivities: _L = 2 x 10^–3 m and _T = 1 x 10^–4 m (Perkins and Johnston, 1963).

To simplify the analysis, it is assumed here that the local K() and () relationships in the stable reference state, are described by the van Genuchten (1980) parametric expressions. Ignoring local hysteresis and local anisotropy, and considering the pressure head, , as the dependent variable, they read as follows:

[7a]

[7b]

where = ( – _ir)/(_s – _ir), is the effective water saturation; _s and _ir are the saturated and the irreducible water contents, respectively; K_s is the saturated conductivity; ß and m are parameters related to the soil PSD, and n = 1/(1 – m).

It is assumed further that each of the parameters of Eq. [7], CEC, A_s, and _b, denoted by p(x), are second-order stationary, statistically anisotropic RSF, characterized completely by a constant mean, <p(x)>, independent of the spatial position, and a covariance, C_pp(x,x'), that, in turn, depends on the separation vector, = x-x', and not on x and x' individually. A three-dimensional, exponential covariance is adopted here for p(x); that is,

[8]

where ' = (x-x')/I_p is the scaled separation vector, ' = |'|; ²_p and I_p = (I_p1,I_p2,I_p3) are the variance and the correlation length scales of p(x), respectively.

Furthermore, axisymmetric anisotropy is adopted for C_pp(); that is, I_pv = I_p1 and I_ph = I_p2 = I_p3 are the characteristic length scales of p(x) in the vertical (longitudinal) direction, and in the horizontal plane, respectively. Statistical parameters of the soil parameters of Eq. [7] were adopted from the Bet Dagan trench (Russo and Bouton, 1992) and are summarized in Tables 1 and 2. The soil in this site is a Hamra Red Mediterranean soil (Rhodoxeralf) whose clay fraction is dominated by montmorillonite. Inasmuch the sampling scheme for CEC, A_s, and _b in this site (Naor, 1987) was not adequate to evaluate the spatial structure of these soil properties, the correlation scales of logK_s were adopted for CEC, A_s, and _b.

View this table: [in this window] [in a new window]   Table 1. Estimates of the statistical parameters of the spatially variable soil properties logKs, logß, n, s, ir, CEC, As, and b.

We considered a flow domain measuring 10 m in each of the horizontal directions and 2.5 m in the vertical direction. The flow domain was discretized to a large number of uniform cells, measuring 0.05 and 0.20 m in the vertical and horizontal directions, respectively. Using the covariance function (Eq. [8]) and the values of ²_p, I_ph, and I_pv (Table 1) the turning bands method (Tompson et al., 1989) was used to generate realizations of logK_s, logß, n, _s, _ir, CEC, A_s, and _b for each of the 50 by 50 by 50 cells of the flow domain, taking into account the cross-correlation coefficients (Table 2) between the various soil properties.

Using Eq. [7], two-dimensional tables of K_r = K/K_s and _r = ( – _ir)/(_s – _ir) were constructed as functions of ß' = ß|| and n. Employing a bilinear interpolation scheme, these tables were used to calculate values of K_r and = _r(_s – _ir) + _s for each cell of the flow domain by means of the generated realizations of ß, n, _s, and _ir and the simulated water pressure head, . For each of the cells of the flow domain, for which the specific soil surface area, A_s, cation-exchange capacity, CEC, bulk density, _b, soil water retention, (), and the hydraulic conductivity, K(), in the reference "stable" state are known, three-dimensional tables of the soil coupling-interaction parameters, K^r, ^r, _ex, and _adm (m = 2,3) were constructed as functions of , R, and C₀ by using the procedure outlined by Russo (1988a). Employing a trilinear interpolation scheme, these tables were used to calculate values of K^r, ^r, _ex, and _adm (m = 2,3) for each cell of the flow domain, by means of the simulated , R, and C₀. The estimated cumulative distribution plots of the resultant relative conductivity, K^r, relative water content, ^r, elution factor for Cl, E_f = 1 – _ex/, and the retardation factor for Na, R_f = 1 + _ad2/, for selected values of the Cl^– concentration, C = 10^–3C₀[1 + (R + 1)^–1] (mM), the SAR = c_Na/%c_Ca = R[10³C₀(R + 1)^–1]^1/2 (mM^1/2), and the pressure head, , are depicted in Fig. 1 .

View larger version (28K): [in this window] [in a new window]   Fig. 1. Cumulative distribution plots of relative conductivity, Kr, relative water content, r, the elution factor for Cl, Ef, and the retardation factor for Na, Rf, for different Cl– concentrations, C, for Na adsorption ratio, SAR = 10, and pressure head, = –2 m.

Boundary and Initial Conditions
For the analyses presented here, we considered the movement of both tracer and Na–Ca solutes in a flow domain, , consisting of an initially tracer-free soil (c_i0 = 0), with prescribed uniform initial concentrations, c_im, m = 1,2,3, and pressure head, _i. Water flow is originated from a time-dependent, periodic potential influx, F(t_j), imposed on the entire soil surface at x₁ = 0. The influx is given by F(t_j) = RI, t^'_j < t_j < t^''_j, F(t_j) = 0, elsewhere, where RI > 0 is a spatially uniform rain or irrigation intensity, j = 1 to N_I, and N_I is the number of the irrigation or rain events. Water may leave the flow system through the lower boundary (x₁ = L₁) by gravity flow. There is no water flow across the vertical planes located at x₂ = 0 and x₂ = L₂ and x₃ = 0 and x₃ = L₃, for which the normal derivatives of vanish.

For the transport, we considered a case in which during each of the irrigation or rain events, Na–Ca salt with concentrations, c_0m, m = 1 to 3, respectively, invades the soil through the entire soil surface, which is orthogonal to the vertical mean flow direction. A tracer solute with concentration, c₀₀, invades the soil via a planar source of dimensions (₂₂ – ₂₁) x (₃₂ – ₃₁) located at the soil surface, during time, t₀ = t^''_k – t^'_k, where k = 1 to N_P, and N_P is the number of the pulse application events. There is no solute transport across the soil surface outside the source, nor across the vertical planes located at x₂ = 0, x₂ = L₂, x₃ = 0 and x₃ = L₃, for which the normal derivatives of c vanish. A zero-gradient-boundary is specified at the bottom (x₁ = L₁). The appropriate boundary and initial conditions for the flow and the transport are given in the Appendix.

	IMPLEMENTATION

TOP ABSTRACT INTRODUCTION THEORY IMPLEMENTATION RESULTS DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
We applied Eq. [1] and [3], subject to Eq. [A1] through [A5] and [A6] through [A12], respectively, to investigate the transport of a tracer solute (bromide) and a Na–Ca salt by a transient, unsaturated water flow. Note that in the case of the tracer solute, the transport is described by Eq. [3] with m = 0 and ^*_m = . The "mixed" form of Richards' equation governing three-dimensional water flow (Eq. [1]) was approximated by a fully implicit Euler scheme with a truncation error of O. The resulting system of nonlinear algebraic equations with respect to the pressure head, , was solved iteratively, by means of the so-called modified Picard method (Celia et al., 1990). The resulting system of linear algebraic equations was solved by line successive overrelaxation (LSOR) method (e.g., Hageman and Young, 1981).

The set of the CDEs (Eq. [3]) governing transport of Na–Ca salt and a tracer solute in a three-dimensional flow domain was approximated by an operator-splitting approach. A numerical solution of the advection equation was obtained by a second-order accurate, explicit finite-difference scheme proposed by Zaidel and Levi (1980), while dispersive fluxes were approximated by a standard central difference scheme. For more details see the appendix in Russo et al. (1994).

In the case of the Na–Ca salt, assuming that the soil CEC is time-invariant and that the CEC is equal to the sum of the exchangeable cations, an iterative scheme was used to couple the m sets of Eq. [3] (m = 1,2,3), subject to two requirements: (i) cation mass balance and (ii) charge balance. In this iterative scheme, for each time step, appropriate exchangeable cation-solution cation adjustments are made to ensure that both cation mass balance and charge balance are preserved—the sum of exchangeable and solution concentration remains constant for each cation, and solution cation charges equal the solution anion charges. The iterative procedure continues until a prescribed error tolerance is met, or until the number of iterations has exceeded a predetermined limit.

Note that the transport equations are strongly coupled to the flow equation because of the dependence on Darcy velocities. The flow equation is itself coupled to the transport equations through the dependence of K() and () on C₀ = c_Na + c_Ca and R = c_Na/c_Ca because of the soil solution–soil matrix interactions. As this coupling is generally mild, a serial solution method for the flow and the transport equations, with a suitable procedure to account for the dependence of K() and () on C₀ and R, provides a feasible and efficient computational procedure. Before the solution of the flow Eq. [1], an estimate of the concentrations, c_m(x,t) (m = 1,2,3) is needed. Using the indices i,j,k and n for the space and the time grid points, respectively, a first estimate of cⁿ⁺¹_m,i,j,k (m = 1,2,3) is obtained (Bresler and Laufer, 1974) from

[9]

Thus, for given values of ⁿ_r,i,j,k, ⁿ_r,i,j,k and cⁿ_m,i,j,k and using the tables of K_r(,R,C₀), _r(,R,C₀), _ex(,R,C₀) and _adm(,R,C₀) and Eq. [9], initial estimates of K^n+1/2_r,i,j,k, ^n+1/2_r,i,j,k, _ex^n+1/2_r,i,j,k and _adm^n+1/2_r,i,j,k are obtained. These, in turn, are used to obtain initial estimates of ⁿ⁺¹_i,j,k and ⁿ⁺¹_i,j,k (Eq. [1]) and cⁿ⁺¹_m,i,j,k (Eq. [3]). By using the initial estimates of ⁿ⁺¹_i,j,k, ⁿ⁺¹_i,j,k and cⁿ⁺¹_m,i,j,k, improved estimates of K^n+1/2_r,i,j,k, ^n+1/2_r,i,j,k, _ex^n+1/2_r,i,j,k and _adm^n+1/2_r,i,j,k are obtained, and, these, in turn, are used again in the solution of Eq. [1] and [3]. The iterative procedure continues until a prescribed error tolerance is met, or until the number of iterations has exceeded a predetermined limit.

Crop, climatic conditions, and agricultural practice typical of the central area of the coastal region of Israel were considered in these simulations to obtain realistic water uptake and surface boundary conditions for the flow. A citrus orchard with a complete cover of the soil surface and a given, time-invariant root effectiveness function (Mantell and Goell, 1972), was considered here. Irrigation dates, rates, and amounts were those used in the agricultural practice in this region, with a total amount of 700 mm yr^–1. Meteorological data from the weather station at Bet Dagan were used to estimate potential evapotranspiration rates, ET_p(t), as well as to provide rain rates and amounts during the winter (typical of Israel), with a total amount of 580 mm yr^–1. The estimated cumulative potential evapotranspiration was 1005 mm yr^–1.

Water uptake by plants was implemented by a maximization iterative approach (Neuman et al., 1975). In this approach, the rate of transpiration per unit area of soil surface at time t, given by T_r(t) = *S_w(x,t)dx, is maximized subject to two requirements. First, the actual transpiration rate is not allowed to exceed the potential rate, T_p. Second, the total pressure head at the root–soil interface, _r, is not allowed to fall below a critical value, _c, equivalent to the so-called wilting point of the soil–plant system (e.g., _c = –150 m). To simplify the problem, soil evaporation was set to zero (i.e., ET_p = T_p). Cumulative seasonal amounts of applied water and potential evapotranspiration used in this study are given elsewhere (Russo et al., 1998).

The initial pressure head, _i, was selected as the pressure head corresponding to <K(,x)>/<K_s(x)> = 0.001, i.e., _i = –3 m, while the rain and the irrigation intensities were selected as RI = 0.75 = 0.25 m d^–1, respectively. To match the behavior of actual rain and irrigation events, both the duration of such an event, t^''_j – t^'_j, and the time interval between successive events, t^'_j – t^''_j–1, were taken as time-dependent, ranging from 0.2 to 0.5 d and from 7 to 14 d, respectively. Water flow and transport of the tracer solute and the Na–Ca salt were simulated for L₁ = 2.5 m, L₂ = 10 m, L₃ = 10 m. The flow scenario started at the beginning of the irrigation season (April 1) and proceeded for four successive years.

For the tracer solute with initial concentration, c_i0 = 0, and inlet concentration, c₀₀ = 10 mol L^–1, the dimensions of the inlet zone are given by ₂₁ = ₃₁ = 3 m and ₂₂ = ₃₂ = 8 m. Two pulses (N_P = 2) of the tracer were applied during the second irrigations in the first year, t^'_j = 7 d, and in the third year, t^'₂ = 727 d; in both applications, the duration of the pulse was the same, t₀ = t^''_k – t^'_k = 0.05 d, k = 1 to N_P. In the case of the Na–Ca salt, during the irrigation season, the transport scenario describes the displacement of a relatively low-Na soil solution (SAR_i = 1.22) by a relatively high-Na solution (SAR₀ = 12.7). In this case the initial concentrations are c_im = 10, 2.5, and 3.75 mol L^–1, and the inlet concentrations are c_0m = 10, 9, and 0.5 mol L^–1, for m = 1, 2, and 3, respectively. During the rainy season, which starts after 240 d and lasts for 90 d, the inlet concentrations are c_0m = 1.5, 0.4, and 0.55 mol L^–1, for m = 1, 2, and 3, respectively (SAR₀ = 0.55).

	RESULTS

TOP ABSTRACT INTRODUCTION THEORY IMPLEMENTATION RESULTS DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
Movement, Spreading, and Breakthrough of the Interacting Solutes
Profiles of the horizontal averages, the mean values, <p(x₁;t)>, and the standard deviations, [<p'(x₁;t)²>]^1/2, where p = C or SAR, are depicted in Fig. 2 and 3 . The profiles in these figures were calculated from the simulated soil Cl^– concentration, C(x,t), and the soil solution Na adsorption ratio, SAR(x,t), obtained from the solution of Eq. [1] subject to Eq. [A1] through [A5], and that of Eq. [3] subject to Eq. [A6] through [A12]. The profiles in Fig. 2 and 3 represent different intervals between two successive rain or irrigation events at the beginning (a,e) and at the end (b,f) of the second rainy season and at the beginning (c,g) and at the end (d,h) of the third irrigation season. It is shown in Fig. 2 that during the rainy seasons (Fig. 3a,b,e,f), the Cl^– concentration, C, is reduced throughout the soil profile, particularly in the upper part. The decrease in the mean Cl^– concentration is associated with a decrease in its standard deviation at depths shallower than the position of the mean Cl^– front, F_Cl, (see Eq. [13] below) and vice versa at depths deeper than the position of mean F_Cl. The reverse tendency occurs during the irrigation seasons (Fig. 2c,d,g,h).

View larger version (34K): [in this window] [in a new window]   Fig. 2. Profiles of the horizontal averages, the mean values (a,b,c,d), and the standard deviations (e,f,g,h) of the Cl– concentration. The profiles in this figure represent different elapsed times (0, 1, 2, 3, 4, 6, 8, and 10 d) following the cessation of a rain or irrigation event, (a,e) at the beginning and (b,f) at the end of the second rain season and (c,g) at the beginning and (d,h) at the end of the third irrigation season.

View larger version (28K): [in this window] [in a new window]   Fig. 3. Profiles of the horizontal averages, (a,b,c,d) the mean values, and (e,f,g,h) the standard deviations of the Na adsorption ratio, SAR. The profiles in this figure represent different elapsed times (0, 1, 2, 3, 4, 6, 8, and 10 d) following the cessation of a rain or irrigation event, (a,e) at the beginning and (b,f) at the end of the second rainy season and (c,g) at the beginning and (d,h) at the end of the third irrigation season.

During the irrigation season, the relatively low-Na, diluted soil solution in the upper part of the soil profile is displaced by a relatively high-Na, concentrated solution. In other words, the preferentially adsorbed Ca is exchanged by Na with less affinity for the exchange complex. In this case, the exchange isotherm is "unfavorable"; that is, Q_Na/c_Na (where Q_Na = _adNac_Na) increases with increasing c_Na. Consequently, the mean SAR profiles (Fig. 3c,d) are smeared, while the SD profiles (Fig. 3g,h) are highly skewed: their peaks decrease and are displaced slightly downwards with increasing time. Note that during the irrigation seasons, the increase in the mean SAR (Fig. 3c,d) is associated with a decrease in the SD of the SAR (Fig. 3g,h) at depths shallower than that of the mean SAR front, F_SAR, (see Eq. [13] below), and vice versa at depths deeper than that of the mean F_SAR. Figures 2 and 3 suggest that the soil solution–soil matrix interactions are intensified during the rain events (soil leaching) and moderated during the irrigation events (soil salinization).

Solute breakthrough for the Na–Ca salt is calculated from the flux-averaged concentration, c_f, of the respective ions, defined (Kreft and Zuber, 1978) by c_f = s· vdA/u·vdA, where s and u are the solute flux and the velocity vectors, respectively, with s = s·v and u = u·v being the mass of solute and the volume of water per unit time moving through a surface element of unit area and water content, , normal to v, and the integration is over a planar area A. Mean and standard deviation of the solute breakthrough curve (BTC), evaluated by averaging c_f over horizontal control planes (CPs) located at three different vertical distances, L, from the inlet zone at the soil surface, are depicted in Fig. 4 for the Cl^– (a,d) Na (b,d), and Ca (c,f) ions.

View larger version (47K): [in this window] [in a new window]   Fig. 4. (a,b,c) Means and (d,e,f) standard deviations of the flux-averaged concentrations of Cl–, Na, and Ca, crossing a control plane (CP) located at three different vertical distances, L, from the injection zone, as functions of time, t.

The SDs of the BTCs (Fig. 4d–f) generally follow the same pattern as the mean BTCs. The former, however, particularly those that were evaluated at the shallower CP, are much more sensitive to the temporal, transient flow conditions at the soil surface than the latter. The uncertainty in the BTCs, expressed in terms of the coefficient of variation, CV[BTC] = SD[BTC]/E[BTC], decreases with increasing soil depth. Below the root zone, CV[BTC] is essentially insensitive to the cyclic recharge at the soil surface.

The seasonal effects of the various ions on the BTCs are clearly shown in Fig. 4. Over the 4-yr period, at soil depths below the root zone, a slight soil salinization occurs during the irrigation seasons (Fig. 4a), while the alkalinization of the soil occurs at a rate that decreases with increasing time and depth (Fig. 4b).

Response of the Concentration-Dependent Flow System
Profiles of the mean, <p(x₁;t)> and the standard deviation, [<p'(x₁;t)²>]^1/2 of p(x), where p = , y = logK, J_i = (– – x₁)/x_i and u_i (i = 1,2) are depicted in Fig. 5 through 8 . The profiles in these figures represent different intervals between successive rain or irrigation events at the end of the second rain season (a,c) and at end of the third irrigation season (b,d). Note that in Fig. 7 and 8 only the transverse components of J and u in the direction parallel to the x₂ axis, J₂, and u₂, respectively, are shown; their counterparts, J₃ and u₃ (not shown here) exhibit the same pattern. It is demonstrated in Fig. 5 through 8 that water application influences the water regime mostly in the upper part of the soil profile. Note that because of the temporal dependency of the events that occur at the soil–atmosphere interface, the amounts and application rates of the applied water, and, concurrently, the distances to the wetting fronts associated with the different seasons (Fig. 5) differ from each other.

View larger version (29K): [in this window] [in a new window]   Fig. 5. Profiles of the horizontal averages, (a,b) the mean and (c,d) the standard deviation of the pressure head, , (c,d) at the end of the second rain season, and (b,d) at the end of the third irrigation season. The profiles in this figure represent different elapsed times (0, 1, 2, 3, 4, 6, 8, and 10 d) following the cessation of a rain or irrigation event.

View larger version (37K): [in this window] [in a new window]   Fig. 8. Profiles of the horizontal averages, (a,b) the mean and (c,d) the standard deviation of the longitudinal and the transverse (parallel to the x2 axis) components of the velocity vector, u, (a,c at the end of the second rain season and (b,d) at the end of the third irrigation season. The profiles in this figure represent different elapsed times (0, 1, 2, 3, 4, 6, 8, and 10 d) following the cessation of a rain or irrigation event.

View larger version (27K): [in this window] [in a new window]   Fig. 6. Profiles of the horizontal averages, (a,b) the means and (c,d) the standard deviations of y = logK (a,c) at the end of the second rainy season and (b,d) at the end of the third irrigation season. The profiles in this figure represent different elapsed times (0, 1, 2, 3, 4, 6, 8, and 10 d) following the cessation of a rain or irrigation event.

View larger version (38K): [in this window] [in a new window]   Fig. 7. Profiles of the horizontal averages, (a,b) the means and (c,d) the standard deviations of the longitudinal and the transverse (parallel to the x2 axis) components of the head-gradient vector, J, (a,c) at the end of the second rainy season and (b,d) at the end of the third irrigation season. The profiles in this figure represent different elapsed times (0, 1, 2, 3, 4, 6, 8, and 10 d) following the cessation of a rain or irrigation event.

In response to a rain or irrigation event, the mean values of , logK, J₁, and the absolute deviation of the mean J_i (i = 2,3) from zero increase, while their standard deviations decrease; on the other hand, the mean value of u₁, the absolute deviation of the mean u_i, (i = 2,3) from zero and the standard deviation of u_i (i = 1,2,3) increase. The reverse tendency occurs during the redistribution period following the termination of the rain or irrigation event.

The mean values of the transverse components of the head gradient, J, and the velocity, u, vectors (Fig. 7 and 8, respectively) are essentially time-independent. On the other hand, the longitudinal component of the mean head gradient vector, <J₁>, is very sensitive to the temporal variations of the influx at the soil surface. In response to the water application, <J₁> increases considerably in the upper part of the soil profile (Fig. 7). It decreases with increasing time following the termination of the rain or irrigation event and may change its sign in the upper part of the soil profile as the profile dries. At a sufficiently large vertical distance from the soil surface; however, <J₁> fluctuates about unity.

Consider the response of the concentration-dependent flow system during the rainy season. As is evident from Fig. 2 and 3, after the first few rain events the displacement of the relatively high-Na, concentrated soil solution by a relatively low-Na, diluted solution is confined to the upper part of the soil profile. After a series of rain events the leached zone is extended to a deeper soil depth, but the decrease in the soil solution SAR is still confined to the few centimeters near the soil surface. The combination of low-concentration, high-Na soil solution results in an increase in the magnitude of the physicochemical interactions between the soil solution and the soil matrix, and, concurrently, in a rearrangement of the soil PSD. The latter results in a considerable decrease in both the mean and the standard deviation of logK (Fig. 6a,c). The decrease in logK, however, is compensated for by considerable increases in the mean value of J₁ (associated with a slight decrease in its SD), in the absolute deviation of the mean J₂ (and J₃, not shown here) from zero, and in their standard deviations (Fig. 7a,c). Consequently, the decrease in the mean value of u₁ and in the standard deviations of u_i, (i = 1,2,3) and the increase in the absolute deviation of the mean u₂ (and u₃, not shown here) from zero are moderate (Fig. 8a,c).

The reverse tendency occurs during the irrigation season. The relatively low-Na, diluted soil solution in the upper part of the soil profile is displaced by a relatively high-Na, concentrated solution. After the first few irrigation events, the increase in both soil solution concentration and SAR is confined to the upper part of the soil. After a series of irrigation events, however, the salinized–alkalinized zone is extended to a deeper soil depth (Fig. 2 and 3). Consequently, during the irrigation season, the increases in the mean and the standard deviation of logK (Fig. 6b,d) are accompanied by decreases in the mean value of J₁ and in the absolute deviation of the mean values of J₂ and J₃ from zero, and in their standard deviations (Fig. 7b,d). Consequently, the mean value of u₁ and the standard deviations of u_i, (i = 1,2,3) increase, while the absolute deviation of the mean u₂ and u₃ from zero decrease during the irrigation season (Fig. 8b,d).

Movement, Spreading, and Breakthrough of the Tracer Solute
To analyze the movement and spreading of the tracer solute, we focus our attention on the moments of the spatial distribution of the point values of the tracer resident concentration, c(x,t), given (Dagan, 1989) as

[10a]

[10b]

[10c]

where M is the total mass of the solute, R(t) = (R₁,R₂,R₃) is the coordinate of the centroid of the solute plume, and S^'_ij (i,j = 1,2,3) are second spatial moments, proportional to the moments of inertia of the plume.

The total mass of the solute, M, the principal components of R, and the principal components of S_ij = S^'_ij – S^'_ij (i,j = 1,2,3) are depicted in Fig. 9 as functions of time, for the first and the third years. Note that both R₁ and S₁₁ fluctuate with time; these fluctuations stem from the temporal and spatial variations of the unsaturated flow. The time fluctuations in R₁, however, tend to smooth out as time passes and the center of mass of the plume is displaced to a sufficiently large distance from the inlet zone.

View larger version (26K): [in this window] [in a new window]   Fig. 9. (a) Total mass, M, of the tracer solute, (b) the principal components of the coordinate location, R, of the tracer center of mass, and (c) the principal components of the spatial covariance tensor, Sij, as functions of time. Results are depicted for the first and the third years.

The decrease in R₁ in the third year compared with the first year reflects the decrease in the tracer solute velocity, V₁ = dR₁(t)/dt, induced by the reduction in the hydraulic conductivity caused by the rearrangement of the soil PSD. Similarly, the changes in the structure of S₁₁ with travel time reflect the changes of the concentration distribution that occur because of the heterogeneity induced in the tracer solute velocity field by small-scale heterogeneity in the soil hydraulic properties. The decreases in both the magnitude of S₁₁ and the temporal fluctuations of S₁₁ in the third year compared with the first year (Fig. 9c) therefore suggest a reduction in the small-scale heterogeneity of the soil hydraulic properties because of the alkalinization process.

The means and the standard deviations of the tracer solute BTCs, evaluated by averaging the tracer flux concentration, c_f, over CPs located at three different vertical distances L from the inlet zone at the soil surface, are depicted in Fig. 10 for the first and the third years. In general, the mean BTCs, E[BTC] (Fig. 10a), exhibit skewness that is characterized by an earlier breakthrough, by an enhanced peak arrival, and by a longer tailing. However, this skewness decreases with increasing distance to the CP. The SD of the BTCs, SD[BTC] (Fig. 10b), also exhibit skewness, with longer tailings than those of the respective E[BTC] (Fig. 10a). Figure 10 suggests that the uncertainty in the mean BTCs, expressed in terms of CV[BTC] = SD[BTC]/E[BTC], is relatively large around the first arrival time, decreases with increasing t, reaches a minimum in the vicinity of the E[BTC] peaks, and increases with increasing t as t continues to increase, since E[BTC] 0 faster than SD[BTC]. Comparison of the means and standard deviations of the tracer BTCs associated with the first and the third years suggests that soil alkalinity, which was shown to moderate the variability in the water velocity (Fig. 8) and concurrently in the tracer solute spreading (Fig. 9c), is shown to decrease the skewing of the tracer pulse also.

View larger version (28K): [in this window] [in a new window]   Fig. 10. (a) Means and (b) standard deviations of the flux-averaged concentration of the tracer solute crossing a control plane (CP) located at three different vertical distances, L, from the injection zone, as functions of time, t. Results are depicted for the first and the third years.

[11]

where T₁ and T₂ are the first two normalized travel time moments of the horizontally averaged flux concentration, given as

[12]

The normalized PDFs for one-particle travel time, g₁L/V₁, calculated from the simulated flux concentrations and the fitted PDFs are depicted in Fig. 11 as functions of T = R₁(t)/L, for CPs located at three different vertical distances, L, from the inlet zone at the soil surface, and for the first and the third years. Note that for both years, the fitted PDFs at the shallow CP located within the root zone cannot reproduce the skewing in the simulated travel time PDFs. For CPs located at greater depths, however, the goodness of fit of the Fickian PDFs (with parameters optimized at the shallow CP) to the simulated PDFs improves, particularly in the third year.

View larger version (30K): [in this window] [in a new window]   Fig. 11. Normalized probability density functions for one-particle travel time (g1L/V1) for the tracer solute, as a function of scaled travel time, T = R1(t)/L, for control planes located three different vertical distances, L, from the injection zone. R1(t) is the longitudinal component of the coordinate location of the tracer center of mass and t is time. Results are depicted for (a) the first and (b) the third years.

Exchange Front and Effective Retardation Factor
For the Na–Ca salt, the vertical position of the front of each of the ions, or of the exchange (expressed in terms of the SAR) front, at a given point in the horizontal x₂x₃ plane, may be defined as

[13]

where = c_Cl, c_Na, c_Ca or SAR = c_Na/%c_Ca.

Horizontal averages of the exchange front, F_SAR, the mean, E[F_SAR], and the standard deviation, SD[F_SAR], are depicted in Fig. 12a as functions of travel time. Both the mean and the standard deviation of the SAR front essentially grow linearly with increasing travel time. Note that the rate of advance of the mean position of the exchange front (Fig. 12a) is much slower than the rate of displacement of the center of mass of the tracer in the longitudinal direction (Fig. 9b). This stems from the highly nonlinear Na–Ca exchange isotherm. During the alkalinization process, because the preferentially adsorbed Ca is substituted by Na that has less affinity for the exchange complex, the adsorption isotherm is "unfavorable"; that is, both the mean and the variance of the slope of the exchange isotherm increase with increasing c_Na. This enhances the effect of the retention of the solute in the exchange phase in retarding the rate of advance of the mean exchange front.

View larger version (19K): [in this window] [in a new window]   Fig. 12. (a) Horizontal averages, the mean and the standard deviation of the position of the SAR front (Eq. [13]), and (b) the effective retardation factor as functions of elapsed time.