Published in Vadose Zone Journal 3:570-591 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
REVIEWS AND ANALYSES
Rainfall-Induced Soil Surface Sealing
A Critical Review of Observations, Conceptual Models, and Solutions
S. Assouline*
The Department of Environmental Physics, Institute of Soil, Water and Environmental Sciences, A.R.O. Volcani Center, Bet Dagan 50250, Israel
* Corresponding author (vwshmuel{at}agri.gov.il).
Contribution of the Agricultural Research Organization, Institute of Soil, Water and Environmental Sciences, Bet Dagan, Israel, No. 617/03.
Received 30 July 2003.
 |
ABSTRACT
|
|---|
Rainfall-induced soil surface sealing can have severe agricultural, hydrological, and environmental effects. Seal formation is a complex phenomenon dominated by a wide variety of factors involving soil properties, rainfall characteristics, and flow conditions. It has been studied through extensive experimental investigations as well as simulation models. This study reviews some of the main issues, in terms of morphology, phenomenology, and both conceptual and empirical modeling approaches to improve our perception of the phenomenon and our ability to simulate its effects on flow processes. The effects of different factors on infiltration during seal formation and in sealed soil profiles are highlighted, including seal representation, soil and rainfall properties, and field heterogeneity. New research opportunities toward a generalized formulation of the processes involved in soil sealing and a reliable quantitative prediction of its effects on flow processes are identified. These are related mainly to the ability to quantify and predict the relationships between soil hydraulic properties and physical, chemical, and biological factors that affect the soil resistance to destruction.
Abbreviations: EC, electrical conductivity • ESP, exchangeable Na percentage • SAR, Na adsorption ratio
|
INTRODUCTION
|
|---|
SOIL SURFACE SEALING can have severe agricultural, hydrological, and environmental effects. This phenomenon decreases the infiltration rate, reduces the available water to the plant in the root zone, diminishes the natural recharge of aquifers, increases runoff and soil erosion, affects seedlings and plant growth, and decreases crop yields. Therefore, in arid and semiarid regions, soil surface sealing can induce a vicious cycle that can intensify desertification since it affects vegetation cover, enhances overgrazing, and finally leaves the soil surface bare and exposed to rainfall. However, if the seal effect is considered negative when maximum infiltration is desired, it can be beneficial when the objective is to maximize runoff, such as in water harvesting systems.
Soil surface sealing may result from different causes, such as rainfall, fire, biological activity, or mechanical or chemical treatments. In the category of rainfall-induced soil sealing, two types of seals can be identified: (i) structural seals directly related to rainfall through the impact of raindrops and sudden wetting and (ii) depositional or sedimentary seals indirectly related to rainfall that result from the settling of fine particles carried in suspension by runoff in soil depressions. This review focuses solely on structural seals.
The formation of structural seals at the surface of bare soils exposed to the direct impact of raindrops is dominated by a wide variety of factors involving soil properties, rainfall characteristics, and flow conditions. The complexity of the phenomenon has challenged many scientists, and it has been studied during the last four decades through extensive experimental investigations as well as simulation models. A large number of experimental studies have investigated the effect of the factors involved in seal formation, either during its dynamic stage, or after it has already reached its final stage when the seal layer is fully developed. Laboratory as well as field experiments were performed on various soil types under either saturated or unsaturated flow conditions, and for a wide range of natural and simulated rainfall intensities and kinetic energies. Different conceptual as well as empirical models of seal formation, properties, and effects on infiltration were suggested.
Dealing with soil surface sealing requires the ability to understand, characterize, and quantitatively represent the phenomenon and the major factors involved in it. This review does not pretend to be an exhaustive survey of all studies performed on this subject. Its aim is to review some of the main issues, in terms of morphology, phenomenology, and both conceptual and empirical modeling approaches, while highlighting recent studies, to improve our perception of the phenomenon and our ability to simulate its effects on flow processes.
The review begins with the description of the processes involved in seal formation. It then presents observations on seal morphology and properties and describes some conceptual models suggested to represent the seal layer and its hydraulic properties. The review then addresses the issue of the dynamics of seal formation, the factors identified as affecting it, and the different models suggested to quantify it. Finally, the approaches developed to model and simulate infiltration in sealed soils and during seal formation are examined, with emphasis on the effects of different factors like seal representation, soil and rainfall properties, and field heterogeneity.
|
THE PROCESSES INVOLVED IN SEAL FORMATION
|
|---|
Structural seals are formed at the soil surface by the destruction of the soil aggregates exposed to the direct impact of the rain drops, compaction, slaking, particle segregation, and pore filling and clogging by wash-in of fine material. McIntyre (1958a)(1958b) considered the compaction of the soil by raindrops to be the main factor in formation of a soil seal. He suggested that the impact of raindrops destroys the aggregates structure at the soil surface. The released fine material is partly washed into the underlying soil, reducing its porosity. Another part is washed away with the runoff. The compaction of fine material from the destroyed aggregates forms a thin and highly dense layer, which prevents subsequent transport of fine particles to lower levels unless it is removed from the soil surface by turbulent water. He indicated that deposition of fine particles in suspension on the soil surface can contribute to the formation of that thin layer, but he did not consider this contribution to affect significantly the seal permeability. Epstein and Grant (1973) also considered compaction to be the predominant factor in soil surface sealing. However, Eigel and Moore (1983) claimed that seal formation involves rearrangement and realignment of the primary soil particles in the zone of the surface seal, rather than any substantial particle migration. Their conclusion was based mainly on the fact that neither rainfall duration nor intensity was found to have an effect on the particle size distribution of the upper soil layer. Assouline and Mualem (2000) showed that the frequency of their measurements (every 15 min) could have missed most of the dynamic stage of seal formation, so that they observed a seal that was already completely formed. Radcliffe et al. (1991) also concluded that the seal was mainly due to grain rearrangement, and only potentially due to limited pore plugging.
Morin et al. (1981) identified the suction that develops at the seal–soil interface as the cause of seal densification and stabilization. This suction, which becomes efficient for values >400 Pa, continuously arranges the seal from the reserve of clay particles in suspension at the soil surface.
Farres (1978) considered seal formation to be based solely on aggregate breakdown and filling of interstices by smaller aggregates and particles. The aggregate breakdown occurred after its strength was reduced by wetting to a level where the stress imposed by raindrops was sufficient to disrupt the aggregate. Moore (1981b) and Valentin and Ruiz Figueroa (1987) suggested that slaking resulting from rapid wetting could also contribute to aggregate breakdown and porosity reduction.
Bresson and Boiffin (1990), observing a gradual change in pore shape within the seal layer, suggested the gradual coalescence of aggregates due to compaction by raindrops under plastic conditions as an additional process in seal formation. Consequently, Valentin and Bresson (1992) developed a typology to differentiate structural crusts, depending on the processes involved in their formation, including aggregate breakdown by entrapped air implosion, infilling of interaggregate pores by silt particles and/or micro aggregates detached by raindrop impact, deformation of aggregates under plastic state, and agglomeration of dust (<0.5 mm) fragments—slaking, infilling, coalescence, and agglomeration crusts, respectively.
Agassi et al. (1981) suggested that chemical dispersion that depends on the chemistry of the soil water system could play a significant role in seal formation in addition to the physical dispersion caused by the impact of the raindrops.
Obviously, all these different processes are involved in the formation of a disturbed layer at the vicinity of the soil surface, the relative importance of each process depending on the specific conditions that prevail in the soil–rainfall system during seal formation.
|
THE SEAL LAYER AND ITS PROPERTIES
|
|---|
Observed Morphology and Properties of the Seal Layer
Dealing with the morphology of the seal layer raises a problem of terminology. The seal layer that develops at the soil surface during rainfall was characterized, in most of the cases, through observations made after the rainfall ceased and the soil surface dried. This drying process causes the formation of a hard layer at the soil surface, which has introduced the use of the term crust to describe it. However, this "crust" may differ from the "seal" formed during the rainfall event, which rules the rainfall–infiltration relationships during seal formation.
In most of the experimental studies investigating soil sealing under controlled conditions, soil samples were packed in trays or columns, and exposed to simulated rainfall. The physical properties of the simulated rainfall were, in general, reported. However, each rainfall simulator produces its specific rainfall characteristics, which reproduce only partially natural rainfall properties. Also, different soil types at different initial conditions, in terms of grain or aggregate size, packing density, and soil water content, were used. Therefore, the comparison between observations generated by different conditions may not be straightforward.
Different reviews on seal morphology and properties are available (Mualem et al., 1990a; West et al., 1992; Mualem and Assouline, 1996). Various direct methods were applied to study these aspects, most of them assuming a finite observable crust of measurable thickness. These methods include examination of thin sections of the crusted soil surface by means of light microscope, electron microscope, or X-rays radiography and tomography; mechanical analysis of the upper soil layer; and measurements of hydraulic head distribution in the upper soil layer. Indirect methods such as analysis of data from infiltration tests were also applied. The characteristics of soil surface seals developing on bare soils exposed to the direct impact of rainfall seem to be as widely variable as the conditions prevailing during their formation.
McIntyre (1958a)(1958b) provided the crust description that was the most widely adopted: a structure of two discrete and uniform layers, a 0.1-mm skin overlying a 2.0-mm washed-in layer. Some studies reported that only one of these two layers could be observed (Chen et al., 1980; Gal et al., 1984; Bajracharya and Lal, 1999; Wakindiki and Ben-Hur, 2002). However, for more than four decades, the two-layer seal structure is observed (Tackett and Pearson, 1965; Onofiok and Singer, 1984; Wakindiki and Ben-Hur, 2002), and sometimes described in McIntyre's exact terms (Tarchitzky et al., 1984). It is interesting to note that all these studies relied on visual examination of photos from microscopes. When different methods were applied, like porosity or bulk density estimates, the results showed that gradual changes of structure within the crust were more likely to exist (Epstein and Grant, 1973; Eigel and Moore, 1983; Boiffin, 1984; Roth, 1997; Bresson et al., 1998; Fohrer et al., 1999).
A wide range of seal thicknesses obtained under laboratory conditions was reported: 0.1 mm (Chen et al., 1980; Wakindiki and Ben-Hur, 2002), <1.0 mm (Morin et al., 1981; Pagliai et al., 1983; Valentin and Ruiz Figueroa, 1987; Remley and Bradford, 1989), between 1.0 and 3.0 mm (McIntyre, 1958b; Tackett & Pearson, 1965; Tarchitzky et al., 1984; Bajracharya and Lal, 1999; Wakindiki and Ben-Hur, 2002), between 3.0 and 10.0 mm (Sharma et al., 1981; Fohrer et al., 1999), >10.0 mm (Bresson and Boiffin, 1990; Roth, 1997; Bresson et al., 1998). Under field conditions, researchers reported larger thicknesses (up to 20.0 mm) than those observed in the laboratory (Hadas and Frenkel, 1982; Boiffin, 1984).
A 0.1-mm-thick crust, termed skin, was also observed on top of depositional crusts, indicating that it is more likely the result of deposition of fine particles in suspension after rainfall ceases, as observed by Pagliai et al. (1983), rather than the result of "secondary raindrop impact mechanisms" as suggested by Chen et al. (1980). Therefore, this skin layer is nonexistent during rainfall and cannot account for the effects of seal formation on infiltration as they are observed and measured. A quantitative estimation of the role of the skin as part of the seal showed that it is of secondary importance and that it is the underlying seal layer that causes the main impedance to water flow (Mualem et al., 1990b). This result is in fact opposite to the conclusion of McIntyre (1958a), who found that in soils with a stable structure, only the skin layer is effective in decreasing permeability and that the washed-in region is absent or negligible.
Considering the crust bulk density, it is unanimously accepted that the crust layer is more compacted than the soil underneath. Eigel and Moore (1983) and Roth (1997) measured a gradual increase of the bulk density from a depth of 10 mm below soil surface to the crust surface. The maximal value measured by Eigel and Moore (1983) was 1.85 g cm–3. Roth (1997) measured maximal values of 1.78 to 1.91 g cm–3 for soils with silty sand, loamy sand, or sandy loam textures and values of 1.44 to 1.65 g cm–3 for loess soils. Sometimes, there were incompatibilities between the findings resulting from bulk density estimates and those from visual evaluation of the seal thickness. For example, Tackett and Pearson (1965) measured a bulk density of 1.61 g cm–3 for the upper 25 mm of a crusted soil profile. In parallel, they concluded from visual observations that the crust layer was 3 mm thick. Assuming that the soil below the 3-mm crust layer remained undisturbed, the bulk density of the observed crust should have been 3.64 g cm–3, to meet the measured value for the whole upper 25-mm layer, which is impossible as soil particle density is about 2.65 g cm–3. Therefore, the real crust is considerably thicker than the observed one. In addition, some serious doubts about the validity of the identification of the seal layer formed during rainfall with the crust layer visually observed afterwards are raised.
Fohrer et al. (1999) scanned soil crust samples with a medical X-ray tomograph and derived the bulk density distribution with depth in the crust layer. The maximum bulk density was found at 1 mm below soil surface for all samples. The bulk density rose steeply from the surface to the 1-mm depth where the maximum value was measured and then decreased exponentially until it merged with the initial soil bulk density at approximately the 10-mm depth. This result differs significantly from the previous results indicating that the maximum bulk density is reached at the soil surface itself. Also, it is not compatible with a compaction process where the maximal effect is expected to be at the soil surface. The lower bulk density at the upper 1-mm layer can correspond to a loose deposition layer occurring after rainfall has ceased. But it could also be an artifact resulting from the image reconstruction procedure. Fohrer et al. (1999) concluded that whether this result is a real sealing pattern or a systematic technical error cannot be decided from the image itself.
The wide variability found in crust thickness reflects also in the estimates of its saturated hydraulic conductivity. The ratio between the saturated hydraulic conductivities of the seal and of the undisturbed soil beneath has a wide range of values, from 20% (Tackett and Pearson, 1965) to 1% (Hadas and Frenkel, 1982; Simunek et al., 1998), and even to 0.1 or 0.01% (Morin et al., 1981; Perez et al., 1999). In some cases, seal conductance, integrating both thickness and saturated hydraulic conductivity, is addressed. In the study of Baumhardt et al. (1990), the conductance of the upper 10-mm of a sealed soil profile ranged between 0.013 and 0.071 h–1, depending on the applied rainfall intensity.
Conceptual Models of the Seal Layer
Hillel and Gardner (1969)(1970) suggested the simplified conceptual approach that a thin uniform layer can model the seal. They represented the seal by means of a completely saturated thin layer, 5-mm thick, of constant low permeability. This model was widely adopted and applied in a series of studies dealing with infiltration into sealed soils (Ahuja, 1983; Brakensiek and Rawls, 1983; Parlange et al., 1984; Vandevaere et al., 1998).
Mualem and Assouline (1989) introduced a new conceptual model of the seal layer. They suggested that the seal is a nonuniform layer situated at the soil surface. It results from compaction and rearrangement of the soil particles in the disturbed upper zone due to the raindrop impact, and from fine soil particles percolating in-depth during infiltration. Consequently, the seal bulk density, 
c, is the highest at the surface and decreases exponentially with depth, h, to that of the undisturbed soil, :
 | [1] |
where Z is the elevation taken positive upward, 
o is the maximum change in bulk density at the soil surface (h = 0), and 
is a characteristic parameter of the soil–rainfall interaction. Identifying the seal thickness, dc, with the lower boundary of the seal layer at a depth h = dc where the changes in the hydraulic properties are insignificant, namely, where (dc) 
10–3o, the value of is [–ln(10–3)/dc]. The model of Mualem and Assouline (1989) was a theoretical one, which was tested indirectly from infiltration data using an inverse method because of the poor vertical resolution of available bulk density measurements. The validity of this model has been demonstrated directly in recent studies where more accurate methods to measure soil bulk density distribution with depth were applied (Roth, 1997; Bresson et al., 1998).
Studying the relationships between crust bulk density and texture, Roth (1997) suggested that a sigmoidal function could be more appropriate because, once the maximum compaction at the surface has been attained, further drop impact would likely induce maximum compaction at increasing depth. Accordingly, a mathematical expression for c(h) was suggested:
 | [2] |
where 
and v are constants related to the soil–rainfall system. Roth (1997) tested the exponential and the sigmoidal models (Eq. [1] and [2]) by sampling 4-cm2 subsamples from the soil exposed to the rainfall and carefully thinning them to different thicknesses in the range of 2 to 12 mm. The bulk densities for each thickness were measured using an immersion method, and the bulk density depth functions derived accordingly. He found that both models showed a good to very good fit to measured data, the exponential model appearing to represent best the initial stages of crust formation, and the sigmoidal model, the later stages.
The exponential model of Mualem and Assouline (1989) can be extended to include the case where maximum compaction extends in-depth. The mathematical expression of this new model is
 | [3] |
where ß and n are constants related to the soil–rainfall system. A comparison between the three models (Eq. [1]–[3]) is depicted in Fig. 1
, based on the experimental data of Roth (1997) for sandy silt (Fig. 1a) and loamy silt (Fig. 1b) soils. The corresponding fitted parameters are shown in Table 1. From a practical point of view, the performances of the three expressions are similar. Because of the limits of the experimental method, there are no bulk density data for the upper millimeter, where the differences between the models concentrate. Therefore, one can find here an advantage for the new model suggested (Eq. [3]), as it is flexible enough to cover the trends of both the exponential model (Eq. [1]) and the sigmoidal one (Eq. [2]).
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 1. The modeled distribution with depth of the bulk density within the seal layer according to the models in Eq. [1], [2], and [3], compared with the experimental data of Roth (1997) for (a) a sandy silt and (b) a loamy silt soil.
|
|
Modeling the Seal Hydraulic Properties
The main purpose for seal modeling was the simulation of infiltration into crust-topped profiles or during seal formation. Since measuring seal properties appeared to be rather difficult because of its small thickness and the lack of suitable experimental procedures for its identification, these were usually simulated using various degrees of arbitrary approximation.
The simplest representation of the seal properties was to assume a saturated layer with constant hydraulic conductivity and thickness (Hillel and Gardner, 1970; Ahuja and Swartzendruber, 1973). Sometimes, the seal conductance was applied to avoid the need to assume a hypothetical thickness (Farrell and Larson, 1972; Ahuja and Ross, 1983). But characterizing the seal layer by its saturated hydraulic conductivity is not sufficient to deal with redistribution and drying processes, where information on the seal hydraulic properties is compulsory. Therefore, other studies adopted the same conceptual approach, representing the seal as a uniform layer, but dropped the limiting assumption of saturation and characterized the seal layer by means of approximated hydraulic properties (i.e., water retention curve and hydraulic conductivity function) (Moore, 1981a; Philip, 1998; Smith et al., 1999). Moore (1981b) assumed that the hydraulic properties of the seal are identical to those of the undisturbed soil and differ only in the value of the saturated hydraulic conductivity. Philip (1998) presented illustrative crusts of various thicknesses; the crust properties differed from those of the soil by means of porosity and internal length scale. Smith et al. (1999) assumed that the hydraulic characteristics of the seal differed from those of the underlying soil only by their lower saturated hydraulic conductivities and air-entry values. However, in all these cases, no real soil crusting was considered but rather hypothetical two-layer flow systems with the so-called crust having prescribed properties arbitrarily chosen.
Two attempts were made to derive the seal hydraulic properties from those of the undisturbed soil below (Mualem and Assouline, 1989; Baumhardt et al., 1990). The study of Baumhardt et al. (1990) maintained the conceptual model of the seal as a uniform layer. The method suggested in that study is depicted in Fig. 2 . The seal water retention curve was determined through the following steps: (i) estimate experimentally the seal conductance and evaluate its saturated hydraulic conductivity, assuming a seal thickness of 5 mm; (ii) based on the estimated saturated hydraulic conductivity, compute its saturated water content, 
sc, using the Kozeny–Carman relationship, and its water-entry value, hwevs, using the Poiseuille equation; and (iii) determine the Brooks and Corey (1964) exponent 
c for the seal layer graphically, assuming a proportional increase in the initial water content of the seal, ic, with bulk density.
Mualem and Assouline (1989) considered the seal a nonuniform layer, where the bulk density distributes with depth from a maximum value at the soil surface to the value representing the undisturbed initial soil. The distributed hydraulic properties of the seal layer were related to the bulk density distribution, and represented by two functions, the retention curve, c(
,c) and the conductivity function, Kc(,c). The Brooks and Corey relationship (1964) for the retention curve, and Mualem's model (1976) for the conductivity function were chosen to express the seal hydraulic properties:
 | [4] |
 | [5] |
where sc and rc are the saturated and the residual water content, c is the pore size distribution parameter, ac is the air-entry value, and Ksc is the saturated hydraulic conductivity of the seal layer. All of these parameters were related to the corresponding parameters of the undisturbed soil and to the bulk density distribution, c(h), of the seal layer (Eq. [1]) according to the relationships:
 | [6] |
 | [7] |
 | [8] |
 | [9] |
 | [10] |
where s is the saturated water content, r is the residual water content, is the Brooks and Corey parameter, a is the air-entry value, Ks is the saturated hydraulic conductivity, and and s are the bulk and the solids density of the undisturbed soil. Equation [6] assumes that the saturated water content is equal to porosity. Equation [7] assumes that the residual water content on a weight basis remains unchanged during compaction. Equation [8] is an empirical relationship based on experimental data related to the effect of compaction on a. Equation [10] is obtained by applying principles similar to those suggested by Mualem (1986). The linear relationship between c and c (Eq. [9]), with the fitting parameter C, is the simplest approximation that can be assumed since data upon which a more complex relationship can rely are not available. To illustrate the results of this model, the hydraulic functions at different depths within the seal layer developed on the Atwood silty clay loam (fine-silty, mixed, semiactive, thermic Typic Paleudalf) (Baumhardt et al., 1990) are shown in Fig. 3
. The calibration of the model for that soil was presented in details in Assouline and Mualem (1997). As the sealed soil surface is approached, s, , and Ks decrease, and r, a, d/d, and dK/d increase, according to the relationships in Eq. [6] through [10].
|
THE DYNAMIC STAGE OF SEAL FORMATION
|
|---|
Two processes occur simultaneously during the collision of a raindrop at a bare soil surface: (i) a transfer of kinetic energy from the raindrop to the soil surface and (ii) soil imbibition. The transfer of energy results in mechanical changes of the soil surface, expressed in terms of compaction, particle detachment, and splash. For a given initial condition of the soil surface, the amount of transferred energy determines the magnitude of these changes. Simultaneously, the soil imbibition (i) facilitates particle detachment, the degree of which is closely related to the chemical conditions of the soil-water system; (ii) enhances the collapse of unsaturated soil aggregates, facilitated by the sudden compression of entrapped air within the aggregates; and (iii) decreases the soil resistance to destruction under subsequent striking raindrops. Loose soil particles resulting from detachment, splash, and aggregate collapse can be transported, at least at the early stages of formation when the infiltration rate is still high, by the infiltrating water within the upper soil layer and fill the interaggregate voids of the undisturbed soil below the soil surface. Also, as the exposure to rainfall continues, further compaction by raindrops and changes in the deeper aggregates once wetted (swelling, collapse) can occur (Bresson and Boiffin, 1990; Bresson and Cadot, 1992). These processes have two consequences: an additional increase of the soil bulk density and the thickening of the affected zone, thus creating a disturbed layer at the vicinity of the soil surface.
Factors Affecting Seal Formation
The extensive experimental work, under laboratory and field conditions, performed to investigate rainfall-induced soil surface sealing, has revealed a large number of soil and rainfall properties that affect this phenomenon. Reviews on physical (Bradford and Huang, 1992) and chemical and mineralogical components (Shainberg, 1992) affecting seal formation are available. A concise presentation of the main factors is presented in the following for the sake of completeness.
The physical rainfall characteristics (i.e., intensity and kinetic energy) were found to play a major role in determining the seal properties and the rate of seal formation (Morin and Benyamini, 1977; Romkens et al., 1986; Mohamed and Kohl, 1987; Baumhardt et al., 1990; Betzalel et al., 1995). Lower rainfall intensities led to lower final infiltration rates (Morin and Benyamini, 1977; Baumhardt et al., 1990; Freebairn et al., 1991). Betzalel et al., (1995) found that both the rate of decline in water intake and the final infiltration rate are dependent on the rainfall kinetic energy. Mohamed and Kohl (1987) concluded that a low kinetic energy rainfall causes a thinner and less compacted seal. However, Romkens et al. (1986) found that both rainfall intensity and kinetic energy affected the seal hydraulic conductivity.
In addition to the rainfall characteristics, a wide range of soil properties also affects soil sealing. The common aspect of these effects is the cohesive power between the soil particles, which is an expression of the soil resistance to the destructive action of rainfall. Experimentally, this is supported by the fact that aggregate destruction and seal formation were found to be affected by (i) soil mineralogy and texture (Mannering, 1967; Tarchitzky et al., 1984; Ben-Hur et al., 1985a; Stern et al., 1991; Bradford and Huang, 1993; Wakindiki and Ben-Hur, 2002), (ii) aggregate size distribution (Moldenhauer and Kemper, 1969; Farres, 1978; Luk, 1985; Freebairn et al., 1991), (iii) aggregate stability (Bryan, 1969; Le Bissonnais, 1996), (iv) initial bulk density (Boiffin, 1984; Luk, 1985), (v) initial water content distribution in the soil profile (Francis and Cruse, 1983; Ben-Hur et al., 1985b; Truman et al., 1990; Le Bissonnais and Singer, 1992), (vi) application of phosphogypsum or polymers to the upper soil layer (Helalia et al., 1988; Shainberg et al., 1990; Agassi et al., 1996; Ben-Hur, 2001), (vii) slope (Poesen, 1987; Bradford and Huang, 1993; Fox et al., 1997), (viii) organic matter content (Le Bissonnais and Arrouays, 1997), (ix) soil exchangeable Na percentage (ESP) (Agassi et al., 1981; Kazman et al., 1983; Agassi et al., 1996), and (x) electrical conductivity (EC) of the applied water (Agassi et al., 1981; Hadas and Frenkel, 1982; Shainberg et al., 1990; Agassi et al., 1996). It is thus evident that the number of factors affecting the rate of the seal formation, its properties, or its effects on infiltration is large. The knowledge gained by the intensive experimental effort invested to date yielded highly valuable results. The main factors involved in seal formation were identified, and a large amount of data on soil seal characteristics and on infiltration under sealing conditions was produced. However, considering the complexity of the interactions between the different factors, it is very difficult or rather impossible to define the respective modus operandi of each factor from the qualitative analysis of the experimental observations alone. That can be done only by means of physically based quantitative theoretical models, which allow us to compare results from one experiment with those from another, predict the response of soil-rainfall systems under different conditions from those experimentally studied, and evaluate the relative effect of each factor.
Seal Formation: Observations and Modeling
The study of the dynamics of seal formation relied on both direct observations, mainly of the morphological aspects of the seal layer (Tackett and Pearson, 1965; Edwards and Larson, 1969; Farres, 1978; Boiffin, 1984; Bresson and Boiffin, 1990; Bajracharya and Lal, 1999; Fohrer et al., 1999) and on indirect measurements such as infiltration experiments (Morin and Benyamini, 1977; Romkens et al., 1986; Baumhardt et al., 1990; Bajracharya and Lal, 1999).
Tackett and Pearson (1965) measured a gradual increase in the bulk density of the upper 25 mm of the soil profile during rainfall, from 1.32 to 1.61 g cm–3. The rate of growth of the bulk density decreased exponentially, the maximum compaction being attained after 50 mm of cumulative rainfall. Similarly, Boiffin (1984) measured a hyperbolic decrease of the void ratio of the upper soil layer with rainfall, at a rate that was higher as the initial soil porosity was larger. Farres (1978) found that the observed crust thickness, dc, increases with the cumulative rainfall, R, and suggested the empirical relationship:
 | [11] |
where a and b are parameters related to the aggregate size of the soil.
Edwards and Larson (1969) measured an exponential decrease of the saturated hydraulic conductivity of 5-mm-thick soil samples during the exposure to simulated rainfall. Adopting this result, a number of studies modeled seal formation by means of the expression:
 | [12] |
where Kc, Ks, and Kf, denote the saturated hydraulic conductivity of the seal (crust) during the dynamic stage, the initial soil, and the completely formed seal, respectively, and is an empirical constant of the soil–rainfall system. The independent variable x was rainfall duration, t (Moore, 1981b; Ahuja, 1983) or cumulative kinetic energy, E (Brakensiek and Rawls, 1983; Chu et al., 1986).
Baumhardt et al. (1990) described the relationship between the conductance, B, of the upper 10-mm Atwood soil, and E using a three-stage equation:
 | [13] |
where Ep and Et are the rainfall cumulative kinetic energy until ponding and until the terminal conductance value, Bt, is reached, and a1 to a5 are empirical parameters for the soil–rainfall system. Bt was found to be dependent only on the rainstorm intensity, I, in agreement with the findings of Romkens et al. (1986).
Mualem et al. (1990c) suggested a dynamic model of seal formation based on the conceptual model of the seal layer presented in Eq. [1]. The maximum increase in the soil bulk density at the soil surface, o, and the seal thickness, dc, were considered to be similar functions of E:
 | [14] |
 | [15] |
where *o and d*c are the maximal values reached after a long exposure to rainfall, and 
and 
are soil–rainfall characteristics. Consequently, the distribution of c with depth (Eq. [1]) becomes a function of E:
 | [16] |
with (E) = [–ln(10–3)/dc(E)] and is the undisturbed soil bulk density. The seal hydraulic properties described using Eq. ([4]–[10]) are thus also dynamic functions of the cumulative rainfall kinetic energy, E. The model was calibrated and checked using data from Morin et al. (1981) (Mualem et al., 1990c; Mualem and Assouline, 1991). The trends in the dynamic changes in the upper layer bulk density, and in the seal thickness and saturated hydraulic conductivity were in agreement with the observations of Tackett and Pearson (1965), Farres (1978), and Edwards and Larson (1969).
Assouline and Mualem (1997) developed a dynamic model that relates the formation of a seal at the surface of a bare soil exposed to water drops impacts to the initial soil mechanical and hydraulic properties as well as the physical characteristics of the regional rainfall or the applied irrigation. At the soil surface, the disturbance resulting from raindrop impact and imbibition is expressed in terms of the bulk density increase, o, taken as a function of rainfall intensity, I, and time of exposure to rainfall, t:
 | [17] |
where *o is as defined above, and the soil–rainfall characteristic, 
, is defined by
 | [18] |
where f(d,I) is the raindrop size distribution, dmax is the maximal drop diameter, k is a rainfall parameter interrelating raindrop velocity and drop diameter d, 
is a fitting parameter, and 
(i, i) is the initial soil shear strength. During the early stage of seal formation, the concentration of the fine particles in the percolating suspension is still low and the open pore structure only slightly affected. Therefore, the soil hydraulic conductivity remains relatively high and allows high infiltration rates with high transport capacity of fine soil particles in-depth. Consequently, for a given soil–rainfall system, the maximal thickness of the disturbed layer should be reached shortly after the beginning of rainfall. To simplify the model, a constant dcmax value was assumed during the whole rainfall period. However, as the rainfall kinetic energy producing the fine particles in suspension and the infiltration rates transporting them in-depth are both functions of the rainfall intensity, I, one may expect dcmax and, consequently, , to be dependent on I as well. The distribution of c with depth (Eq. [1)] becomes
 | [19] |
with (I) = [–ln(10–3)/dcmax(I)]. Consequently, the seal hydraulic properties described using Eq. [4] to [10] became also dynamic functions of I, h, and t. The model was calibrated and validated using infiltration data from Morin et al. (1981) and Baumhardt (1985) and performed well in predicting infiltration curves for a wide range of conditions different from those of the calibration case (Assouline and Mualem, 1997). An illustration of the functions c(E,h) and Kc(E,h) characterizing the dynamics of seal formation is depicted in Fig. 4
. It is interesting to note that the theoretical development leading to the definition of (Eq. [18]) has revealed that both E and I play a role in determining the seal properties, as pointed out by the statistical analysis of Romkens et al. (1986).
The assumption of a constant maximal seal thickness during the whole dynamic stage is apparently in disagreement with the models of Farres (1978) (Eq. [11]) and Mualem et al. (1990c) (Eq. [12]), and in accord with the conclusion of Eigel and Moore (1983) that the seal thickness is independent of E. The seal thickness was generally determined experimentally by optical means (light or electron microscopy), and the resolution of the observations was not specified. Recently, Roth (1997) applied the immersion method to measure soil bulk density and assumed that the seal thickness is at the depth where c(h)/ = 1.01. The relative error in estimating changes in the soil bulk density in the X-ray radiography procedure of Bresson et al. (1998) is 0.015. The resolution, 
, of the soil bulk density measurement method applied determines the minimal observable change in bulk density, min = . The seal will be identified with the soil layer where c > + min. Therefore, the thickness of the "observable" seal layer, dobs, will be related to as it will be the depth h = dobs where c(h,t) = min+ = (1+ ). According to Eq. [19], and for a given rainfall intensity, it is possible to express the observable seal thickness as a function of time, dobs(t) for every :
 | [20] |
To illustrate this point, Eq. [20] was applied to the case of the dynamic model calibrated to the data of Morin et al. (1981) for the Sharon (coarse-silty, mixed, active, acid, mesic Oxyaquic Udifluvents) sandy loam soil and I = 70 mm h–1 (Assouline and Mualem, 1997, 2000). The resulting dobs(t) functions, expressed in terms of cumulative rainfall, are presented in Fig. 5
, for three values (0.015, 0.025, and 0.05), in comparison with the modeled dc(t) function (corresponding to = 0.001), and the calibrated dcmax value (12.3 mm). It clearly indicates that the observed seal thickness is time dependent and that this relationship is dominated by the resolution of the soil bulk density measurement method. At the initial stages, the seal is hardly observable, and then it quickly reaches (within 35 mm or 30 min of rainfall) almost its final thickness. As the resolution of the applied measurement technique is higher (lower values), the seal layer is observable earlier and its thickness increases faster toward a higher final value. The rapid increase of dobs(t) within the first 15 min of exposure to rainfall may explain why Eigel and Moore (1983), observing the seal every 15 min from the beginning of rainfall, did not detect any increase of the seal thickness and concluded that it was independent of the cumulative rainfall kinetic energy. The dynamic evolution of the seal layer observed by Farres (1978) for the case of the large aggregates is also depicted in Fig. 5. Although the dynamic model was not calibrated for that case, the observed trend is well represented by the model, considering that 0.015 < < 0.025 could be a fair estimate of the resolution characterizing Farres method (photomicrograph analysis of thin sections).
View larger version (23K):
[in this window]
[in a new window]
|
Fig. 5. The dynamic evolution of the observable seal thickness as calculated for three different resolution levels , compared with the modeled seal thickness (solid line) and the experimental data of Farres (1978).
|
|
The exponential model of Mualem and Assouline (1989), with the condition for the seal limit at /o = 10–3, leads to final seal thicknesses that are much larger than the common values found in the literature. The observable final thicknesses calculated when the resolution of the observation method is accounted for, as presented in Fig. 5, are very close to those reported in the literature, indicating that the contradiction between the reported seal thicknesses and the values resulting from the exponential model is only an apparent one. The observed thin seals in fact represent much thicker disturbed layers, where the unobserved remaining disturbed part can still affect the hydraulic properties of the initial soil and, consequently, the flow processes in the upper soil layer.
|
MODELING FLOW PROCESSES IN SEALED SOILS AND DURING SEAL FORMATION
|
|---|
Sealed soils can be considered as a special case of layered soils. Therefore, pioneering studies dealing with flow processes, and mainly, infiltration, into layered soils (Colman and Bodman, 1945; Hanks and Bowers, 1962; Philip, 1967; Childs, 1969) are pertinent and relevant as far as these processes are considered in sealed soil profiles. However, the basic difference between simulating flow processes in a layered soil and in a sealed one is the fact that, in the latter, the upper soil layer cannot be characterized by an artificial layer with arbitrarily chosen properties. As has been shown in the previous sections, the seal physical and hydraulic properties stem from the hydraulic properties of the soil upon which it develops as well as from the specific soil and rainfall conditions prevailing during formation. Therefore, this should be the key point in differentiating between studies dealing with layered soils or sealed ones.
Infiltration in Sealed Soils and During Seal Formation
A review of the progress in soil seal characterization and representation to model infiltration into sealed soils can be found in Ahuja and Swartzendruber (1992) and Mualem and Assouline (1992)(1996). The problem of calculating the infiltration rate in crusted soils was first addressed by Hillel and Gardner (1969)(1970), who presumed that a sealed, or crusted, soil can be modeled as a uniform soil profile capped with a saturated thin layer of low permeability and prescribed constant properties. Their simplified solution was based on the Green and Ampt (1911) approach and assumed constant water content (or suction) at the interface between the seal and the soil beneath. Variations and extensions of this basic approach were suggested in a series of studies on infiltration through sealed soils with constant seal properties (Ahuja, 1974, 1983; Moore, 1981a; Parlange et al., 1984). This approach was also applied to simulate infiltration during the dynamic stage of seal formation. In this case, time-dependent seal hydraulic conductivity functions were incorporated in the models (Farrell and Larson, 1972; Whisler et al., 1979; Moore, 1981b; Ahuja, 1983; Brakensiek and Rawls, 1983; Chu et al., 1986). These studies considered the sealed soil as a special case of a layered profile and the seal as a thin layer, often prescribed as 0.5 cm thick, that saturates very quickly if not instantaneously. Recently, Vandervaere et al. (1998) applied the Green and Ampt model, assuming that the wetting front potential decreases suddenly as the wetting front leaves the seal and enters the soil. The result is a discontinuous drop in the infiltration rate at that moment. Philip (1998), analyzing the results of some of the studies that applied the Green and Ampt model, concluded: "Green-Ampt is ill-fitted to the analysis of infiltration into crusted soils, and there is no convincing way of patching it" (Philip, 1998, p. 1926). He suggested the flux-concentration method (Philip, 1973) to solve accurately the nonlinear problem of ponded infiltration into crusted soils while avoiding the simplifications and approximations that characterized the previous studies. An additional conceptual model was suggested by Smith et al. (1999), who extended the model of Corradini et al. (1997), which corresponds to infiltration into homogeneous soils regardless of surface saturation, to represent infiltration into layered or crusted soil profiles. However, the seal layer properties were still arbitrarily chosen and assumed to vary from those of the initial undisturbed soil by only the air-entry value and the saturated hydraulic conductivity.
A different approach was developed by Baumhardt et al. (1990), Mualem et al., (1993), and Assouline and Mualem (1997). They addressed the problem of infiltration into sealed soils or during seal formation, while adjoining to the seal layer complete hydraulic properties related to the characteristics of the specific soil–rainfall system. These studies solve numerically the one-dimensional water flow equation, also known as Richards' equation. The model of Baumhardt et al. (1990) considers the seal a homogeneous layer and still maintains the assumption that its thickness is 0.5 cm, as did earlier studies. The performance of this model is illustrated in Fig. 6
, which depicts predicted cumulative infiltration curves during soil surface sealing for two rainfall intensities. The other solutions apply the model of Mualem and Assouline (1989), who considered the seal a nonuniform layer where the bulk density distributes exponentially with depth, although an equivalent homogeneous layer still can represent technically the conceptual nonuniform seal if and when necessary (Mualem et al., 1993). In the following, infiltration into sealed soils and during seal formation will be illustrated using this approach, as it will allow also evaluation of the effect of assuming a uniform seal layer instead of the nonuniform one.
Solution of the Flow Equation in Sealed Soils and During Seal Formation
The solution of the flow problem relies on the numerical solution of Richards' equation governing transient, one-dimensional vertical water flow in saturated–unsaturated soil for specific initial and boundary conditions:
 | [21] |
where t is time, Z is the vertical coordinate taken positive upwards, () is the water retention curve relating water content to capillary head, and K() is the hydraulic conductivity function relating hydraulic conductivity to capillary head. In principle, during seal formation, changes with t, and consequently, depends on both and . Therefore, the expansion of 
/t should include an additional term, (/)/(/t). This second term adds significant complexity to the numerical solution because of the necessity of dealing with a dynamically shrinking grid, which is associated with the bulk density change. This aspect is not addressed in soil seal studies because the seal is rather thin, and the resulting effect of the shrinking matrix domain on the water mass balance is negligible. Even if a much thicker disturbed layer is considered (Mualem and Assouline, 1989), the big change in bulk density is limited to the upper 2 to 3 mm. Therefore, the associated change that could result in the translocation of the soil surface and the change in the solid matrix domain are thus practically insignificant. In other words, the second term(/)/(/t) can be neglected as far as the water mass balance is concerned, and Richards' equation in its simplified form of Eq. [21] can be applied. The finite difference scheme resulting from the Crank–Nicholson approximation (Crank, 1975) has been applied to solve Eq. [21] numerically.
Two different approaches can be followed in applying such a solution to flow processes in sealed soils or during seal formation. One is to consider the flow system a continuous nonuniform soil profile. In this system, there is no interface between the seal layer and the undisturbed soil, as their hydraulic properties merge below Z = –dc. However, there is a gradual but steep change in the seal properties at the vicinity of the soil surface as specific water retention curves, (), and hydraulic conductivity functions, K(), distribute with depth, h, according to c(h), for the final seal, or to c(I,h,t), for the dynamic seal (Eq. [4]–[ 10] and Eq. [19]). The other approach is to replace the nonuniform seal by a uniform equivalent layer, thereby generating a two–uniform layers flow system. For the undisturbed soil (Z –dc), the flow equation (Eq. [21]) is solved for the hydraulic properties of the undisturbed soil. For the uniform seal layer (0 > Z > –dc), it is solved for the equivalent mean functions, c() and Kc(), for which the mean seal hydraulic parameters are computed according to Eq. [6] to [9], where 
= –d–1c
dc0c
dh replaces c(h) and the mean saturated seal hydraulic conductivity, Ksc, is evaluated as the harmonic average of Ksc(h) (Eq. [10]) over the disturbed layer, 
sc = dc
–1. In this case, there is a discontinuity in the hydraulic properties at Z = –dc.
The grid used for the finite difference scheme was irregular. In the upper 8-cm layer, the nodal spacing, Z, was 0.2 cm, while it was 2.0 cm for all the remaining soil profile. To reduce numerical errors, the fine grid domain was set thicker than the nonuniform seal layer, and stretched deeper than the interface between the equivalent seal layer and the underlying undisturbed soil in the uniform seal case, so that it included all the steep changes in hydraulic properties. For that same reason, the internode hydraulic conductivity of the nonuniform seal was evaluated by means of the harmonic average while it was evaluated by means of the arithmetic average when the soil layers were homogeneous.
The flow equation for the different approaches was solved using the following conditions. The lower boundary condition of the flow system, at depth L, was characterized by a constant hydraulic head, L, low enough so that the hydraulic conductivity, K(L), and the downward flux, q(L), can be neglected:
 | [22] |
The upper boundary condition switched from a Neuman condition to a Dirichlet condition during the wetting process. Before ponding, while Z=0 < 0, the upper boundary condition was
 | [23] |
with I being the applied rainfall intensity. At ponding time, t = tp, the boundary condition switched to
 | [24] |
with tp being the ponding time. The initial condition was assumed to be
 | [25] |
This solution, for the uniform equivalent seal layer, was applied to a large range of problems related to flow processes during the dynamic stage of soil surface sealing and when the seal is completely established. Detailed presentation of the results and discussions are available in the literature (Mualem et al., 1993; Assouline and Mualem, 1997, 2000, 2003). For the simulation of infiltration during the dynamic stage of seal formation, the solution was calibrated to the data of Baumhardt (1985), the experimental setup consisting of 30-cm-long perforated soil columns packed with dry Atwood silty clay loam. Large ranges of simulated rainfall intensity (20, 30, 40, 60, and 90 mm h–1), kinetic energy density (0, 11.4, 20, and 27.5 J m–2 mm–1), and duration (25, 60, and 120 min) were applied. The rainfall intensity and kinetic energy were independent variables; the effect of each of these two variables was studied while keeping the other constant. The solution was calibrated for the (60 mm h–1; 27.5 J m–2 mm–1; 120 min) case. The performance of the dynamic model and the numerical solution for the validation phase is shown in Fig. 7
. The predicted infiltration curves during soil sealing at rainfall intensities of 30 and 90 mm h–1 are in good agreement with the measured data, indicating that the model account for the main factors affecting seal formation.
The Effect of Assuming a Uniform Seal Layer
The solution of the flow equation (Eq. [21]) for the uniform equivalent seal layer case was compared with that resulting from the nonuniform exponential seal layer one to estimate the effect of the seal representation on infiltration (Assouline and Mualem, 2001). For the Atwood soil case, the bulk density of the disturbed layer, c, decreases from 1.8 g cm–3 at Z = 0 to 1.4 g cm–3 at Z = –dc, with dc = 3.6 cm. Its corresponding saturated water content, sc(Z) is also exponentially distributed with depth and increases from 0.27 cm3 cm–3 at Z = 0 to 0.42 cm3 cm–3 at Z = –dc. The saturated hydraulic conductivity, Ksc, is a complex nonlinear expression involving the relationships assumed between c and the hydraulic parameters of the soil (Eq. [10]). The resulting Ksc(Z) is a sigmoid function varying from 0.014 mm h–1 at Z = 0 to 7.0 mm h–1 at Z = –dc. In the case of the uniform seal, the constant values for c and sc are the arithmetic average of the corresponding distribution functions c(Z) and sc(Z) within the seal domain, being 1.46 g cm–3 and 0.40 cm3 cm–3, respectively. Ksc is the harmonic mean of Ksc(Z) over the seal thickness, and is equal to 0.42 mm h–1. In this case, the differences between the values corresponding to the soil surface and the bottom of the seal are reduced considerably in comparison with the exponential seal case.
During the dynamic stage of seal formation, the approach adopted to represent the seal layer in the flow system has only little effect on the infiltration curve. The exponential seal representation increases the ponding time and consequently the infiltration rates just after ponding compared with the uniform approximation. However, the difference after a long exposure to rainfall is very small. This result indicates that the dynamic model parameters, determined by a calibration procedure assuming an equivalent uniform layer, are still valid and can be adequate for predicting the infiltration function when the seal is in fact nonuniform. However, a significant difference between the two approaches exists in terms of the soil water content distribution with depth, and especially at the soil surface. The dynamics of the water content changes with time at the soil surface are depicted in Fig. 8
for two rainfall intensities. The curves in this figure illustrate the interaction between the decreasing porosity, on one hand, and the increasing water content, on the other hand, during rainfall and seal formation, as calculated in the two approaches. In both cases, the initial porosity and the initial water content at the soil surface, at the beginning of rainfall, are 0.42 and 0.275 cm3 cm–3, respectively. During the exposure to rainfall, the porosity decreases with time, at a rate that is a function of the soil–rainfall parameter (Eq. [18]), and the water content increases with time, at a rate that depends on the rainfall intensity. At the point at which the water content of the seal surface reaches saturation, it reverses its trend and decreases as it follows the decreasing seal porosity. This pattern is common to both approaches but with big differences in water content values. Applying the uniform seal approximation significantly diminishes the amplitude of the changes in the saturated water content with time and with the rainfall intensity. The reason is that it tends to the mean porosity of the equivalent seal for the uniform seal approximation, while it tends to the much lower minimal porosity reached at the soil surface for the nonuniform seal case.
TOP
ABSTRACT
INTRODUCTION
