VZJ sign up for etocs
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


Published online 20 November 2006
Published in Vadose Zone J 5:1246-1256 (2006)
DOI: 10.2136/vzj2005.0144
© 2006 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 Google Scholar
Google Scholar
Right arrow Articles by Fujimaki, H.
Right arrow Articles by Nakane, K.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Fujimaki, H.
Right arrow Articles by Nakane, K.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Fujimaki, H.
Right arrow Articles by Nakane, K.
Related Collections
Right arrow Soil Salinity
Right arrow Evapotranspiration Models
Right arrow Solute Transport Models

ORIGINAL RESEARCH

Effect of a Salt Crust on Evaporation from a Bare Saline Soil

Haruyuki Fujimakia,*, Takahiro Shimanoa, Mitsuhiro Inoueb and Kazurou Nakanec

a Univ. of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8572, Japan
b Arid Land Research Center, Tottori Univ., 1390 Hamasaka, Tottori 680-8550, Japan
c National Research Inst. for Earth Science and Disaster Prevention

* Corresponding author (fujimaki{at}sakura.cc.tsukuba.ac.jp)

Received 9 December 2005.



    ABSTRACT
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 THEORY
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Evaporation from the soil surface is a major cause of the salinization of irrigated soils in arid and semiarid regions. To optimize irrigation scheduling under saline conditions, it is essential to be able to accurately predict solute transport in soils and the evaporation rate. We conducted laboratory column experiments under constant meteorological conditions, except for radiation, which was automatically regulated such that the temperature of the soil remained the same as that of the air. The concentration of the initial soil solution and that of the inflowing water from the bottom of the 5.2-cm-long column were set at 3000 g m–3. The evaporation experiments were performed with three combinations of soil and solute. Although the soil surface was kept wet by maintaining a low suction at the bottom, the evaporation rate was found to decrease considerably with time. This decrease could not be explained by a decrease in osmotic potential alone, but rather was due also to the formation of a salt crust near the surface. The bulk transfer equation for evaporation was therefore modified to include a resistance to water vapor diffusion caused by the salt crust. The dependence of the salt crust resistance on the amount of accumulated salt was evaluated experimentally and theoretically. In our numerical analysis, we used independently estimated soil hydraulic and solute transport parameters. Results show that the convection–dispersion equation (CDE) tends to overestimate backward diffusion near an evaporating soil surface, thus significantly delaying salt accumulation at the soil surface and decreasing the evaporation rate. Since the CDE uses an analogy of Fick's law to describe mechanical dispersion, the dispersion term overestimated the downward transport of solutes against upward convective transport.

Abbreviations: CDE, convection–dispersion equation


    INTRODUCTION
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 THEORY
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
IRRIGATION-INDUCED SALINIZATION is a primary cause of reduced agricultural productivity in arid and semiarid lands. To optimize irrigation scheduling under saline conditions, it is essential to be able to accurately predict solute transport in soils. Of particular importance is the prediction of the evaporation rate from saline soil surfaces, since this is the main driving force for salt accumulation. Conversely, accumulated salts near the soil surface may affect the evaporation rate. Although the mechanisms of soil salinity affecting the evaporation rate have not been studied in great detail, three key factors appear to be involved: (i) a reduced osmotic potential at the soil surface, (ii) increased albedo, and (iii) increased resistance to water vapor diffusion due to the formation of a salt crust. While recently several models have been developed that consider the effects of osmotic potential on evaporation from a bare saline soil (e.g., Noborio et al., 1996; Yakirevich et al., 1997), these models do not consider the effects induced by a salt crust that is often present in arid lands with shallow groundwater tables.

Frequent and uniform irrigation with a sprinkler system may impede the precipitation or crystallization of excess salts on the soil surface. Continuous salt accumulation at the wetted soil surface is unavoidable, however, for furrow irrigation or drip irrigation with saline water (Ayers and Westcot, 1985; Abrol et al., 1988; Mmolawa and Or, 2000). Salt buildup hence often leads to the formation of a salt crust on the soil surface. Fujimaki et al. (2003) investigated the effect of a salt crust on soil albedo and presented an empirical equation describing the dependence of albedo on both water content and salinity. A salt crust would also impede water vapor transport into the atmosphere by blocking air flow much like a straw or gravel mulch.

To date, few studies have been performed on the resistance of a salt crust to water vapor transfer. Based on both field and laboratory studies, Chen (1992) reported that a very thin (several millimeters) salt crust can reduce evaporation to only a few percentage points of the potential evaporation rate, even when the underlying soil is saturated, but did not mathematically quantify the relationship between the degree of evaporation and salinity. Shimojima et al. (1996) investigated the effects of solute crystals formed in the dry surface layer of a soil on evaporation. Their experiments were limited to sand with a dry surface layer and very low evaporative conditions.

To accurately predict the relative humidity at the evaporating soil surface and the extent of the salt crust being created, it is essential that solute transport near the soil surface be predicted precisely. While several studies have analyzed concentration profiles that developed during evaporative conditions using the CDE (Todd and Kemper, 1972; Elrick et al., 1994; Shimojima et al., 1996), few tested results using independently measured values of the mechanical dispersion coefficient, Dm. The CDE is valid only after a certain average travel distance has been attained (Jury et al., 1991). The applicability of the CDE to areas near evaporating soil surfaces is questionable, since these areas constitute a no-flow boundary condition for the upward moving solutes. Todd and Kemper (1972) proposed a method to determine Dm from a nearly steady-state concentration profile having a steep gradient near the evaporating soil surface; however, they did not compare Dm with results from a downward flow experiment conducted at a similar water content. If the CDE is to be valid near an evaporating soil surface, Dm obtained by their method must agree with the value obtained from a downward flow experiment. Elrick et al. (1994) presented analytical solutions for salt accumulation applicable to soil profiles having a water content gradient. Instead of using independently measured Dm values, they compared the analytical solutions with measured concentration profiles during upward flow using best-fit values for Dm. In a related study, Shimojima et al. (1996) neglected Dm because of the very low evaporation rate in their experiments.

The objectives of this study were (i) to quantify the resistance to evaporation caused by salt crusts and (ii) to test the applicability of the CDE to transport conditions near evaporating soil surfaces. To achieve this, we used independently obtained Dm values.


    THEORY
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 THEORY
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
The evaporation rate, E (g cm–2 s–1), from a soil surface may be described by the bulk transfer equation (Daamen and Simmonds, 1996; Noborio et al., 1996; Yakirevich et al., 1997):

Formula 1[1]
where {rho}v* is the saturated water vapor density (g cm–3), hr is the relative humidity, ra is the aerodynamic resistance (s cm–1), and where the subscripts s and a denote the soil surface and air at the reference height, respectively. We tested the validity of the aerodynamic equation in a laboratory column experiment and found the equation to be valid for our experimental conditions. While the evaporation rate ranged from 15 to 38 mm d–1 in the test, the root mean square error between the measured and predicted evaporation rates was only 0.9 mm d–1. Although the resistance formulation is commonly used, some researchers have used a conductance variant of Eq. [1]. For example, Salhotra et al. (1985) analyzed evaporation rates from evaporation pans with saline water using the conductance form. Since a salt crust could provide additional resistance, rsc, to evaporation, we added this resistance to the denominator of Eq. [1] to give

Formula 2[2]
A similar resistance term was used previously to incorporate the effect of a dry surface layer (van de Griend and Owe, 1994) or of straw or gravel mulches (Farahani and DeCoursey, 2000) on evaporation.

Assuming thermodynamic equilibrium between the liquid and gaseous phases, the relative humidity in a soil or at the soil surface (hrs) can be calculated using the following equation (Philip and de Vries, 1957):

Formula 3[3]
where hre is the equilibrium relative humidity, {psi}w is the water potential (cm), Rv is the gas constant for water vapor (4697 cm K–1), and T is temperature (K). The water potential, {psi}w, is the sum of matric ({psi}m) and osmotic ({psi}o) potentials as follows:

Formula 4[4]
The osmotic potential, {psi}o (cm), of a solution can be estimated from (Campbell, 1985)

Formula 5[5]
where {varpi} is a unit-conversion factor (10.2 cm kg J–1), {nu} is the number of ions per molecule (i.e., two for NaCl), {chi} is the osmotic coefficient, C is the molar concentration of the solute (mol kg–1), and R is the universal gas constant (8.31 J mol–1K–1). The osmotic coefficient, {chi}, was calculated using the empirical equation

Formula 6[6]
where c is the concentration of the solute (mg cm–3), and ax, bx and px are fitting parameters. Values of these parameters as obtained by fitting data from the literature are listed in Table 1.


View this table:
[in this window]
[in a new window]
 
Table 1. Fitted parameter values of the dependence of the osmotic coefficient on the concentration of each solute.

 
Equation [4] indicates that both the water flow and solute transport equations must be solved to calculate {psi}w. The one-dimensional water balance equation for the combined liquid and gaseous phases is given by

Formula 7[7]
where {theta} is the volumetric water content (cm3 cm–3), t is time (s), ql is the liquid water flux (cm s–1), qv is the water vapor flux (cm s–1), z is depth (cm), and Sw is a sink term (assumed to be zero in this study). The liquid water flux, ql, is described using Darcy's law:

Formula 8[8]
where K is the hydraulic conductivity (cm s–1), {rho}w is the density of the soil solution (g cm–3), and {rho}w0 is the density of pure liquid water at the reference temperature (g cm–3). By using {rho}w/{rho}w0 instead of unity, the effect of increased density on the gravitational potential can be considered. We have incorporated the change of {rho}w with concentration. We neglected the contribution of the osmotic potential to liquid water flow since osmotic potential will be a driving force only during limited conditions where salt sieving in a dry clayey soil restricts solute diffusion (Letey et al., 1969; Hillel, 1998). Note that K does not contain a water vapor component. When the viscosity of a soil solution changes with temperature or the solute concentration, K may be expressed as

Formula 9[9]
where µ is the viscosity of the soil solution, while the subscripts 0 indicate reference values. In this study, we did consider the dependence of µ on solute concentration. The value of K may also depend on the exchangeable sodium percentage (ESP) and salinity (Hillel, 1998). Since we used a coarse-textured soil with little clay, we neglected the dependence of the hydraulic conductivity on ESP and salinity.

When the temperature gradient is negligible, the water vapor flux, qv, can be expressed as follows (Campbell, 1985; Mehta et al., 1994):

Formula 10[10]
where a is the air-filled porosity, {tau}g is the tortuosity for gas transport, Dva is the diffusion coefficient of water vapor in free air (g cm–2s–1), and {rho}w is the density of water (0.997 g cm–3 at 25°C). We used a value of 0.66 for the tortuosity (Penman, 1940; Cass et al., 1984, Daamen and Simmonds, 1996). The dependence of Dva on temperature was described as (Kimball et al., 1976)

Formula 11[11]

The saturated vapor density, {rho}v*, is also temperature dependent as follows (Kimball et al., 1976):

Formula 12[12]

Since our study focuses on the effect of salinity on evaporation and evaporation-driven solute transport, we conducted the experiments under isothermal conditions. For this reason, we did not analyze thermally induced water (vapor) flow or heat movement.

One-dimensional nonreactive solute movement in homogeneous soil is commonly analyzed with the CDE, which results from combining equations for the mass balance and the solute flux as follows:

Formula 13[13]
where qs is solute flux density (mg cm2 s–1), s is the mass of crystal per unit volume (mg cm–3), cmax is the saturated concentration (mg cm–3), and D is the dispersion coefficient (cm2 s–1), being the sum of the effective ionic diffusion coefficient, Di, and the mechanical dispersion coefficient, Dm:

Formula 14[14]

The effective ionic diffusion coefficient, Di, of a solute is the product of the temperature-dependent aqueous diffusion coefficient, Di0, and the tortuosity factor for ionic diffusion, {tau}s, the latter being a function of water content:

Formula 15[15]
The mechanical dispersion coefficient, Dm, is commonly assumed to be proportional to the pore water velocity, v (cm s–1) as follows:

Formula 16[16]
where {lambda} is dispersivity (cm). In unsaturated soils, Dm may also depend on the water content (Toride et al., 2003). Since our experiments were conducted at water contents similar to those used for determining the dispersivities, we neglected the dependency of Dm on water content.

Although some ion exchange may have occurred in our experiments, the retardation factor may be set equal to unity since we measured concentrations using the electrical conductivity. Released or exchanged ions would contribute nearly equally to the electrical conductivity as the original solutes.

Assuming that crystallization occurs instantaneously as the concentration reaches saturation (cmax), and that supersaturation does not persist for any significant period, the increase in the mass of crystal per unit volume, s (mg cm–3), and unit time is given by

Formula 17[17]


    MATERIALS AND METHODS
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 THEORY
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
In this study, we used Masa loamy sand (Udorthent) and Toyoura sand taken from the Tottori and Yamaguchi prefectures, respectively, in Japan. Results of particle size, bulk density, and saturated hydraulic conductivity measurements are presented in Table 2. Soils were all thoroughly leached and air dried before packing. All experiments were conducted at 25 ± 1°C. Saturated hydraulic conductivities were measured using the falling head method.


View this table:
[in this window]
[in a new window]
 
Table 2. Characteristics of the two soils used in this study.

 
Evaporation Experiment
Each soil was packed as uniformly as possible to predetermined bulk densities in soil columns of 3.8-cm diameter and 5.2-cm height. The walls of the columns consisted of acrylic rings, either 0.5 or 1.0 cm high, while a porous plate (5.0 mm thick) was used at the bottom. A thermocouple was inserted horizontally at a depth of 0.2 cm. Because slight shrinkage due to both desorption and evaporation was observed in several preliminary experiments, soil was added so that the soil surface were always about 1 mm higher than the top of the column wall. Figure 1 illustrates the experimental setup. The porous plate was hydraulically connected to a Mariotte flask to supply solution while maintaining a certain pressure head at the bottom. The hydraulic conductivity of the 0.5-cm-thick porous plate was 1.2 cm h–1. The soil samples were initially saturated from the bottom with a 3000 g m–3 NaCl or KCl solution. After reaching nearly saturation, the pressure head at the bottom, {psi}bot, was decreased to –75 cm (loamy sand) or –40 cm (sand) by lowering the level of the Mariotte flask. We conducted three runs with different combinations of soil and solute: Masa loamy sand and NaCl (A), Masa loamy sand and KCl (B) and Toyoura sand and NaCl (C).


Figure 1
View larger version (24K):
[in this window]
[in a new window]
 
Fig. 1. Experimental setup. The lamp to heat the surface was intermittently switched off automatically.

 
About 2 h after decreasing the pressure head at the bottom of the column, when outflow from the column had practically ceased, the soil surface was uncovered to allow evaporation to occur under nearly constant meteorological conditions, except for radiation, which was automatically regulated using a thermostat such that the soil temperature remained more or less constant at 25°C. Hereafter, time = 0 refers to the initiation of the evaporation phase. Evaporation was accelerated by blowing air across the soil surface with an electric fan. Experiment A was conducted in a relatively small laboratory in which the location of electric fans and columns could not be flexibly adjusted. The aerodynamic resistance, ra, of the soil columns for this reason was assumed constant for a given location, but different among the columns. Experiments B and C were performed in a wind tunnel. Values for ra and the relative humidity of the ambient air are listed in Table 3. The wind speed at a height of 2 cm for Exp. C was about 0.7 m s–1. The value of ra for each run or column was estimated from the evaporation rate during 1 h from a similar soil column nearly saturated with distilled water, at the same location and subject to the same wind speed as the salinized column. Rearranging Eq. [1] allows ra to be estimated from known values of E, hra, hrs, Ta, and Ts. The outside of each column was wrapped with 2.0-cm-thick polystyrene foam for heat insulation, while a doughnut-shaped brim was placed around the top of each column so that air would flow smoothly across the soil surface. The brims were made of white polyvinyl chloride, being 1 mm thick and 45 mm wide.


View this table:
[in this window]
[in a new window]
 
Table 3. Aerodynamic resistances and relative humidities for the evaporation experiments

 
In a study comparing evaporation rates from salinized and solute-free sand and glass bead columns containing water tables, Shimojima et al. (1996) reported 70 and 30% lower evaporation rates from the salinized sand and glass-bead columns, respectively, than from the solute-free columns. Nassar and Horton (1999) also compared evaporation rates from salinized and solute-free columns with impermeable bottom boundaries. They found that the salinized columns had 5 to 22% lower rates of evaporation than the solute-free columns. In such cases, the effect of salinity may be masked by a decrease in the evaporation rate due to a decrease in the matric potential during the desorptive phase. To avoid this uncertainty, we maintained a relatively large (–75 or –45 cm) {psi}bot to keep the soil surface wet enough so that the effect of the matric potential would be negligible, thus isolating the effect of salinity. The concentration of the supplied solution was always 3000 g m–3.

Water loss due to evaporation from the columns was manually measured with an electronic balance at 12-h intervals. We expected that the relative humidity just above the soil surface would not be uniform, since the water vapor evaporated from the upstream side is transported downstream. For this reason we turned the columns manually 180° after each weighing. Three replicates were prepared for each condition, with the soil columns being sectioned one by one to determine the water content and solute concentration profiles. Water contents were determined by oven drying at 105°C, while solute concentrations were estimated from the electrical conductivity of 1:5 soil/water extracts.

Hydraulic Properties
The soil water retention data for the soils were measured using the hanging water column and vapor equilibrium methods. Since the evaporation experiments involved very little wetting after the initial saturation phase, we do not present wetting retention data here.

For Masa loamy sand, the retention data were fitted with the empirical equation (Fujimaki and Inoue, 2003a)

Formula 18[18]
where {theta}sat is {theta} at {psi}m = 0, and {alpha}, {zeta}, n, and m are fitting parameters. The parameter {psi}0 is the pressure head where the water content becomes nearly zero (i.e., oven dry). In this study, {psi}0 was set to –107 cm, while m was handled as an independent fitting parameter. Fitted parameter values were: {theta}sat = 0.43, {alpha} = 0.042, {zeta} = 0.031, n = 1.37, and m = 0.32.

The retention data for Toyoura sand were fitted with the bimodal retention function:

Formula 19[19]
where w, {alpha}1, {alpha}2, n1, and n2 are the fitting parameters. We obtained the following parameter values: {theta}sat = 0.443, w = 0.748, {alpha}1 = 0.026, {alpha}2 = 0.021, n1 = 33.6, and n2 = 1.84. Data and fitted curves are shown in Fig. 2 .


Figure 2
View larger version (24K):
[in this window]
[in a new window]
 
Fig. 2. Soil water retention curves for the two soils used in this study.

 
The unsaturated hydraulic conductivity of Masa loamy sand was determined using a transient evaporation method (Fujimaki and Inoue, 2003a). For this soil, we assumed applicability of the empirical function (Campbell, 1974)

Formula 20[20]
where Ksat is the saturated hydraulic conductivity ({psi}m = 0) and {omega} is a fitting parameter whose value was estimated to be 7.2 using inverse parameter estimation.

For Toyoura sand, a flux-controlled steady-state evaporation method was used for determining Ksat (Fujimaki and Inoue, 2003b). Measured hydraulic conductivity data were fitted with the empirical equation:

Formula 21[21]
where ak and bk are fitting parameters. Analysis of the data produced ak = 1.57 and bk = 0.023. Figure 3 shows plots of the estimated hydraulic conductivity functions.


Figure 3
View larger version (21K):
[in this window]
[in a new window]
 
Fig. 3. Hydraulic conductivity for the two soils used in this study.

 
Solute Transport Properties
The dependency of the tortuosity factor for ionic diffusion, {tau}s, on the water content for each soil was determined using a transient-state method with two soil blocks (Kemper, 1986; Mehta et al., 1995; Olesen et al., 2000). A finite difference numerical solution of the diffusion equation (CDE with v = 0) was used for the inverse parameter estimation approach. The objective function, O(D), to be minimized for this purpose was assumed to be of the form

Formula 22[22]
where i is the sample number, N is the number of samples, and ccal is the calculated concentration. The golden section method (Press et al., 1989) was used to minimize the objective function.

The inversely estimated values of {tau}s for each water content are plotted in Fig. 4 . The S-shaped distribution for Toyoura sand is consistent with results by Mehta et al. (1995), who fitted an S-shaped distribution involving two power functions to their data. In our study, we fitted the data with the single equation

Formula 23[23]
containing four unknown parameters (as, bs, cs, and ds). For Masa loamy sand, ds was fixed at zero. Fitted parameter values are listed in Table 4.


Figure 4
View larger version (20K):
[in this window]
[in a new window]
 
Fig. 4. Tortuosity factor for ionic diffusion as a function of water content for the two soils used in this study.

 

View this table:
[in this window]
[in a new window]
 
Table 4. Fitted parameter values in Eq. [22] for each soil.

 
The mechanical dispersion coefficients were also determined by numerical inversion of the measured concentration profiles at fixed times. Air-dried Masa loamy sand for this purpose was packed into a 20-cm-long column having an inside diameter of 2 cm. A relatively small inner diameter was used to minimize three-dimensional flow near the soil surface induced by the inevitably nonuniform application to the soil surface. The bottom boundary was effectively that of a seepage face, consisting of a metal mesh along which air was blown using an electric fan to promote evaporation of the out-flowing water. We applied a constant water flux at the soil surface such that the resulting pore water velocity and water content were comparable to those of the evaporation experiment. When the wetting front reached a depth of about 15 cm, the solute concentration of the applied water was stepwise increased from zero to 5000 g m–3. Figure 5 shows that relatively uniform water content profiles were present in the columns once the wetting front reached the bottom. When the estimated solute front had reached a depth of about 10 cm, the column was sectioned to measure the concentration profile.


Figure 5
View larger version (19K):
[in this window]
[in a new window]
 
Fig. 5. Concentration and water content profiles for the displacement experiment for (a) Masa loamy sand and (b) Toyoura sand.

 
Measured and fitted concentration profiles for Masa loamy sand are shown in Fig. 5a. The slight increase in concentration near the bottom may have been due to solute initially contained in the air-dry soil sample, although the soil was well leached before air drying. For this reason, we used only data from the top 16 cm for the inverse analysis. In the numerical inversion, the inter- and extrapolated water content profile was used assuming that water flow was in a steady state.

Shiozawa and Fujimaki (2004) previously showed that infiltration into air-dried sand may lead to unstable flow, in which case fingering could occur if the cross-sectional area is relatively large. For this reason, we initially saturated the column from the bottom such that the bottom pressure head, {psi}bot, can be controlled by means of a ceramic plate and associated tubing. After saturation, {psi}bot was decreased to –40 cm by lowering the outflow point of the tube, while the inflowing water flux at the surface was maintained at the prescribed rate. When outflow became constant, the concentration was stepwise increased to 5000 g m–3. Figure 5b shows that the water content somewhat increased with depth, especially near the bottom. This nonuniform distribution probably was caused by having an insufficient {psi}bot to attain a uniform profile combined with the relatively low hydraulic conductivity of the ceramic plate. Since the water content varied little vs. depth except near the bottom, we still could assume a constant dispersivity. Our numerical solutions of the CDE using dispersivities optimized so that the objective function given by Eq. [22] became minimal agreed well with the measured profiles (Fig. 5).

Table 5 lists the inversely estimated dispersivities obtained with the above procedure (which considered also Di in Eq. [14]). We also fitted concentration profiles using an analytical solution assuming uniform water content distributions. An analytical solution of the CDE for a uniform initial value, constant flux at the upper boundary, and zero concentration gradient at infinite depth was derived by Lindstrom et al. (1967) and Parker and van Genuchten (1984). Fitting the analytical solution to the data gave similar dispersivities as the numerical inversion, as shown by the results in Table 5. We used the numerically determined dispersivities in the following simulations.


View this table:
[in this window]
[in a new window]
 
Table 5. Optimized dispersivities for each soil.

 
Numerical Simulation
The water flow equation (Eq.[7]) including isothermal vapor movement was solved using a finite difference method based on the mass-conservative iterative scheme proposed by Celia and Bouloutas (1990). Space increments, {Delta}z, were set at a constant value of 0.05 cm. Time steps, {Delta}t, were adjusted automatically such that, at each time step, the number of iterations was approximately five and the maximum change in ln({psi}m) was <0.693 [=ln(2)]. The initial pressure heads corresponded to an equilibrium pressure profile for a bottom boundary condition of {psi}bot = –75 cm (Masa) or –40 cm (Toyoura). The lower boundary condition was set at a constant pressure head of {psi}bot, while an atmospheric condition was used at the soil surface, with the evaporation rates at each time step being calculated using Eq. [1] or [2]. The water potential of the top element was used to calculate the surface relative humidity, hrs, using Eq. [3].

The CDE was solved with an implicit finite difference scheme using the same spatial and time discretization as for the water flow analysis. Numerical dispersion, which results in an overestimation of the spreading of the solute due to discretization of the convection term, was corrected using a scheme proposed by Bresler (1973) in which the numerical dispersion coefficient, Dn, is subtracted from D at each time step and for each element as follows:

Formula 24[24]

Formula 25[25]
where the subscript j denotes the discretized time level, and Dc is the corrected dispersion coefficient (which is actually multiplied by the concentration gradient in the numerical scheme). The initial condition was a uniform concentration of 3000 g m–3. The lower boundary condition was set at a constant concentration of 3000 g m–3, while the upper condition was that of a zero flux.

Computations of water and solute movement were performed sequentially at each time step.


    RESULTS AND DISCUSSION
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 THEORY
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Evaporation Experiment
Observed evaporation rates are shown in Fig. 6 . For each run, the evaporation rate initially decreased rapidly, followed by a more gradual decrease. Since the soil surfaces were wet enough to maintain the matric potential, {psi}m, at the soil surface within a range not affecting the relative humidity, hr, any reduction in the evaporation rate may be interpreted as being caused by salinity.


Figure 6
View larger version (19K):
[in this window]
[in a new window]
 
Fig. 6. Temporal change in the evaporation rate for (a) Masa loamy sand and NaCl (b) Masa loamy sand and KCl, and (c) Toyoura sand and NaCl. In (a), numerical solutions are for Column 3; rsc is the resistivity of the salt crust; half and quarter dispersivity is in reference to the dispersivity in the top 2 cm compared with the independently measured value.

 
Figure 7 shows a comparison of measured and simulated water content profiles. Different settings in the numerical simulation described below gave almost identical results for the water content profiles. The numerical solutions were also mostly invariant with time. Thus, only the solutions of the reference case (with the measured dispersivity and rsc > 0) for two examples are shown in Fig. 7. Matric potentials at the soil surface calculated from measured water content and retention curves were always larger than –100 cm, at which value a reduction of the evaporation rate does not occur under solute-free conditions. The numerical solutions were in good agreement with the measured profiles.


Figure 7
View larger version (21K):
[in this window]
[in a new window]
 
Fig. 7. Comparisons of measured and simulated water content profiles.

 
Measured concentration profiles are shown in Fig. 8 . Since the topmost sample contained crystals, the saturated concentration is plotted at the soil surface. For every run, solute mostly accumulated in the top 0.5 cm. A slight increase in concentration with time after about 48 h was observed for each run, but the increase was restricted to the upper 1.0 cm.


Figure 8
View larger version (14K):
[in this window]
[in a new window]
 
Fig. 8. Concentration profiles for (a) Masa loamy sand and NaCl, (b) Masa loamy sand and KCl, and (c) Toyoura sand and NaCl.

 
Simulated evaporation rates are plotted in Fig. 6. Dotted curves are the numerical solutions using Eq. [1] (i.e., with rsc neglected); these curve substantially overestimated measured evaporation rates. Since osmotic potential does not decrease further after reaching the saturation concentration of the solute, the predicted evaporation rates become constant. Any further decrease can therefore be attributed to the resistance to evaporation due to the presence of a salt crust.

By rearranging Eq. [2], rsc can be calculated from the evaporation rate, the aerodynamic resistance, ra, and the equilibrium relative humidity at the saturation concentration assuming {psi}m is negligible. Evaporation rates (Fig. 6) at the time when the columns were dismantled were used in the calculation. Calculated salt crust resistances are plotted in Fig. 9 as a function of accumulated salt, with each column representing one point in the figure. Because of the difficulty in extracting crystallized salt from the soil–solution–crystal mixture, we related rsc to the mass of salt above a depth of 0.25 cm, {Gamma} (mg cm–2):

Formula 26[26]


Figure 9
View larger version (28K):
[in this window]
[in a new window]
 
Fig. 9. Salt crust resistance as a function of accumulated salt in depth z <0.25 cm.

 
The maximum height of the salt crust was observed to be about 0.3 cm. The limit of –1 cm is a rather arbitrary value that makes allowances for the height. To incorporate the effect of rsc into the models, some relationship with {Gamma} must be established. We fitted the available data with a logarithmic function:

Formula 27[27]
where ar and br are fitting parameters. Having a threshold {Gamma} is reasonable, since a salt crust should form only after a certain amount of solute has accumulated near the soil surface. Figure 9 shows good agreement between the measured data and fitted relationship. Using an empirical function such as Eq. [27] also allows a more a systematic discussion of the differences in the distributions. The rsc({theta}) function for the KCl crust was consistently larger than that for the NaCl crust. Differences in the geometry of the crystals, which were visibly distinguishable, probably caused the different rsc({theta}) curves. The KCl crust appeared to be flatter and with fewer openings than the NaCl crust. On the other hand, no significant differences were apparent between the two soils for the same solute (NaCl). This lesser dependency on soil type seems reasonable since the salt crust is formed mainly at or above the soil surface. Using a magnifying glass, we observed the surface (z = 0.25 cm) that appeared just after slicing the top ring and found no crystals. Tanner and Shen (1990) investigated the dependency of the resistance to vapor transfer of flat surface residues on wind velocity. It is probable that rsc will also depend on wind velocity since the salt crust seems thin enough to allow air to flow through the crust. Further experiments are needed to evaluate the dependence of rsc on wind velocity. In our numerical simulation, we kept rsc at zero until crystallization started. When the calculated concentration of the top element reached saturation, rsc often already had a large value, leading to an unstable numerical solution of the evaporation rate. To stabilize the numerical analysis, {Gamma} was approximated using

Formula 28[28]
This approximation gradually increased rsc from zero when crystallization first started.

Evaporation rates simulated using the obtained rsc({theta}) relationship showed better agreement with observations in the later stages of the experiments, as shown in Fig. 6. Still, the evaporation rate was consistently overestimated before reaching the threshold value of {Gamma}. These results indicate that the numerical solution underestimated solute concentrations at the soil surface. In addition to the salt crust resistance itself, precision in the numerical solution of the surface solute concentration is critical for accurately predicting evaporation rates from saline soils.

Figure 10 shows numerical solutions of the concentration profiles for rsc > 0. The CDE in each case overestimated the downward movement of solute against the upward movement of water. This may be caused by the fact that the CDE uses the Fickian law of diffusion for describing mechanical dispersion. The only mechanism for backward (downward in this case) movement of solute is ionic diffusion since solute does not move against flow as a result of mechanical dispersion. A serious limitation of the CDE is that this equation macroscopically cannot distinguish between diffusion and mechanical dispersion processes, and hence by necessity these two different processes are lumped together into one single term. A fictitious example of how CDE can predict physically unrealistic backward mixing is shown in Fig. 11 . Assume that water is flowing at a uniform and constant rate of 1.0 cm h–1 to the right in an infinite medium whose dispersivity is 0.5 cm. If a Dirac delta is initially present at x = 5 cm and the diffusion coefficient is negligible compared with mechanical dispersion, then the concentration should not increase upstream from the initial condition. The solution of the CDE, however, in this case equivalent to a normal distribution, predicts the physically odd situation as portrayed in the figure. Toride et al. (1993) similarly pointed out that the CDE can predict physically inappropriate distributions just after starting from a stepwise initial distribution.


Figure 10
View larger version (17K):
[in this window]
[in a new window]
 
Fig. 10. Comparison of measured and simulated concentration profiles for (a) Masa loamy sand and NaCl, (b) Masa loamy sand and KCl, and (c)Toyoura sand and NaCl; rsc is the resistivity of the salt crust; half and quarter dispersivity is in reference to the dispersivity in the top 2 cm compared with the independently measured value.

 

Figure 11
View larger version (17K):
[in this window]
[in a new window]
 
Fig. 11. Solution of the convection–dispersion equation assuming solute-free input to a soil having a Dirac delta initial distribution at x = 5 cm; {lambda} = dispersivity, Di = effective ionic diffusion coefficient, v = pore water velocity.

 
The mechanical dispersion and diffusion coefficients at a depth of 0.05 cm, and the times displayed in Fig. 10, are listed in Table 6. In soils where water flows at a rate of approximately 10 mm d–1, the mechanical dispersion coefficient is dominant or of a comparable magnitude to the diffusion coefficient. This shows that Dm is crucial for accurate transport predictions.


View this table:
[in this window]
[in a new window]
 
Table 6. Comparison of diffusion and dispersion coefficients at depth z = 0.05 cm.

 
The osmotic potential gradient may induce some vaporization and vapor diffusion near the soil surface. The osmotically driven vapor flux can be calculated using Eq. [3]GoGo to [6] and [10]Go to [12]. If internal vaporization, which occurs when a vapor concentration gradient exists, is overestimated, concentrations above the source of the water vapor would also be overestimated. A sensitivity analysis with qv = 0, however, indicated that for our experimental conditions, internal vapor movement had little effect on the numerical solution of the concentration profiles, since the vapor flux was small (<2.0 mm d–1) compared with the evaporation rate, and vaporization occurred only in an extremely shallow layer (z < 0.2 cm).

When we reduced the dispersivity by a factor of two or four (i.e., to {lambda}/2 or {lambda}/4) for transport in the top 2 cm of the columns, significantly better predictions were obtained for both the evaporation rates and the concentration profiles, as shown in Fig. 6 and 10, respectively. Using the reduced dispersivity values (even if merely an initial guess), the model predicted earlier crystallization, which corresponded to the actual visual observation. Assuming Dm = 0 at or near the soil surface may be inappropriate since for relatively broader variations in the pore water velocity, ions may have more opportunity to move backward. This does not mean, however, that the same dispersivity as in the deeper soil should be used near the evaporating soil surface. Fujimaki et al. (1997) reported a similar discrepancy for transport in Tottori dune sand, and found good results when the dispersivity was decreased by a factor of 3 compared with the independently measured value. Our results now confirm overestimation of backward diffusion using the CDE for two other soils.

Even when the dispersivity was reduced by a factor of four, the model still overestimated the evaporation rate during the early stage of evaporation for the Masa soil. Since the soil surface was wet, the sharp decrease in the evaporation rate cannot be explained by a reduction in the matric potential. Possible reasons may be (i) underestimation of the solute concentration at the surface due to overestimation of the tortuosity factor for ionic diffusion, and (ii) nonequilibrium in the relative humidity (a violation of Eq. [3]) under conditions of high salinity. Other unknown mechanisms may also exist, however. Further studies are needed to more accurately predict evaporation rates from wet saline soils.


    CONCLUSIONS
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 THEORY
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
We conducted laboratory column transport experiments under constant meteorological conditions, except for radiation, which was automatically regulated so that the soil temperature remained constant and the same as that of the air. The concentrations of the initial soil solution and the water flowing into the bottom were set at 3000 g m–3, while the height of the column was 5.2 cm. Evaporation experiments were performed with three combinations of soil and solute: sand and NaCl, loamy sand and NaCl, and loamy sand and KCl. Even though the soil surface was kept wet by keeping the suction at the bottom low (75 or 45 cm), the evaporation rate was found to considerably decrease with time. This decrease could not be explained by osmotic potential alone. The bulk transfer equation for evaporation was therefore modified to include an additional resistance to water vapor diffusion caused by the formation of a salt crust. The dependence of the salt crust resistance on the mass of accumulated salt was evaluated. Using the salt crust resistance, we calculated solute transport and the evaporation rate. Water vapor movement and the effect of osmotic potential were also considered. Results using independently obtained hydraulic and solute transport parameters showed that the CDE overestimated backward diffusion near the evaporating soil surface, which resulted in a significant delay in salt accumulation at the soil surface, as well as a delay in the evaporation rate decrease. Since the CDE uses an analogy of Fickian law to describe mechanical dispersion, the dispersion term overestimated downward movement of solute against the upward convective transport process. When the dispersivities for transport in the top 2 cm were reduced by a factor of two compared with the independently measured values, substantially better predictions were obtained for both the evaporation rates and concentration profiles.


    ACKNOWLEDGMENTS
 
We thank the National Research Institute for Earth Science and Disaster Prevention, Japan, for the use of their wind tunnel.


    REFERENCES
 TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 THEORY
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 





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