Search

Treated wastewater (TWW) is an important water resource, especially in semiarid and arid regions. However, there are concerns that irrigation with TWW could lead to degradation of soil physical and hydraulic properties. The objective of our study was to determine the effects of long-term (≥ 15 years) irrigation with secondary TWW on some basic and hydraulic soil properties of a clay soil. Undisturbed soil samples (cores) were taken to a depth of 4.5 m (in sections of 0.5 m) over a diagonal cross section of a five year old orchard irrigated with TWW. Samples were taken from five sites within the tree rows (i.e., representing soil directly affected by TWW; referred to as “within rows”) and four sites between the rows of trees (i.e., the control treatment representing soil that was not directly subjected to the irrigation water; referred to as “between rows”). Soil analyses included an array of basic properties, determination of a continuous particle size distribution and measurement of the saturated hydraulic conductivity (HC). The latter two were used for the computation of soil characteristic curve, Θ(ψ), and the unsaturated HC curve, K(ψ). Similar bulk density, moisture content, cation exchange capacity, pH and exchangeable sodium percentage (ESP) levels were observed for the TWW irrigated samples and the control ones. However, irrigation with TWW caused a significant reduction in the saturated HC, Ks. The computed Θ(ψ) curve at a given soil depth, averaged over the different sites, was similar for the TWW-irrigated samples and the control ones. On the contrary, the computed K(ψ) curve at a given soil depth, averaged over the different sites, for the TWW-irrigated samples were lower than those for the control samples at matric potential ≳− 100 cm (= pF ≲ 2); similar K(ψ) values were noted at pF > 2 for the two treatments. The observed differences in the hydraulic properties between the TWW-irrigated samples and the control ones in this specific matric potential range, albeit the similarity in their ESP, suggest that long term irrigation with TWW affected structural porosity via narrowing macro- and mesopores (> 70 and 30–70 μm, respectively). It is further postulated that this adverse impact of irrigation with TWW on structural porosity might be associated with previously reported effects of TWW on the composition of the dissolved organic matter in the soil solution.

Detailed numerical analyses of flow and transport were used to investigate the effect of spatially connected features on transport in three-dimensional (3-D), spatially heterogeneous, combined vadose zone-groundwater flow systems. Formations with spatially connected fine-textured and coarse-textured features (F-formation and C-formation, respectively), representing the10th and the 90th percentiles of the distributions of the formation's hydraulic parameters, respectively, were considered here. Results of the analyses suggest that in steady state flow, when the unsaturated zone of the combined flow domains is relatively wet, as compared with a Multivariate-Gaussian (MG) formation, spatially connected features may reduce the solute first arrival time, particularly in the C-formation, and may enhance the spreading of the solute breakthrough, particularly in the F-formation. The effect of the spatially connected features on the hydrological response, however, decreases as the unsaturated zone becomes drier. The latter result stems from the decrease in the fraction of the water-filled, pore-space occupied by the connected structures, with decreasing water content. The latter finding also explains the result that the response of more realistic, combined flow systems, whose unsaturated zone is associated with relatively low, intermittent water contents, is essentially independent of the spatially connected features of the formations, regardless of their soil texture.

We study the standard deviation of water saturation SDS as function of the mean saturation (Formula presented.) by a stochastic model of unsaturated flow, which is based on the first-order solution of the three-dimensional Richards equation. The model assumes spatially variable soil properties, following a given geostatistical description, and it explicitly accounts for the different scales involved in the determination of the spatial properties of saturation: the extent L, i.e., the domain size, the spacing Δ among measurements, and the dimension (Formula presented.) associated to the sampling measurement. It is found that the interplay between those scales and the correlation scale I of the hydraulic properties rules the spatial variability of saturation. A “scale effect” manifests for small to intermediate L/I, for which SDS increase with the extent L. This nonergodic effect depends on the structural and hydraulic parameters as well as the scales of the problem, and it is consistent with a similar effect found in field experiments. In turn, the influence of the scale (Formula presented.) is to decrease the saturation variability and increase its spatial correlation. Although the solution focuses on the medium heterogeneity as the main driver for the spatial variability of saturation, neglecting other important components, it explicitly links the spatial variation of saturation to the hydraulic properties of the soil, their spatial variability, and the sampling schemes; it can provide a useful tool to assess the impact of scales on the saturation variability, also in view of the several applications that involve the saturation variability.

The yield of avocado trees planted on clayey soils decreases due to irrigation with treated-waste water (TWW). We hypothesized that the main cause for this yield reduction is the reduction in water uptake by the trees roots. The aim of this numerical study was to identify the main soil factors that control the reduction in water uptake by the trees roots, and to test various soil substrate-based and water-based management schemes design to counterweigh the water uptake reduction. The study relies on physically based, three-dimensional (3-D) simulations of flow and transport in variably saturated, spatially heterogeneous, flow domain, conducted for three successive years. The main findings of this study suggest that: (i) the long-term effect of irrigation with TWW on the response of the flow system is attribute to the salinity of the TWW, and not to its sodium adsorption ratio, SAR; (ii) with respect to improving water uptake by the trees' roots, the water-based scheme that alternates irrigation water quality between TWW and desalinized water, DSW, (ADW) performed better than the water-based scheme that uses fresh water only (FW). The soil substrate-based schemes, TUFp, that used trenches with highly coarse-textured soil material and pulse irrigations, and, particularly, SAop, that used trenches with finer soil texture, performed substantially better than the soil substrate-based scheme that used trenches with highly coarse-textured soil material only (TUF); (iii) with respect to minimizing solute leaching below the root zone, the water-based schemes, FW, and, particularly, ADW, performed substantially better than the soil substrate-based schemes.

Pollution of groundwater by nitrate originating from irrigated fields was considered for this study. We hypothesized that under cropped conditions, low-salinity irrigation water (e.g., desalinated water) could reduce nitrate leaching below the root zone, due to two possible mechanisms: (i) decreased vertical water fluxes and (ii) increased nitrogen uptake by plant roots due to chloride&ndash, nitrate competition. The main goal of this study was to investigate this hypothesis. Considering a citrus grove, the investigation relied on three-dimensional (3-D) simulations of flow and transport in a variably saturated and spatially heterogeneous flow domain performed for three successive years. Results of the analyses suggest that the main mechanism responsible for the reduction in the nitrate leached below the root zone under irrigation with low-salinity water is the effect of the latter on the spatial distribution of the rate of water uptake by the roots. The latter, in turn, significantly reduces water content, hydraulic conductivity, and vertical velocity, and, consequently, solute mass fluxes along the soil profile. On the other hand, chloride&ndash, nitrate interaction has only a relatively small effect on the nitrate mass fluxes at relatively deep soil depths, far below the root zone, particularly when the irrigation water salinity decreases.

Nitrate contamination of groundwater under land used for intensive-agriculture is probably the most worrisome agro-hydrological sustainability problem worldwide. Vadose-zone samples from 0 to 9 m depth under citrus orchards overlying an unconfined aquifer were analyzed for variables controlling water flow and the fate and transport of nitrogen fertilizers. Steady-state estimates of water and NO3-N fluxes to groundwater were found to vary spatially in the ranges of 90-330 mm yr(-1) and 50-220 kg ha(-1) yr(-1), respectively. Calibration of transient models to two selected vadose-zone profiles required limiting the concentration of NO3-N in the solution that is taken up by the roots to 30 mg L(-1). Results of an independent lysimeter experiment showed a similar nitrogen-uptake regime. Simulations of past conditions revealed a significant correlation between NO3-N flux to groundwater and the previous year's precipitation. Simulations of different nitrogen-application rates showed that using half of the nitrogen fertilizer added to the irrigation water by farmers would reduce average NO3-N flux to groundwater by 70%, decrease root nitrogen uptake by 20% and reduce the average pore water NO3-N concentration in the deep vadose zone to below the Israeli drinking water standard; hence this rate of nitrogen application was found to be agro-hydrologically sustainable. Beyond the investigation of nitrate fluxes to groundwater under citrus orchards and the interesting case-study aspects, this work demonstrates a methodology that enables skillful decisions concerning joint sustainability of both the water resource and agricultural production in a common environmental setting.

First-order analysis was used to analyze the probability density function (PDF) of the logarithm of conductivity, log K , in variably saturated, bimodal, heterogeneous formations, and to investigate the implications of the results with respect to solute transport. The bimodal formations were viewed as mixtures of two populations (background soil and embedded soil) of differing spatial structures. Two distinct cases were considered: in the first case, the texture of the embedded soil is finer than that of the background soil, while the second case is the reverse situation. Results of this study suggest that, because for a given volume fraction of the embedded soil, both the mean and the variance of log K depend on the mean pressure head, H , in a manner that depends on the texture of the embedded soil relative to that of the background soil, so does the log K PDF. Furthermore, because of the inherent concave behavior of the log K variance–pressure head relationships in these formations, the shape of the bimodal log K PDF may change with increasing H , such that at a certain critical mean pressure head it may degenerate into an “equivalent,” unimodal PDF. One of the main findings of this study suggests that even when the two subdomains of the formation are characterized by mild heterogeneity, a relatively small volume fraction of coarse-textured embedded soil may lead to a highly skewed log K PDF, which, in turn, may exhibit an exceedingly large tail associated with the small K values, particularly when the contrast between mean conductivities of the two subdomains of the formation at saturation is relatively large and when the formation is relatively dry.

One of the alternative solutions to the problems associated with the agricultural use of stony desert soils is to apply imported soil material into trenches that are aligned along the crop rows. The purpose of the present study was twofold: first, to analyze the effect of two imported soil materials available in the Arava Valley of Israel on the movement and spread of water and solute originating from multiple trickle line laterals, taking into account the trench geometry, the discharge of the trickle line laterals, the pattern of the plant root distribution, and the spatial heterogeneity of the hydraulic properties of the local soil; and, second, to extend the analyses to hypothetical imported soil materials characterized by different sets of hydraulic parameters (saturated conductivity, K s , and the van Genuchten soil parameters α and n , related to the soil9s pore size distribution). Results of the present study suggest that in the case of relatively low-conductive imported soil materials, a wide and shallow trench may leach the solute below a horizontal control plane more efficiently than a narrow and deep trench, while the opposite is true for high-conductive imported soil materials. Combinations of the threshold values of K s and α above which a narrow and deep trench is more efficient in solute leaching than a wide and shallow trench calculated here suggest that, for a given n , a relatively large threshold value of K s may compensate for a relatively small threshold value of α, particularly when n is relatively large.

First-order analysis, based on a stochastic continuum presentation of flow and a general Lagrangian description of transport, was used to investigate the effects of several characteristics of a bimodal, spatially heterogeneous, variably saturated formation, on solute breakthrough curves, under steady-state, gravity-dominated unsaturated flow conditions. The bimodal formation was viewed as a mixture of two populations (background soil and embedded soil) of differing hydraulic properties and differing spatial structures. Results of the present analysis, restricted to relatively small inclusions' volume fraction, suggest that as in saturated flow, for a formation of given statistics and mean water saturation, both the spreading of the mean solute breakthrough curve and the uncertainty in its prediction may be appreciable, especially in bimodal formations in which: (i) the contrast between the mean properties of the two subdomains is relatively large, (ii) the inclusions' volume fraction is relatively large, and (iii) the characteristic length scales of the inclusions are relatively large as compared with those of the heterogeneity of the background soil. Results of this study, unique to unsaturated flow conditions, stem from the inherent concave nature of the log-conductivity varianceimean pressure head relationships in bimodal, variably saturated formations. Consequently, under unsaturated flow conditions, both the spreading of the mean solute discharge and the uncertainty in its prediction are concave functions of mean saturation. Starting with saturated formation, they initially decrease with decreasing mean saturation, reach their minima, and then increase as mean saturation further decreases, exceeding their counterparts in saturated flow conditions. Implications of the results with respect to the problem of groundwater contamination where risk assessment is required, are briefly discussed.

The purpose of the present study was to analyze movement and spreading of water and a passive solute (chloride) on a hillslope under surface drip irrigation, taking into account the texture of the soil, the slope of the terrain, the spatial heterogeneity of the soil hydraulic properties, and water extraction by plant roots. Results of the present investigation suggest that under surface drip irrigation, the movement and spread of water and solutes are affected mostly by the soil texture and less by the terrain slope. Increasing terrain slope is shown to increase the deflection of the trajectories of the centroid of the solute mass from the vertical axis, particularly in fine-textured soils associated with low saturated conductivity and considerable capillary forces, in which the interaction between adjacent drip line laterals may be appreciable. Transient flows originating from periodic water application and water uptake by plant roots are shown to enhance the effect of the terrain slope on water flow and solute movement, particularly in fine-textured soils. Implications regarding the problem of sensor placement with respect to drip irrigation management are briefly discussed.

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.

The transport of a pulse of trifluralin (2,6-dinitro-N,N-dipropyl-4-[trifluoromethyl] benzenamime) applied via a subsurfacedripper was analyzed numerically. Results of the analyses suggest that the movement and spread of trifluralin in the soil is considerably retarded by its strong adsorption to the solid phase of the soil. This is particularly so in soils which contain a considerable fraction of organic C and in fine-textured (clayey) soils with low hydraulic conductivity and high water retentivity. Water uptake by plant roots and the resultant rapid decrease of water velocity with increasing distance from the dripper restricts further the downward movement of trifluralin and its potential to pollute the groundwater. The presence of dissolved organic matter (DOM) in the irrigation water may enhance both the movement and the spread of trifluralin in the soil, particularly in coarse-textured soils with a relatively small fraction of organic C. Because of the strong adsorption of trifluralin to the soil, its concentration in the aqueous phase of the soil is very low and it decreases further with increasing time because of degradation and nonequilibrium sorption. This is particularly so in coarse-textured soils with a relatively large fraction of organic C. Nevertheless, results of the analyses suggest that for soils of quite widely differing textures and organic C contents, a trifluralin concentration of c t = 10 -6 kg m -3 may persist in the vicinity of the dripper for a relatively long period of time (90 d) even with relatively small applied mass (3 x 10 -5 kg).

Numerical simulations of water flow and solute transport were used to investigate the effect of the recharge (rain/irrigation) at the soil surface on solute spreading and breakthrough in a realistic, three-dimensional, heterogeneous, vadose zone–groundwater flow system. Results of the analyses suggest that smaller recharge at the soil surface increases the relative variability in the response of the unsaturated flow system, decreases the groundwater recharge from the unsaturated zone, and, concurrently, decreases the degree of nonuniformity of the groundwater flow. Consequently, smaller recharge at the soil surface enhances the longitudinal and, especially the transverse spreading of the solute plume and increases the spreading of the expected solute BTC at the lower boundary of the unsaturated zone. Furthermore, smaller recharge at the soil surface decreases the extent of the vertical penetration of the solute plume into the groundwater, slows the peak arrival of the expected BTC, and considerably increases the spreading of the expected solute BTC at a given vertical control plane in the saturated region. It was shown that for the combined flow system considered here in which the solute plume is displaced to a sufficiently large vertical distance from the inlet zone at the soil surface, the variable that controls the transport is the cumulative recharge at the soil surface, and not the recharge rate.

Many details about the flow of water in soils in a hillslope are unknowable given current technologies. One way of learning about the bulk effects of water velocity distributions on hillslopes is through the use of tracers. However, this paper will demonstrate that the interpretation of tracer information needs to become more sophisticated. The paper reviews, and complements with mathematical arguments and specific examples, theory and practice of the distribution(s) of the times water particles injected through rainfall spend traveling through a catchment up to a control section (i.e., "catchment" travel times). The relevance of the work is perceived to lie in the importance of the characterization of travel time distributions as fundamental descriptors of catchment water storage, flow pathway heterogeneity, sources of water in a catchment, and the chemistry of water flows through the control section. The paper aims to correct some common misconceptions used in analyses of travel time distributions. In particular, it stresses the conceptual and practical differences between the travel time distribution conditional on a given injection time (needed for rainfall-runoff transformations) and that conditional on a given sampling time at the outlet (as provided by isotopic dating techniques or tracer measurements), jointly with the differences of both with the residence time distributions of water particles in storage within the catchment at any time. These differences are defined precisely here, either through the results of different models or theoretically by using an extension of a classic theorem of dynamic controls. Specifically, we address different model results to highlight the features of travel times seen from different assumptions, in this case, exact solutions to a lumped model and numerical solutions of the 3-D flow and transport equations in variably saturated, physically heterogeneous catchment domains. Our results stress the individual characters of the relevant distributions and their general nonstationarity yielding their legitimate interchange only in very particular conditions rarely achieved in the field. We also briefly discuss the impact of oversimple assumptions commonly used in analyses of tracer data.

Kinetically controlled volatilization and dissolution of nonaqueous phase liquids (NAPLs) may play an important role in the transport of volatile compounds in the unsaturated (vadose) zone. In this study, some one- and two-dimensional steady state transport problems are solved analytically. The one-dimensional case is pertinent to pollution by a relatively long, mainly horizontally spread leak of NAPL. The two-dimensional case corresponds to situations in which the pollution spreads primarily vertically, originating at the ground surface and migrating to the top of the capillary fringe, and in which the solution domain may be represented by a cross-sectional model. Solutions of the steady state transport problems are used to investigate effects of several parameters, characterizing the advective-dispersive and purely diffusive transport regimes, on the NAPL concentration distribution for the one- and two-dimensional cases, respectively. Results of this analysis indicate that the mass exchange between NAPL and other phases may not reach equilibrium, even for relatively large mass transfer rate coefficients and small water infiltration rates, if this zone has relatively small vertical or horizontal extent and is located close to the fully open ground surface. Analysis of local volatilization and dissolution fluxes shows that, under equilibrium conditions, the main losses of the organic phase take place at the upper part of the NAPL zone.

Simulations of transient, unsaturated vertical flow in different initially dry homogeneous soils using Richards' equation suggest that, for a relatively coarse numerical grid, there is no universal weighting scheme to estimate interblock conductivities capable of giving an accurate resolution of the resulting steep wetting fronts, except for the scheme associated with the Kirchhoff transformation (KT). The principal drawback of KT is that it can hardly be applied to heterogeneous soils. A new asymptotic weighting (AW) scheme, based on the asymptotic behavior of the hydraulic conductivity function in the vicinity of residual saturation, is proposed. Numerical tests demonstrate that the AW scheme describes wetting fronts in homogeneous soils nearly as accurately as a scheme based on KT. Applicability and robustness of the AW scheme for simulating water infiltration into heterogeneous soils having substantial vertical variation of hydraulic properties are also demonstrated. Extension of the proposed scheme to higher dimensions is straightforward.

Core scale estimates of soil parameters of the Gardner-Russo and van Gemachten models of the hydraulic conductivity and water retention functions were obtained for 417 undisturbed soil cores taken from a wall of a trench (20 m long, 2.5 m deep), using a procedure based on inverse problem methodology. These estimates were used to evaluate the first two statistical moments of the underlying random space functions (RSFs), using the restricted maximum likelihood estimation procedure, coupled with the weighted least squares procedure, to estimate parameters of models of the covariance and the drift functions of the pertinent RSFs. The fitted models were used to evaluate the mean and covariance functions of the hydraulic conductivity and water rentention functions for given water saturations. Covariance functions of log-saturated conductivity (log Ks) and the “shape” parameters of the Gardner-Russo and van Genuchten models exhibited statistical anisotropy characterized by aspect ratios that vary between 3 to 4. Correlation scales of log Ks were larger than those of the “shape” parameters. Consequently, the product of the variance of log unsaturated conductivity and its correlation scale remained essentially invariant for a considerable range of water saturation. The implications of these results regarding stochastic modeling of transport in heterogeneous porous formations and possible applications of the results of this study are discussed briefly.

[1] The use of saline water for irrigation is based on the concept that increasing irrigation water quantity may offset in part the negative effects associated with higher salinity. The main goals of this paper were (1) to investigate the impact of the controllable variables, the relative quantity (Qr), and the quality (C0) of the irrigation water on water uptake, active root volume, and solute leaching under daily drip irrigation in fine (clay) and coarse (sandy) textured soils; (2) to test the hypothesis that the reduced active root zone associated with salt accumulation may further reduce water uptake by plant roots; and (3) to assess consequences of using (C0-Qr) substitutions based on the FAO recommendations on salinity control with respect to agricultural and environmental aspects. Novel findings of this study suggest that in the clay soil mean water extraction rate, Sw is quite robust to C0-Qr substitutions, and the decrease in transpiration, T, induced by increasing C0 is mainly due to the decrease in the soil volume providing most (90%) of the root water uptake, Vs. In contrast, in the sandy soil, Vs may increase with C0-Qr substitutions, and the decrease in T is mainly due to the decrease in Sw, induced by osmotic pressure head reduction. Consequently, in the clay soil, a simultaneous increase of Qr and C0 based on the FAO substitution curve may lead to an increase in root uptake efficiency, T/Vs, while it may cause the opposite effect in the sandy soil. Furthermore, compensation for the adverse effects of increasing C0 on relative transpiration, Tr, by C0-Qr substitutions is more efficient in the clay soil than in the sandy soil. The C0-Qr substitution approach may lead to an increase in the salt load leached below the root zone, thus increasing the potential for polluting groundwater resources. Results of the present study suggest that for both soils, from both agricultural and, particularly, environmental points of view, irrigation with high-quality water which may comply with economical criteria is desirable.

[1] A stochastic model is developed for the analysis of the traveltime distribution fτ in a hillslope. The latter is described as made up from a surficial soil underlain by a less permeable subsoil or bedrock. The heterogeneous hydraulic conductivity K is described as a stationary random space function, and the model is based on the Lagrangian representation of transport. A first-order approach in the log conductivity variance is adopted in order to get closed form solutions for the principal statistical moments of the traveltime. Our analysis indicates that the soil is mainly responsible for the early branch of fτ, i.e., the rapid release of solute which preferentially moves through the upper soil. The early branch of fτ is a power law, with exponent variable between −1 and −0.5; the behavior is mainly determined by unsaturated transport. The subsoil response is slower than that of the soil. The subsoil is mainly responsible for the tail of fτ, which in many cases resembles the classic linear reservoir model. The resulting shape for fτ is similar to the Gamma distribution. Analysis of the fτ moments indicates that the mean traveltime is weakly dependent on the hillslope size. The traveltime variance is ruled by the distribution of distances of the injected solute from the river; the effect is coined as source zone dispersion. The spreading due to the K heterogeneity is less important and obscured by source zone dispersion. The model is tested against the numerical simulation of Fiori and Russo (2008) with reasonably good agreement, with no fitting procedure.

[1] First-order analysis was used to investigate the effect of a few characteristics of the vadose zone and the solute source at the soil surface on the expected breakthrough curves (BTCs) of tracer solutes in a combined unsaturated-saturated spatially heterogeneous formation with a three-dimensional, statistically anisotropic structure. The analysis is based on a stochastic continuum representation of the stationary and nonstationary Eulerian velocity vectors in the unsaturated and saturated domains, respectively, and a general Lagrangian description of the solute mass flux in steady state flow. Results of the present analysis suggest that the length scales and the flow characteristics of the spatially heterogeneous, unsaturated zone that affect the travel time probability density function between the soil surface and the water table also affect the expected solute BTC at a control plane located downstream in the groundwater. The flow in the unsaturated zone might enhance the early arrival time and the spreading of the expected solute BTC in the groundwater, especially when the vertical extent of the unsaturated zone is relatively large; when the unsaturated zone consists of relatively dry, coarse-textured, stratified soil material; and when the horizontal extent of the solute source in the direction of the mean groundwater flow is relatively large. The effect of the characteristics of the vadose zone and the solute source at the soil surface on the expected solute BTC in the groundwater decreases as the distance to the control plane increases and as the process of the solute breakthrough is more advanced. The relatively good agreement between the results of the presentfirst-order,stochasticanalysiswiththeresultsofdetailed,numericalanalysesofflow and transport in a three-dimensional heterogeneous, combined vadose zone–groundwater flow domain suggests that the first-order stochastic analyses may be sufficiently reliable to showappropriatetrendswhichmightoccurinrealisticcombinedvadosezone–groundwater flow systems. Citation: Russo, D., and A. Fiori (2009), Stochastic analysis of transport in a combined heterogeneous vadose zone – groundwater flow system, Water Resour. Res., 45, W03426, doi:10.1029/2008WR007157.

The problem of transport of a conservative nonreactive solute in a vertical cross section of a hypothetical partially saturated, scale-heterogeneous soil under transient water flow was analyzed here. It was assumed that locally the water flow and the solute transport can be described by the Richards' equation and by the one-component convection dispersion equation, respectively. The simulated water content and the solute concentration distributions in the vertical cross section of the soil at different elapsed times were quantified in terms of space averages and two-point autocorrelation functions. The time evolution of the solute plume was quantified in terms of its first two normalized spatial moments, from which the time dependence of the longitudinal and the transverse components of the solute velocity vector, and the spatial covariance tensor, were estimated. The results of this study, which are relevant to solute transport at the local or the plume scale, demonstrated the considerable variability in the solute concentration in space and time, due to the complex heterogeneity of the soil hydraulic properties in both the vertical and the horizontal directions. Consequently, the movement of the solute plume was characterized by a compression-expansion phenomenon, attributed to the decrease in the effective solute velocity through the zones of relatively fine-textured soil material. It was concluded that existing stochastic vadose zone transport models may be applicable to shallow depths but may fail to describe the actual spread of a solute plume when the transport takes place at relatively large depths, mainly because of the neglect of the significant vertical heterogeneity in the soil hydraulic properties.

Although volatile organic compounds located in buried waste repositories or distributed through the unsaturated soil zone have the potential to migrate to the atmosphere by vapor diffusion, little attention has been paid in the past to estimating the importance of volatilization losses. In this paper a screening model is introduced which evaluates the relative volatilization losses of a number of organic compounds under standard soil conditions. The model is an analytic solution to the problem wherein the organic chemical is located at time zero at uniform concentration in a finite layer of soil covered by a layer of soil devoid of chemical. The compound is assumed to move by vapor or liquid diffusion and by mass flow under the influence of steady upward or zero water flow while undergoing first-order degradation and linear equilibrium adsorption. Loss to the atmosphere is governed by vapor diffusion through a stagnant air boundary layer. Calculations are performed on 35 organic compounds in two model soils with properties characteristic of sandy and clayey soil. The model identifies those compounds with high potential for loss during 1 year after incorporation under 100 cm of soil cover and also is used to calculate the minimum soil cover thickness required to reduce volatilization losses to insignificant levels during the lifetime of the compound in the soil. From the latter calculation it was determined that certain compounds may volatilize from deep subsurface locations or even groundwater unless the soil surface is sealed to prevent gas migration.

[1] Field-scale solute transport in a three-dimensional, heterogeneous, variably saturated soil, originating from multiple planar sources (MS) is analyzed and compared with its counterpart originating from a single planar source (SS). The case under consideration is a citrus grove planted on a Hamra Red Mediterranean soil (Rhodoxeralf) in the central part of the coastal region of Israel, with a distinct rainy period during the winter and irrigations during the rest of the year. Results of the analyses show that for both the MS and the SS cases, solute transport during the irrigation season is characterized by a restricted downward movement and spread and by a considerable increase in concentration, while the opposite situation occurs during the rain season. In addition, results of the analyses suggest that as compared with the SS case, the MS case is characterized by a larger uncertainty in the concentration point values, by a slower solute convection and a smaller longitudinal solute spread, by a larger transverse solute spread, and by larger skewness and uncertainty in the solute breakthrough curve (BTC). Regarding the solute sampling design problem, our findings suggest that for both cases a pair of sampling points located in the horizontal plane of the field may be sufficient in order to provide relatively accurate estimates of characteristics of the transport (i.e., displacement and spread of the solute plume in the direction of the mean flow and mean solute BTC) for the entire spatially heterogeneous domain. To achieve these desirable results, in the MS case the first sampling point must be located within the wetted area in the vicinity of one of the planar sources, while in the SS case these results are independent of the location of the first sampling point. For both cases, however, an effort to quantify the uncertainty in the mean solute BTC requires few additional sampling points.

Six boreholes were drilled during the course of a year to a depth of 2 m beneath the water table, located at a depth of about 28 m, under agricultural land sprinkler irrigated with treated sewage effluents in the Coastal Plain aquifer of Israel to determine the extent of penetration of 20 pharmaceuticals and personal care products (PPCPs) into the unsaturated zone. The fields were planted to turf and had different histories of effluent irrigation. From each borehole, 7 to 21 samples were taken for analysis of PPCPs, as was the underlying groundwater. Nine PPCPs (carbamazepine and its metabolite 10-hydroxy-10,11-dihydrocarbamazepine, acridone and acridine, venlafaxine, sulfamethoxazole, oxcarbazepine, O -desmethylvenlafaxine, and caffeine) were detected in the vadose zone of the study area to a depth of 27 m. For example, the detected concentrations of carbamazepine were up to 109 ng/kg, of caffeine up to 36,700 ng/kg, and of venlafaxine up to 50 ng/kg. Only five of the compounds (carbamazepine, acridone, venlafaxine, sulfamethoxazole, and caffeine) were found in the underlying groundwater with concentrations in the nanogram per liter range. The results of this work show that significant amounts of PPCPs can penetrate even a thick vadose zone of 27 m with sections containing up to 50% clay and up to 0.40% soil organic C. Venlafaxine, for example, penetrated the vadose zone at an average velocity of 2.8 to 4 m/yr. Irrigation with treated sewage effluents or contaminated water should be carefully considered if the penetration of PPCPs into groundwater is undesirable.

[1] First-order analysis, based on a stochastic continuum presentation of the flow and a general Lagrangian description of the transport, was used to investigate the effects of a few characteristics of a bimodal, heterogeneous, partially saturated formation, on the principal components of the time-dependent macrodispersion coefficient, under steady state, unsaturated flow conditions. Results of the present analysis show that generally, a bimodal distribution of the formation properties, characterized by a relatively complex spatial correlation structure, contributes to the variability in water velocity and, consequently, may considerably enhance solute spreading. This applies especially in formations in which (1) the correlation length scales and the variances of the soil properties associated with the embedded soil are much larger than those of the background soil, (2) the contrast between mean properties of the two subdomains is large, (3) mean water saturation is relatively small, and (4) the volume fraction of the flow domain occupied by the embedded soil is relatively large. One of the major findings of this study is the important role played by the formation parameter α under unsaturated flow conditions. Consequently, the macrodispersion in a bimodal, unsaturated formation may be larger than its counterpart in a saturated formation, especially when the degree of saturation decreases, when the volume fraction of the flow domain occupied by the embedded soil increases, when the disparity between the correlation length scales of the two subdomains is relatively small, and when the texture of the embedded soil is coarser than that of the background soil.

In a heterogeneous field in which the soil water properties vary, the yield of a crop may also differ from point to point. Using actual field variability data of irrigated peanut (Arachis hypogaea L.) yield, of soil water properties, and stochastic approaches, effects of spatial variability of soil water on the crop yield distribution and on the spatial correlation structures are evaluated. Values of total dry matter yield (YT) [together with two additional crop yield components pods (Yₚ) and hay (Yₕ)], soil water content [before (θb) and after (θₐ) irrigation], and saturated conductivity, Kₛ, were used to estimate distribution types (normal or lognormal) and their moments. The two-point autocorrelation functions of YT, Yₚ, Yₕ, θₐ, θb, and Kₛ, were calculated and used to estimate their integral scales (J). The calculated J-values for Kₛ, θₐ and θb, and YT are 25, 21, and 13 m, respectively. These results and calculations of the cross correlations between crop yield and soil water components were used to establish an approximate relationship between the crop yield and the soil water variables regarded as random spatial variables.

This review presents and examines relevant information from existing spatial variability studies of soil water and solute transport properties. Although most of the information available allowed only a conventional statistical analysis (mean and variance) of the pertinent properties, the field studies of (Nielsen, Biggar, and Erh (1973)) and (Russo and Bresler (1981)) were also suitable for spatial structure analysis. Detailed structural analysis of the saturated hydraulic conductivity (Ks) of these two fields demonstrated how this type of analysis may reveal field characteristics that are not apparent from conventional statistical analysis. Using the Akaike Information Criterion for model discrimination, the three-dimensional spatial distributions of lnKs of both fields were shown to be described best by a spherical covariance function and a linear drift function. The Hamra field of (Russo and Bresler (1981)) had a much larger deterministic drift component and a smaller stochastic component than the Panoche field of (Nielsen, Biggar, and Erh (1973)). The stochastic component of lnKs in the Bet-Dagan field possessed a large nugget variance (40 percent of total) and was characterized by an integral scale of J = 14.5 m, as compared with J = 8.1 m and a small nugget variance (13 percent of total) in the Panoche field.

Transient one-dimensional vertical transport of mixed Na/Ca solutions in a saturated-unsaturated homogeneous soil during infiltration was studied using a modified version of the transport model of Bresler (1973a), taking into account physicochemical interactions between the soil solution and the soil matrix expressed in terms of changes in the hydraulic conductivity and retentivity functions, anion exclusion, and cation exchange. These latter, in turn, were derived from theoretical considerations based on the mixed ion diffuse double-layer theory, the structure of the clay particles, the pore size distribution of the soil, and hydrodynamic principles. The solute transport was analyzed for three soils of different textures and for different sets of boundary and initial conditions. Results of the analyses suggested that the effect of soil-matrix-soil-solution interactions on the transport of water and solutes may be significant and generally increases as the wetting zone soil water content and soil solution sodium adsorption ratio increase, as its solute concentration decreases, as the clay fraction of the soil increases; and as the soil texture becomes finer. For a given soil the magnitude of these interactions and their effect on the transport process are affected by the surface boundary conditions and the initial conditions. For a given set of boundary and initial conditions the retardation of the water and solute movement due to the soil-solution-soil-matrix interactions, relative to a reference inert state, depends on the soil texture. The implication of this finding with respect to solute transport in spatially variable fields is discussed briefly.

In heterogeneous fields in which the soil water properties and the initial salinity vary under a deter- ministic water application rate uniform throughout the field, the salinity during leaching might also differ from place to place. A geostatistical approach was used to investigate the spatial variability of three soil properties: the saturated hydraulic conductivity Ks, the soil characteristic α = d log (K)/d(ψ) (K being the soil hydraulic conductivity and ψ being the pressure potential), and the dispersivity λ, as well as the initial salinity EC0, by using actual measured field data. These properties were used as input parameters of a simplified water and salt flow model, which in turn was coupled with the conditional simulation method to analyze the salinity profile and its spatial distribution during leaching in a 187-ha plot of land. Analysis of the results showed that 107 hours of continuous leaching (6527 m3 ha−1) were required to obtain an average salinity of EC* = 5 dS/m for the field layer between the soil surface and the 40-cm depth. By considering the leaching of the different sites in the field, rather than that of the entire field, it was shown that, theoretically, the amount of water for leaching required to obtain EC* = 5 dS/m uniformly throughout the field can be reduced to 4038 m3 ha−1 (a reduction of 38%). Practically, since the field had to be subdivided into small subplots and because of engineering design requirements for the water supply system, the amount of water for leaching required to obtain the above-mentioned value of EC* for the entire field (5193 m3/ha) could be reduced by only 20%.

The effect of physicochemical interactions between the soil solution and the soil matrix on the spatial variability of the soil-dependent variables during transient transport of mixed Na/Ca−Cl salts through the unsaturated zone of a large field-scale soil is analyzed. The spatial variability of the soil water and solute transport properties evaluated in part I were used as inputs to a conceptual stochastic model describing one-dimensional vertical transport of water and soil-interacting solutes through the unsaturated zone of a spatially variable soil, viewed as a collection of vertically homogeneous and independent soil columns. Using data of the spatial distributions of pertinent soil properties from the Bet Dagan field, the transport process was simulated for boundary and initial conditions pertinent to the application of low-salinity and low-alkalinity waters (rain) to the surface of a saline-alkaline soil, considering effects of the soil solution concentration and composition on the soil water and solute transport properties. The spatial variability (in terms of field averages and coefficients of variation CV) of the profiles of the soil water content θ, the chloride concentration C, and the sodium adsorption ratio SAR of the soil solution and the soil hydraulic conductivity K at given elapsed times t were presented and compared with those evaluated from simulations of the same boundary value problem where the effect of the soil solution concentration and composition on the water and solute transport properties was disregarded. Results of the analyses suggested that because of soil solution-soil matrix interactions the field-average movement of both the water and the solutes may be retarded and their spatial variability may be increased relative to the case where these interactions had not been considered. In the Bet Dagan field, after t = 5 hours of continuous infiltration, the field averages of the positions of the wetting front and the chloride front were retarded by 10 and 15%, respectively; the relative variabilities of the positions of these fronts were increased by 16 and 38%, respectively; and the equivalent effective dispersivity was increased by 18%, relative to the reference inert case.

A geostatistical approach was used to investigate the spatial variability of two measured soil hydraulic parameters, the saturated hydraulic conductivity Ks and α = d[log K(h)]/[d(h)], where K is the hydraulic conductivity and h is the soil water pressure head. Using a linearized solution to steady infiltration from a circular shallow pond, the spatial distribution of the midway soil water pressure head between emitters, hc (for a given trickle discharge rate Q and spacing between emitters, d), was calculated using the spatially variable soil parameters Ks and α. The spatial distribution of the crop yield Y was calculated from the spatially variable hc by using a quadratic expression for the Y(hc) relationships. It was found that the spatial variability of each of the measured input parameters (Ks and α) or the calculated output variables (hc and Y) is expressed through the variogram γ and can be characterized in terms of the coefficient of variation CV and the integral scale J. In the case of log Ks, for which the sill is reached at relatively small range, two different γ were tested. Results showed that the kriged values of log Ks as well as the resultant values of hc = hc(Ks α; d, Q) were relatively insensitive to the selected γ. The spatial variability of hc reduced the average (over the field) crop yield that would have been obtained when hc was uniform throughout the field. The spatially variable soil parameters Ks and α were used to calculate the spatial distribution of the spacing between emitters (d) required to achieve (for a given trickle discharge) a prespecified value of hc, uniform throughout the field. It was found that the use of a spatially variable d increased the field crop yield by 3% relative to the yield that was obtained when a single-valued d uniform throughout the field was used. The spatial distribution of d was combined with hydraulic principles to derive the spatial distribution of lateral diameter and length for engineering design requirements.

Leaching a gypsiferous-sodic soil under different soil conditions and water qualities was investigated using the transport model of E. Bresler (1973) combined with the salinity model of C. W. Robbins et al. (1980b). Incorporating field-measured values of the soil hydraulic properties and the dispersivity coefficient, profiles of the soil solution concentration of Cl, SO4, Ca and Na, and soil solution electrical conductivity (EC) obtained under intermittent leaching were simulated for a set of representative cases typical to gypsiferous-sodic soils. For a given soil depth and amount of applied water the simulated soil solution concentration of chloride was affected only by the initial soil salinity and the salinity of the applied water. Simulated soil solution concentrations of calcium, sodium, and sulfate were also affected by the initial soil exchangeable sodium percentage (ESP), initial gypsum content, and the cation exchange capacity (CEC) of the soil. Sensitivity analysis of the model results suggest that soil CEC, initial soil salinity and ESP, and gypsum content have a significant impact on the soil salinity profiles during leaching. Possible applications of the combined model results for salinity control were demonstrated. Expressing the soil system response in terms of the soil solution EC, rather than the specific ion concentrations, the combinations of water quantity and quality (salinity) which yield a given soil solution EC at a given soil depth for different initial conditions and soil CEC were calculated. For a given salinity of the applied water, the amount of water required to maintain a given EC at a given soil depth increased as the initial soil salinity, soil ESP, gypsum content, and CEC rose. The calculated isosoil salinity lines can be used for economic evaluation of the least cost combination of the applied water quantity and quality which will ensure that a critical salinity at a given soil depth is not exceeded.

Transient one-dimensional vertical transport of interacting solutes through a saturated-unsaturated stratified soil profile during infiltration was studied using a modified version of the transport model described in part 1 [Russo, this issue], taking into account the depth variations of the coupling interacting soil input parameters. The solute transport was analyzed for two sequences of soil layering: a fine-textured soil overlying coarser textured soils and a coarse-textured soil overlying finer textured soils. Each sequence was analyzed for two different initial conditions and a given surface boundary condition pertinent to a common irrigation practice using nonsaline water. Results of the analyses suggested that the effect of soil-matrix-soil-solution interactions on the transport of water and solutes may be significant and generally depends on the sequence of the soil layering and on the shape of the initial soil solutions's sodium adsorption ratio (SAR) profile. For a given initial soil solution's SAR profile the effect of soil-solution-soil-matrix interactions on the transport process is greater when a fine-textured soil is overlying coarser textured soils. For a given sequence of layering the movement of the interacting solutes may be retarded or accelerated relative to a reference inert state, depending on the shape of the initial soil solution's SAR profile relative to the sequence of the layering. The possible effects of soil-solution-soil-matrix interacting on the spatial variability of solutes on a field scale are discussed briefly.

In a heterogeneous field in which the soil water properties vary under a “deterministic” uniform trickle irrigation system, the midway soil-water pressure head hc and the yield of a crop also differ from place to place. These differences may, in turn, reduce the average (over the field) yield relative to the yield that would be obtained if the soil was uniform throughout the field. A field experiment was conducted to test the hypothesis that this yield reduction may be eliminated by using a spatially variable trickle irrigation system. Twenty-five plots (200 m2 each) were established on a 30-m2 grid. Half of each plot was equipped with a standard trickle irrigation system with constant spacing between emitters of d = 50 cm (control plots), and the other half was equipped with a trickle irrigation system for which the spacing between the emitters was selected by using the pertinent hydraulic properties (the saturated hydraulic conductivity Ks and the soil parameter α) according to the procedure of Bresler (1978) as described in paper 1 (Russo, 1983b). Values of hc measured at different times, as well as the total fruit yield Y of bell pepper (Capsicum frutescens var. “Maor”), were used to estimate the seasonal and the spatial distributions of hc and the spatial distribution of Y and their moments. The variograms of hc and Y were calculated and used to estimate their integral scales. It was found that the use of a spatially variable d relative to the use of a uniform d did not change the seasonal behavior of hc but reduced the spatial variability in hc and Y by 35% and 11%, respectively, and increased the integral scale of hc and Y by 30% and 10%, respectively, but increased the average total fruit yield by only 1.9%. The use of a spatially variable d reduced the dependence of Y on hc. This indicates that when the emitters are properly spaced, it is not the water but other factors that most influence yield. When a constant d was used, the dependence of Y of hc decreased with time. This and the relatively good agreement between the values of hc measured at the initial stages of the growing season and those calculated in paper 1 demonstrate that the concept of hc is important in the early stages of the plant's growth, when the root system is not fully developed. Both the theoretical (paper 1) and the experimental results showed that although Ks and α, as well as hc, varied considerably in the field the spatial variability of the crop yield was relatively small. This explains why the use of a spatially variable d essentially was not an improvement over the fixed d. It is suggested that this study will be considered as a methodological one, which can be adapted to solve practical problems associated with field spatial variability.

Parameter estimation procedures involving solution of the inverse flow problem pertinent to a transient experiment constitute a powerful method of determining the soil hydraulic properties. One of the problems associated with these procedures is the selection of an appropriate model to describe the soil hydraulic properties. Three different models were considered: those of van Genuchten (VG) and Brooks and Corey (BC) model, and the exponential model of Gardner for the hydraulic-conductivity-soil-water-pressure relationships, coupled with a new expression for the soil-water-content-soil-water-pressure relationship (GR model). For a given model the soil hydraulic properties were determined from simulated and measured outflow experiments supplemented with water content at a soil water pressure head of −15,000 cm H2O, using the parameter estimation procedure of Kool et al. (1985a) and data from two soils: a hypothetical sandy loam soil (assumed to be described by the VG model) and a silt loam soil. Model validation tests were performed, and the most appropriate model was selected from the candidate models by discrimination tests using the Akaike Information Criterion (AIC). In the case of the hypothetical sandy loam, the VG model was found to be the most accurate and most consistent with the data. The performance of both the BC and the GR models (in terms of the AIC) was less good but very similar to each other. In the case of the silt loam soil, the VG model with additional parameter m, which accounts for the correlation between pores and for the flow path tortuosity (m=2.02) rather than a constant m=0.5, was found to be most accurate and most consistent with the data. The performance of the BC and the GR models was similar when m was considered as an unknown parameter instead of a constant.