SSSAJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


Published online 27 October 2005
Published in Soil Sci Soc Am J 69:1943-1954 (2005)
DOI: 10.2136/sssaj2005.0051
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (6)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Crescimanno, G.
Right arrow Articles by Garofalo, P.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Crescimanno, G.
Right arrow Articles by Garofalo, P.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Crescimanno, G.
Right arrow Articles by Garofalo, P.
Related Collections
Right arrow Vadose Zone Processes and Chemical Transport
Right arrow Inverse Procedures/Parameter Estimation
Right arrow Water Flow Models

Soil Physics

Application and Evaluation of the SWAP Model for Simulating Water and Solute Transport in a Cracking Clay Soil

Giuseppina Crescimanno* and Paolo Garofalo

Università degli Studi di Palermo, Dipartimento ITAF–Sezione Idraulica, Viale delle Scienze, 13 90128 Palermo, Italy

* Corresponding author (gcrescim{at}unipa.it)


    ABSTRACT
 TOP
 ABSTRACT
 INTRODUCTION
 THEORY OF SWAP
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
In Sicily, the increasing scarcity of quality water is leading to irrigation with saline water in soils having a considerable susceptibility to cracking. Irrigation systems involving high application rates are used in these irrigated areas, and bypass flow during irrigation is thus prevalent. Adoption of management practices accounting for cracking is therefore necessary to prevent salinization and land degradation. In this paper, water flow and solute transport in a Sicilian cracking soil irrigated with saline water was simulated by using the soil-water-atmosphere-plant environment (SWAP) model, and the simulated results compared with measured values of soil moisture and salinity. The soil hydraulic parameters were obtained by inverse method based on multi-step outflow experiments, adopting two different sets of hydraulic parameters/functions, that is, (i) the van Genuchten-Mualem, (VGM model) and (ii) the Brutsaert retention equation coupled with the hydraulic conductivity model proposed by Gardner (B-G model). The results obtained using field measurements from four soil profiles of a cracking clay soil showed that SWAP provided accurate predictions of water content, {theta}, when the soil hydraulic properties were expressed according to the B-G model. Using the B-G hydraulic parameters/functions, the model was calibrated with reference to the dispersivity (Ldis). A calibration value of about 20 cm was found for the four different profiles. In the conditions occurring in the Sicilian area where we are focusing our attention, the predictive errors associated with the simulated ECsat, can be considered acceptable if the purpose of application is to predict the influence of salinity on crop yield.

Abbreviations: ADE, advection-dispersion equation • B-G, Brutsaert–Gardner • COLE, coefficient of linear extensibility • ECsat, electrical conductivity of saturated soil extract • ECw, electrical conductivity of irrigation water • ESP, exchangeable sodium percentage • MSTEP, multi-step • RMSR, root mean squared residual • SAR, sodium adsorption ratio • SCIM, suction crust infiltrometer method • SWAP, Soil-Water-Atmosphere-Plant environment • VGM, van Genuchten-Mualem


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 THEORY OF SWAP
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
THE AVAILABLE AMOUNT of fresh water for agriculture, and specifically for irrigation, is decreasing all over the world. The quality of irrigation water is also deteriorating, and saline/sodic waters are increasingly used in many arid and semiarid regions of the world (FAO, 1992). From a global perspective, irrigated agriculture makes an essential contribution to the food needs of the world. However, vast areas of irrigated land are increasingly threatened by salinization and land degradation (Szabolcs, 1992).

The effects of salinity are manifested in loss of land, reduced rates of plant growth, reduced yields and, in severe cases, total crop failure (Rhoades and Loveday, 1990). The salt composition of the soil water also influences the composition of cations of the exchangeable complex of soil particles, affecting soil structural and hydraulic properties (Crescimanno et al., 1995). All these factors should be considered in developing sustainable irrigation strategies finalized to protect soil from salinization (Crescimanno et al., 2004).

In Sicily, the increasing scarcity of good quality water coupled with intensive use of soil under semiarid to arid climatic conditions, is leading to irrigation with saline water on soils having a high shrink-swell potential and susceptibility to cracking (Crescimanno and Provenzano, 1999). These soils are irrigated in the summer season, when the cracks open up, by sprinkler systems, which involve high application rates. Because of these high application rates, bypass flow, that is, the rapid transport of water and solutes via macropores or cracks to subsoil or to groundwater (Bouma, 1991; Crescimanno, 2001), is prevalent during irrigation.

Laboratory investigations performed on undisturbed soil columns sampled from these areas showed that salinization or leaching occurred during bypass flow depending on the concentration of the applied solution compared with the concentration of the pore solution (Crescimanno and De Santis, 2004); the efficiency of salt-leaching was found to depend on crack volume (Crescimanno et al., 2002). The low values of the sodium adsorption ratio (SAR) of irrigation water, and the low values of the exchangeable sodium percentage (ESP) measured in these soils, indicated no risk of sodication under current conditions. These results suggest that management strategies accounting for cracking and bypass flow should be adopted to prevent salinization and land degradation in these or similar areas (Crescimanno et al., 2004). Application of physically based models simulating water and solute transport, and predicting soil salinity, expressed by measurement of soil electrical conductivity (U.S. Salinity Laboratory Staff, 1954) or by concentration of the pore solution, represents an essential tool for developing management scenarios suitable to prevent salinization, as these models can be used to analyze and compare different options.

Enormous advances have been made during the last decades in modeling flow and transport processes in the vadose zone between the soil surface and the groundwater table. Simunek et al. (2003) made a complete review of the different approaches developed for modeling preferential and nonequilibrium flow and transport in the vadose zone, indicating the need for calibration/validation of the different models by comparison with field data. However, although numerical models have become more and more sophisticated, their success and the reliability of predicted values are critically dependent on accurate information of soil hydraulic parameters (Wagner et al., 1998). The water retention-hydraulic conductivity model proposed by VGM model (Mualem, 1976; van Genuchten, 1980) has been widely used to characterize soil hydraulic properties. Many databases and pedotransfer functions have been obtained using the VGM parameters (Nemes et al., 2003). However, the mathematical constraints in the VGM model may determine a poor performance of this model to represent the unsaturated hydraulic conductivity of heavy clayey, structured soils (Fuentes et al., 1992; Schaap and Leij, 2000).

Vereecken et al. (1989) investigated the applicability of different retention equations to a large number of Belgian soils having different texture. They found that the equation proposed by Brutsaert (1966), (B), which is the same as the van Genuchten (VG) equation when condition m = 1 applies instead of m = 1 – 1/n, provided a better estimation of the water retention curve compared with the VG equation. Vereecken et al. (1990) also investigated the suitability of the hydraulic conductivity function proposed by Gardner (G) to fit hydraulic conductivity measurements obtained by the crust method (Booltink et al., 1991). Working on 127 cores exploring a wide texture range, they found that this model provided the best prediction of the k(h) function compared with application of different hydraulic conductivity equations.

Measurement of soil hydraulic properties is difficult, especially unsaturated hydraulic conductivity (Schaap and Leij, 2000). Estimation of soil hydraulic parameters by means of inverse modeling (Kool et al., 1987) has become an attractive alternative to traditional methods (Bitterlich et al., 2004). Popular laboratory approaches for the inverse estimation of soil hydraulic parameters have been one-step or multi-step (MSTEP) outflow methods (van Dam et al., 1994; Crescimanno and Iovino, 1995).

Whereas the general flexibility of outflow experiments for identifying hydraulic functions has been demonstrated in several studies (Hopmans and Simunek, 1999), the applicability and success of this method have been shown to depend, among other possible factors, on suitable parameterization of the hydraulic functions (Vrugt et al., 2003).

Crescimanno and Baiamonte (1999) investigated the results of using (i) the VGM model or (ii) the B retention equation coupled with the k(h) Gardner equation (Gardner, 1958) (B-G model) in a parameter estimation procedure based on MSTEP outflow experiments. The investigation was performed on some Sicilian structured clay soils, in which the suction crust infiltrometer method (Booltink et al., 1991) was used to obtain independent k(h) measurements. They found that using the B-G model made it possible (i) to estimate hydraulic conductivity functions in close agreement with the independently measured hydraulic conductivity values and (ii) to estimate a water retention curve that was consistent with the water content-pressure head values obtained during the MSTEP experiments. On the contrary, application of the same parameter estimation procedure using the VGM model led to a very poor agreement between measured and estimated k(h) and {theta}(h) values.

Van Dam et al. (1997) developed a model for fine-textured clay soils containing shrinkage cracks. This model, named SWAP, takes into account shrinking and swelling as a function of water content. The model assumes that water and solutes can move instantaneously to specified bypass depths once the infiltration capacity of the soil matrix is exceeded by rainfall rate and a critical depth of water has formed at the soil surface (Verburg et al., 1996). Since water flow in the matrix is described by the Richards equation, while water in the cracks moves from the soil surface to some specified depths, this model can be viewed as a subgroup of the dual permeability models (Jarvis, 1994). However, the advantage of SWAP compared with the dual permeability models is that only one set of hydraulic parameters/functions is needed to characterize the soil instead of the larger number of parameters necessary to characterize the different pore regions. In addition, shrinkage and cracking are specifically accounted for and flow through the cracks can be calculated.

The SWAP model provides as output the water content, {theta}, (and pressure head, h), as well as the electrical conductivity of the saturated extract, ECsat, (Rhoades, 1996). Reduction in crop yield is calculated as a function of ECsat (Maas and Hoffman, 1977).

With reference to solute transport, SWAP predicts solute concentration by using the advection-dispersion (ADE) equation, assuming a constant dispersivity (Ldis) and neglecting the process of cationic exchange. Due to these as well as to other simplifying assumptions, the accuracy of the predicted ECsat in clay soils needs to be checked by comparison with measured ECsat values.

Only a limited number of applications of the SWAP model relating to application of saline/sodic waters in agricultural soils can be found in the current literature. Smets et al. (1997) calibrated and validated SWAP to use the model to evaluate the effect of various irrigation practices on salinization and crop transpiration in Pakistan. Using measurements of water content ({theta}), pressure head (h), and ECsat, in four profiles, they found that the model accurately predicted both h and {theta}, but underestimated the predicted ECsat. Previously determined VGM hydraulic parameters were used as initial values and slightly modified in a trial-and-error process to obtain optimal calibration results. Kelleners et al. (1999) investigated the influence of spatially variable hydraulic properties on solute transport simulated by SWAP. Hydraulic properties derived from outflow experiments and expressed according to the VGM model were used as input in SWAP. However, no comparison with measured values was reported in their paper. Tedeschi and Menenti (2002) applied SWAP to develop management scenarios in Southern Italy. Calibration was performed with reference to {theta} values, using hydraulic parameters measured under field conditions and expressed according to the VGM model. However, no comparison between measured and predicted ECsat was reported in their paper. Singh (2004) used SWAP to develop guidelines for irrigation planning in India. However, no information about the hydraulic parameters adopted in their simulations was given in their paper.

The objective of this paper was to evaluate the ability of SWAP to accurately predict both {theta} and ECsat in four profiles of a Sicilian cracking soil located in a field (Mazara del Vallo, TP) where saline water is used for irrigation, with a risk of salinization identified in the course of previous investigations (Crescimanno, 2001; Crescimanno and De Santis, 2004).

The soil hydraulic parameters were determined by the inverse method based on multi-step outflow experiments by representing the soil hydraulic functions by the VGM model and by the B-G model (Crescimanno and Baiamonte, 1999). Using measurements of water content obtained in the field over a 2.5-yr period, the ability of SWAP to accurately predict the water content was checked using both the VGM model and the B-G model. After selection of the soil hydraulic models/parameters providing the best prediction of {theta}, the ability of SWAP to predict ECsat was tested against measurements of ECsat obtained on soil samples collected in the field at the same sampling dates as for the water content. The model was calibrated with reference to the Ldis parameter, representing the dispersivity in the ADE (Warrick, 2003). This equation remains the foundation on which most analyses of solute transport in porous media have been based.


    THEORY OF SWAP
 TOP
 ABSTRACT
 INTRODUCTION
 THEORY OF SWAP
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
One-dimensional, vertical, transient, unsaturated flow in the SWAP model (van Dam et al., 1997) is described by the Richards equation, which is solved numerically.

The shrinkage characteristic is expressed by the model proposed by Kim (1992):

[1]
where e = Vp/Vs is the void ratio, (cm3 cm–3), u = Vw/Vs is the moisture ratio, (cm3 cm–3) and {alpha}sh, ßsh, {gamma}sh are dimensionless fitting parameters; Vp is the total pore volume, Vs is the solid volume and Vw is the water volume. The shrinkage characteristic allows the calculation, at a certain soil depth or node i, of the relative cross-sectional area of the cracks at the soil surface, Ac (cm2 cm–2), according to the following steps (Bronswijk, 1989):

Vs = 1 – {theta}s, is the solid volume, where {theta}s, is the volumetric water content at saturation (cm3 cm–3);

u = {theta}i/Vs is the moisture ratio, with {theta}i is the water content of node i, following from the solution of the Richards equation at this time step;

e is calculated from Eq. [1];

Vp = e x Vs is the total pore volume;

{Delta}V = {theta}sVp is the shrinkage soil volume with respect to maximum soil volume;

— vertical subsidence {Delta}z follows from 1 – {Delta}V/V = (1 {Delta}z/z)rs, where rs is the geometry factor and V is the volume (cm3) of soil cube with side z (cm);

Vc = {Delta}V{Delta}z is the volume vertical crack;

Ac = Vc/(1 – {Delta}z). The model calculates a crack volume if {theta}i is lower than the water content corresponding at the beginning of the structural shrinkage (Crescimanno and Provenzano, 1999).

The matrix and crack infiltration at a given rainfall intensity, P (cm h–1), are calculated as follows:

[2a]

[2b]

[3a]

[3b]
where Imax is the maximum infiltration rate of the soil matrix, (cm h–1), Ic is the infiltration rate into the cracks, (cm h–1), and Am (cm2 cm–2) is the relative area of soil matrix.

Solute transport in the model is described with the ADE, (Warrick, 2003):

[4]
where {theta} is the volumetric water content, c is the solute concentration (g cm–3) in soil water, t (h) is time, q is the soil water flux (positive upward) (cm h–1), z is the vertical coordinate with the origin at the soil surface (positive upward) (cm), and D is the apparent diffusion coefficient (cm2 h–1).

The apparent diffusion coefficient D accounts for both mechanical dispersion and effective ionic or molecular diffusion according to:

[5]
where Ddif is the effective ionic or molecular diffusion coefficient (cm2 s–1) and Ddis is the dispersion coefficient (cm2 s–1). Under laminar flow conditions, Ddis is proportional to the pore water velocity, v (cm s–1), according to the following equation:

[6]
where Ldis is the dispersivity (cm), which reflects the complexities of the flow pathways and the heterogeneity in local fluid velocities in the direction of flow (Beven et al., 1993).

Although some laboratory (Vanderborght et al., 2000) and field investigations (Forrer et al., 1999) have shown that Ldis depends on flow rate, the simplifying assumptions that (i) Ldis is invariant with time and depth and that (ii) D is linearly related to Ldis in Eq. [6] are made in SWAP. Another simplifying assumption concerning solute transport in SWAP is that no distinction is made between different cations and anions, and only the total amount of salts is considered.

The root water uptake is semi-empirically described by a sink term, which is a function of the maximum root water uptake, the soil water pressure head, and the salt concentration. The maximum root water uptake may be uniform or linearly declining with depth (Feddes et al., 1998). At the bottom of the system, boundary conditions can be described with various options, for example, water table depth, flux to groundwater, or free drainage.

Saline soil moisture conditions reduce root water uptake owing to an increased osmotic head of the soil water. In SWAP, the osmotic head is added to the matrix head, and the total head is used to derive the reduction factor for root water uptake. SWAP reduces crop transpiration owing to water and salinity stress only, while all other conditions are assumed to be optimal.

According to Maas and Hoffman (1977), ECsat is related to relative crop yield, Yr, by the following relationship:

[7]
where i is the threshold salinity value, and p is the percentage decrement value for unit increase of salinity in excess of the threshold. For grapes, i is equal to 1.5 dS m–1 and p = 9.6% (FAO, 1992).


    MATERIALS AND METHODS
 TOP
 ABSTRACT
 INTRODUCTION
 THEORY OF SWAP
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Data Collection
Data collection was performed in a 25 by 25 m plot located in Sicily (37°40'55'' N latitude; 12°38'50'' E Longitude) where irrigation with saline waters is performed on grapes by a sprinkler system, which allows high application rates at the soil surface.

The soil physical and chemical characteristics, together with the pedological classification, are reported in Table 1. The electrical conductivity, ECw, and SAR of irrigation water were 2.1 dS m–1 and 4.0 (mmol L–1)0.5, respectively. Gravimetric water content, U, was determined on undisturbed soil cores sampled at 30 and 45 cm in the selected profiles at different dates (from 14 July 1998–31 Dec. 2000). Bulk density ({rho}b) was determined from the measured shrinkage curves; U and {rho}b(U) were used to calculate volumetric water content, {theta}, which therefore accounted for a variable soil volume. Soil saturated extracts were prepared using the soil collected in the field. Soil electrical conductivity, ECsat, was determined on these extracts by a conductivimeter (Crison, Micro CM 2002, Crison Instruments, Barcelona, Spain).


View this table:
[in this window]
[in a new window]
 
Table 1. Classification, physical, and chemical properties, COLE and shrink-swell potential of the considered soils.

 
Soil Shrinkage and Hydraulic Characteristics
Replicated soil cores having different sizes according to the physical and hydraulic characteristics to be measured were sampled from the different horizons in the four profiles called Baglio1, Baglio2, Baglio3, and Baglio4. Undisturbed soil cores (diameter d = 8.5 cm, height H = 11.5 cm) were sampled to measure the soil shrinkage curve, which was determined by measuring vertical and horizontal shrinkage (Crescimanno and Provenzano, 1999). The model proposed by Kim (1992) was fitted to the measured u, e values. The coefficient of linear extensibility, COLE (Grossman et al., 1968), indicating the shrink-swell potential (Parker et al., 1977), was also calculated.

Soil columns (d = 20 cm, H = 20 cm) were sampled to measure the saturated hydraulic conductivity of the soil matrix, ksat, and some k(h) values close to saturation, by the suction crust infiltrometer method, SCIM (Booltink et al., 1991). Replicated soil cores (d = 8.5 cm, H = 5 cm) were sampled to determine soil hydraulic parameters/functions by inverse method based on MSTEP outflow experiments. The saturated water content, {theta}s, was assumed to be equal to the water content value measured by a hanging water column apparatus (Burke et al., 1986) at a pressure head value of –2 cm. This {theta}s value was found to be consistent with the calculated porosity [{theta}s = us/(1+es)]. The MSTEP experiments were performed in pressure cells by applying three successive steps with pneumatic pressures ranging from 10 to 40 cm, from 40 to 70 cm, and from 70 to 800 cm (Crescimanno and Iovino, 1995). After the MSTEP experiments, the cores were put in a pressure plate apparatus to measure the water content at h = –15300 cm, that is, wilting point, {theta}wp. Independent measurement of {theta}s and {theta}wp was necessary because only the pressure range from 10 to about 1000 cm can be explored in the pressure cells used for the MSTEP.

Parameter estimation was performed according to Crescimanno and Baiamonte (1999), representing the soil hydraulic functions by the VGM model (Mualem, 1976; van Genuchten, 1980):

[8]

[9]
where h (cm) is the pressure head, {theta}s is the volumetric water content at saturation, {theta}r is the residual water content, k is the unsaturated hydraulic conductivity (cm h–1), ksat is the saturated hydraulic conductivity (cm h–1), {alpha} (cm–1), and n, m, and {gamma} are empirical parameters, with m = 1 – 1/n; and, alternatively, by the equation proposed by Brutsaert (1966), for the water retention curve:

[10]
coupled with the model proposed by Gardner for the hydraulic conductivity function k(h):

[11]
ß and {lambda} are empirical parameters.

The hydraulic model represented by Eq. [10] and [11] (B-G model) couples a {theta}(h) function with a closed-form equation not derived from the {theta}(h) function. Instead, the VGM model used the statistical pore-size distribution model of Mualem (1976) to obtain a predictive equation for the unsaturated hydraulic conductivity function in terms of water retention parameters.

Optimization was performed on the outflow volumes, supplemented by four {theta}(h) values obtained during the MSTEP experiments ({theta} values at –10, –40, –70, and –800 cm) and by the k(h) values obtained by the SCIM method. Parameter estimation was performed by fixing both the saturated water content and the saturated hydraulic conductivity, ksat, at the measured values. Optimized parameters were therefore {theta}r, {alpha}, {gamma}, and n in the VGM model, and {theta}r, {alpha}', ß, {lambda}, and n' in the B-G model.

Model Application and Calibration
Climatic data (rain intensity, maximum and minimum temperature, rainfall height) recorded daily from 8 July 1998 to 31 Dec. 2000 by a rain gauge located in the field were used as input in SWAP. Annual rainfall in 1998, 1999, and 2000 was 390 mm in average and the mean annual reference evapotranspiration in 1998, 1999, and 2000 was 1450 mm.

Irrigation was simulated according to the scheduling adopted in the years considered, during which, because of water scarcity, very limited irrigation amounts were supplied (66 mm in 1998, 48 mm in 1999, and 24 mm in 2000). The annual irrigation amount supplied under normal conditions is about 120 mm.

Aerial photography was used to determine the soil cover fraction in the field considered; a root distribution characterized by 60% roots in the 30- to 70-cm layer, and by 20% both in the 0- to 30-cm and in the 70- to 100-cm layers, was assumed, on the basis of the results of previous investigations (Barbagallo et al., 2004).

Simulations were performed by using a bottom boundary condition of freely draining profile. The VGM and the B-G hydraulic parameters were used to simulate water transport, and the accuracy of predicted {theta} values was evaluated by calculating the Root Mean Squared Residual, RMSR{theta}, between measured and predicted {theta}:

[12]
where N is the number of measurements.

To check systematic errors between measured and predicted {theta}, the predicted {theta} values were regressed against the measured values, and the hypotheses that the slope (b) of the regression line was not significantly different from 1, and that the intercept (a) of the regression line was not significantly different from 0, were checked. The Durbin-Watson (DW) test was used to check whether the random errors in the regression line exhibited autocorrelation.

After selection of the hydraulic parameters/functions providing the best agreement between measured and predicted {theta}, the model was calibrated with reference to the Ldis parameter. On the basis of field investigations (Warrick, 2003), values ranging from 5 to 25 cm with a 1-cm step were used for model calibration. For each considered Ldis, the RMSR associated to the predicted ECsat, RMSRECsat, was calculated by Eq. [12]. Owing to the assumption of a unique Ldis for the whole profile, RMSRECsat was calculated by considering the whole set of ECsat measured at both 30 and 45 cm.

To check systematic errors between measured and predicted ECsat, the predicted ECsat was regressed against the measured values, and the hypotheses that the slope (b) of the regression line was not significantly different from 1 and that the intercept (a) of the regression line was not significantly different from 0, were checked. The DW test was again used to check whether the random errors in the regression line exhibited autocorrelation.


    RESULTS AND DISCUSSION
 TOP
 ABSTRACT
 INTRODUCTION
 THEORY OF SWAP
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Shrinkage and Hydraulic Characteristics
A good fit of the Kim model to the experimentally measured values of void ratio, e, and moisture ratio, u, Fig. 1 , was obtained for the Ap horizon of the Baglio1 profile. The same good fit, not shown for brevity's sake, was found for the other soil horizons and profiles. This indicated the suitability of the Kim model to accurately represent the soil shrinkage curve.



View larger version (20K):
[in this window]
[in a new window]
 
Fig. 1. Shrinkage characteristic obtained for Ap horizon of Baglio1 profile. The continuous line represents the Kim model fitted to the measured (u, e) values.

 
The COLE values calculated for the different horizons made it possible to classify the soils as having a shrink-swell potential (Parker et al., 1977) from medium to very high (Table 1). This indicated that the soils had a considerable susceptibility to cracking at decreasing water content.

The parameter estimation procedure based on the B-G model provided an estimated k(h) function in close agreement with the k(h) values measured by the SCIM (Fig. 2b , Ap horizon of the Baglio1 profile). Good agreement was also observed between the predicted {theta}(h) function and the measured ({theta}, h) values obtained by the MSTEP experiments (Fig. 2a). However, an unsatisfactory estimation of the k(h) function was obtained for this soil, Fig. 2b, using the VGM model. The same result, not shown in Figures for brevity's sake, was found for the other profiles and horizons.



View larger version (17K):
[in this window]
[in a new window]
 
Fig. 2. Water retention curves (a) obtained by the parameter estimation method using the van Genuchten–Mualem (VGM) ({blacktriangleup}{blacktriangleup}{blacktriangleup}) and Brutsaert-Gardner (B-G) (__) models including the SCIM measurements in the optimization procedure and fixing the ksat; and using the VGM model (---) but not including the suction crust infiltometer method (SCIM) measurements in the optimization procedure and estimating the ksat. {theta}wp is the water content independently measured at h = –15300 cm. Hydraulic conductivity functions (b obtained from the parameter estimation method using the VGM ({blacktriangleup}{blacktriangleup}{blacktriangleup}) and B-G (__) model including the SCIM measurements in the optimization procedure and fixing the ksat; and using the VGM model but (---) not including the SCIM measurements in the optimization procedure and estimating ksat.

 
Analysis of the VGM and the B-G hydraulic parameters obtained for the four profiles, Table 2, indicated that higher {theta}r values than expected for clay textures were obtained for almost all the soils using both the VGM and the B-G model. However, the estimated {theta}r, should be considered a fitting parameter rather than a parameter with its physical meaning (Nimmo, 1991). The water retention curve predicted by the B-G model accurately matched the water content independently measured at –15300 cm, that is, wilting point, {theta}wp (Fig. 2a). This demonstrated the good prediction of the water retention function at the lowest {theta} values. A good matching between the measured {theta}wp and that predicted using the B-G parameters was observed for all the other profiles and horizons.


View this table:
[in this window]
[in a new window]
 
Table 2. Hydraulic parameters determined according to the van Genuchten-Mualem (VGM) model and to the hydraulic conductivity equation proposed by Gardner coupled with the Brutsaert retention equation (B-G model).

 
In contrast, the predicted VGM water retention curve did not match {theta}wp (Fig. 2a). Inclusion of {theta}wp in the optimization procedure did not improve prediction of the water retention curve. Since almost the same {theta}r values were obtained using either the B-G or the VGM model (in Baglio2 alone, the {theta}r optimized using B-G was slightly lower than that estimated using VGM), the poor fit of the VGM model to the water retention curve did not depend on the estimated {theta}r, leading to {theta} values, Table 2, which were higher than those expected for these textures and than those obtained using the B-G model.

The poor prediction of the VGM water retention curve was due to the mathematical constraints in the VGM closed form model (Bitterlich et al., 2004), which prevented the VGM model from simultaneously fitting the {theta}(h) and k(h) measurements. This is demonstrated by the fact that this prediction can be improved (Fig. 2a) if the k(h) values are not included in the parameter estimation procedure, and if the saturated ksat is an optimized parameter (Fig. 2b). In this case, a reasonable {alpha} value is found ({alpha} = 0.024 cm–1), and the VG water retention is very similar to the B function (Fig. 2a). However, even in this case an unsatisfactory match can be observed between the {theta} predicted at 15300 cm and {theta}wp. In addition, the estimated k(h) function was lower than that previously estimated and lower than the measured k(h) values.

We can conclude that in our clayey, structured soils, use of the B-G model made it possible to obtain a better estimation of the {theta}(h) and k(h) functions compared with application of the VGM model. The poor flexibility of the VGM model to represent the hydraulic properties of clayey, structured soils was indicated by the results of previous investigations (Fuentes et al., 1992; Crescimanno and Baiamonte, 1999; Schaap and Leij, 2000; Bitterlich et al., 2004).

Prediction of Water Content
The predicted {theta}, using SWAP with the B-G parameters as input, was in close agreement with those measured (Fig. 3 , Baglio1 profile). The good match between the predicted {theta} and the lowest measured {theta} ({theta} of about 0.25 cm3 cm–3 in the Ap horizon, and {theta} of about 0.30 cm3 cm–3 in the A1 horizon, Fig. 3), confirmed that a reliable prediction of the water retention curve was obtained using the B-G model in the {theta} range from saturation to wilting point.



View larger version (35K):
[in this window]
[in a new window]
 
Fig. 3. Daily volumetric water content, {theta}, predicted by soil-water-atmosphere-plant environment (SWAP) at (a) 30 cm and at (b) 45 cm using the van Genuchten-Mualem (VGM) and the Brutsert-Gardner (B-G) parameters.

 
On the contrary, a poor agreement between simulated and measured {theta} was obtained when the VGM hydraulic parameters were used as input (Fig. 3, Baglio1 profile). A better prediction of {theta} obtained by running SWAP with the B-G parameters was also indicated by the lower RMSR{theta} obtained using the B-G model compared with those associated with the VGM parameters (Table 3). The a and b parameters of the equation found by regressing the predicted {theta} against those measured (Table 3) were not significantly different from 0 and 1 respectively at the 0.05 probability level when the B-G model was used. This indicated a satisfactory match between measured and predicted {theta}.


View this table:
[in this window]
[in a new window]
 
Table 3. Parameters indicating agreement between measured and predicted {theta} using the van Genuchten-Mualem (VGM) and the Brutsaert-Gardner (B-G) hydraulic models.

 
On the contrary, the condition that a and b were not significantly different from 0 and 1 was not verified when the VGM parameters were used, except for the Baglio4 profile. However, even in this case, a better match between measured and predicted {theta} was obtained using the B-G parameters.

The DW values (Table 4) showed that there was neither positive nor negative autocorrelation and that it was therefore possible to exclude internal dependence of errors. This result indicated that the accuracy of {theta} predicted by SWAP depended on using soil hydraulic properties that reflected the soil hydraulic behavior of the soils considered. This is consistent with previous results (Schaap and Leij, 2000) indicating that parameters in soil hydraulic functions characterizing water retention and hydraulic properties are the most important input variables for models based on numerical solutions of the variably saturated flow (Richards) equation. Since no calibration was performed to adjust the estimated hydraulic properties, the good match between the measured and the simulated {theta}, obtained using the B-G soil hydraulic properties, proved that the parameter estimation procedure adopted provided a reliable estimation of the {theta}(h) and k(h) functions. Our results also suggest that soil hydraulic parameters derived by pedotransfer functions from textural data and based on VGM parameters (Nemes et al., 2003), for clayey, structured soils may lead to inaccurate simulation of water transport when simulation models based on the Richards equation are used. Trial and error processes are often used to adjust the VGM hydraulic parameters until a good match between measured and predicted {theta} is obtained (Smets et al., 1997). However, if this procedure is used, validation of the adjusted hydraulic parameters should be preliminary to any further applications of models for predictive and/or management purposes.


View this table:
[in this window]
[in a new window]
 
Table 4. Parameters indicating agreement between measured and predicted saturated electrical conductivity of saturated soil extract (ECsat) using the Brutsaert-Garner (B-G) hydraulic model.

 
Water Storage in Cracks
Crack volume as a percentage of the volume at saturation, (Fig. 4) , calculated from the shrinkage curve, was maximum during the summer season, when the soil was drying and irrigation was necessary. Water storage in the cracks calculated in the whole profile during the simulation period (Fig. 5) represents the fraction of the applied water (rainfall and/or irrigation) not flowing directly into the matrix, but stored in the cracks before infiltrating into the matrix. The cumulative flow from the cracks to the matrix was rather low as a consequence of the very limited amount of irrigation supplied during the simulation period. However, this value would considerably increase with the application of a greater amount of irrigation. An increase in water storage in the cracks and in this flow could be useful to promote leaching of salts accumulated in the upper profile when alternating waters of different salinities are used in irrigation management (Crescimanno et al., 2002).



View larger version (22K):
[in this window]
[in a new window]
 
Fig. 4. Daily volumetric water content, {theta}, predicted by soil-water-atmosphere-plant environment (SWAP) at 30 cm using the B-G parameters; crack volume at 30 cm as percentage of volume at saturation, Vcr/V.

 


View larger version (15K):
[in this window]
[in a new window]
 
Fig. 5. Daily water storage in cracks and cumulative flux from cracks to matrix vs. time

 
Calibration of Ldis
The ECsat, predicted by SWAP for different Ldis values (Fig. 6a and 6b , Baglio1 profile) decreased at decreasing Ldis. This behavior is justified by the fact that at decreasing Ldis, all of the solutes tend to move at the same velocity v (cm s–1) and to arrive at the bottom profile after a time which approaches tb = L/v, where L (cm) is the length of the pathway. At decreasing Ldis, solute transport is by convection and tends to mimic "piston flow." The effect of dispersion is to spread the solute front, causing some early (t < tb) and late (t > tb) arrival of solutes, with respect to time tb. Ideally, for a pulse input of solute at the soil surface, at increasing Ldis, the cumulative quantity of solute leaving the bottom profile for t > tb decreases, and the storage of solute in the profile increases (mass conservation principle). This explains why the cumulative solute flux (mg cm–2 d–1) leaving the bottom profile, calculated by SWAP, ranged from a value of 71.08 mg cm–2 d–1 at Ldis = 5 cm, to a value of 45.43 mg cm–2 d at Ldis = 25 cm, and ECsat increases at increasing Ldis (Fig. 6). As can be seen in the figures, Ldis is a scale factor and does not influence the time evolution of the predicted ECsat.



View larger version (43K):
[in this window]
[in a new window]
 
Fig. 6. Daily electrical conductivity of saturated soil extract, ECsat, predicted at (a) 30 cm and at (b) 45 cm using the Brutsert-Gardner (B-G) parameters at different values of the dispersivity, Ldis.

 
The RMSRECsat calculated on the whole set of ECsat values measured at the depths of 30 and 45 cm, as a function of Ldis (Fig. 7) , indicated that the minimum RMSRECsat corresponded with Ldis equal to 22 cm for Baglio1, to 21 cm for Baglio2, and to 18 cm for both Baglio3 and Baglio4. The Ldis values obtained for the four profiles ranged in a very narrow interval (18–22 cm). As can be seen in Fig. 7, the RMSRECsat values within this range of Ldis are almost constant, and an average value of 20 cm could be assumed for the four profiles without affecting the accuracy of the predicted ECsat.



View larger version (32K):
[in this window]
[in a new window]
 
Fig. 7. Root Mean Squared Residual, RMSRECsat, associated with the predicted electrical conductivity of saturated soil extract, ECsat, vs. dispersivity, Ldis. RMSRECsat refers to all the measurements collected at 30 and 45 cm.

 
It is interesting to notice that although different mean ECsat values were measured in the four profiles during the simulation period (Table 4), almost the same Ldis value (Ldis = 20 cm) was found as the calibration value. If confirmed for other soils, this calibration procedure provides an "effective" dispersion coefficient, which reflected the complexities of the flow pathways and heterogeneity in local fluid velocities in the flow direction (Beven et al., 1993; Forrer et al., 1999). It would also mean that the ADE (Eq. [4]) is applicable in a functional sense in which the mean transport velocity reflects the mass flux of water averaged over some unit areas in the system (Beven et al., 1993).

Prediction of the Electrical Conductivity of Saturated Extract
The a and b parameters of the equation found by regressing the ECsat predicted at the calibrated Ldis against those measured (Table 4) were not significantly different from 0 and 1 respectively, at the 0.05 probability level both at 30 and at 45 cm. This indicated that no systematic errors were associated with the predicted ECsat. However, visual observation of the predicted ECsat (Fig. 6) showed that in some cases the predicted ECsat did not match the temporal evolution of the measurements. It was therefore possible that local nonequilibrium conditions existed during the simulation period. This would invalidate some of the assumptions of the ADE equation and cause inaccurate prediction of ECsat. However, the DW statistic (Table 4) proved that the random errors in the estimated ECsat were independent (the significance level was always 0.05), excluding internal dependence of errors.

To further evaluate the prediction of ECsat provided by SWAP, we compared the RMSRECsat with the ECsat variation determining a variation of 5% in crop yield, which is equal to 0.521 dS m–1 according to Eq. [7]. The RSMRECsat values, reported in Table 4, were always lower than this value. As a consequence, we concluded that the predictive errors associated with the simulated ECsat could be considered acceptable if the purpose of application is to predict the influence of salinity on crop yield. However, our results were obtained for a nonsodic soil, under a condition of SAR of irrigation water close to the soil ESP, and consequently sodium in the solution and in the exchange complex should be almost in equilibrium. Under this condition, the simplifying assumption that the salts are not adsorbed to soil solids could be considered acceptable. Prediction of ECsat provided by SWAP should be carefully checked when irrigation is performed on sodic soils, or when sodication can be the consequence of using irrigation waters with SAR higher than soil ESP (Crescimanno and De Santis, 2004).


    CONCLUSIONS
 TOP
 ABSTRACT
 INTRODUCTION
 THEORY OF SWAP
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
The SWAP model applied to four Sicilian cracking profiles provided a satisfactory prediction of {theta} when the soil hydraulic characteristics were represented using the B-G model. A less accurate prediction of {theta} was obtained using the VGM model to represent soil hydraulic properties. These results confirm previous results indicating that the B-G model is suitable to represent the soil hydraulic functions of clayey, structured soils, also indicating that accurate simulation of water transport with SWAP depends on the use of hydraulic parameters and functions which adequately represent soil hydraulic behavior.

The narrow range of variation of the calibrated Ldis in four soils having different hydraulic parameters, with an average value of 20 cm for the four considered profiles, seemed to indicate that the calibrated Ldis was a lumped parameter representing the irregularities in the flow pathways. The possibility of using the same Ldis value to predict ECsat for additional soil profiles located in the same area will be further checked.

For this Sicilian area where salinization is the main consequence of irrigation, the predictive errors associated with the simulated ECsat, can be considered acceptable if the purpose of application is to predict the influence of salinity on crop yield. However, the temporal evolution of the measured ECsat was not always satisfactorily predicted, probably because of locally occurring chemical processes not accounted for by the model. Further investigation will be performed to test the ability of SWAP to predict ECsat in other soil profiles and to develop management scenarios optimizing irrigation and preventing salinization.


    ACKNOWLEDGMENTS
 
Research supported by MURST (Rome, Italy), PRIN 2003: "Strategie per l'utilizzazione di acque salino/sodiche in terreni argillosi con crepacciature". The helpful comments and suggestions of two anonymous reviewers are gratefully acknowledged.

Received for publication February 17, 2005.


    REFERENCES
 TOP
 ABSTRACT
 INTRODUCTION
 THEORY OF SWAP
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
Vadose Zone JHome page
J. C. van Dam, P. Groenendijk, R. F.A. Hendriks, and J. G. Kroes
Advances of Modeling Water Flow in Variably Saturated Soils with SWAP
Vadose Zone J., May 27, 2008; 7(2): 640 - 653.
[Abstract] [Full Text] [PDF]


Home page
Soil Sci.Home page
G. Crescimanno and P. Garofalo
Management of Irrigation with Saline Water in Cracking Clay Soils
Soil Sci. Soc. Am. J., August 22, 2006; 70(5): 1774 - 1787.
[Abstract] [Full Text] [PDF]