Study of the effect of wind speed on evaporation from soil through integrated modeling of the atmospheric boundary layer and shallow subsurface.

In an effort to develop methods based on integrating the subsurface to the atmospheric boundary layer to estimate evaporation, we developed a model based on the coupling of Navier-Stokes free flow and Darcy flow in porous medium. The model was tested using experimental data to study the effect of wind speed on evaporation. The model consists of the coupled equations of mass conservation for two-phase flow in porous medium with single-phase flow in the free-flow domain under nonisothermal, nonequilibrium phase change conditions. In this model, the evaporation rate and soil surface temperature and relative humidity at the interface come directly from the integrated model output. To experimentally validate numerical results, we developed a unique test system consisting of a wind tunnel interfaced with a soil tank instrumented with a network of sensors to measure soil-water variables. Results demonstrated that, by using this coupling approach, it is possible to predict the different stages of the drying process with good accuracy. Increasing the wind speed increases the first stage evaporation rate and decreases the transition time between two evaporative stages (soil water flow to vapor diffusion controlled) at low velocity values; then, at high wind speeds the evaporation rate becomes less dependent on the wind speed. On the contrary, the impact of wind speed on second stage evaporation (diffusion-dominant stage) is not significant. We found that the thermal and solute dispersion in free-flow systems has a significant influence on drying processes from porous media and should be taken into account.

1. Introduction

A critical component of the water cycle at local, regional, and global scales is evaporation. Many researchers focus on evapotranspiration (ET), to include both bare soil evaporation and the vegetative portion (transpiration). Nonetheless, what is fundamental to the process of mass transfer across soil/atmospheric interfaces happens in bare soil; vegetation is a medium that cuts across the soil (roots) and atmosphere (air). Bare soil evaporation in arid or semiarid settings can account for more than half of ET and is therefore critical to its understanding [Huxman et al., 2005] and has been found to be very difficult to project in modeling studies [Seager et al., 2007]. Even though decades of research have improved our understanding of bare soil evaporation, many knowledge gaps still exist in the current science on how the soil water in the shallow subsurface close to the land surface interacts with the air in the atmosphere. Understanding this interaction is paramount to our understanding of many emerging problems to include climate change, water and food supply, leaking of geologically sequestered CO2 from soil [Oldenburg and Unger, 2004], the accurate detection of buried objects such as landmines [Das et al., 2001], and the remediation of contaminated soil in the shallow subsurface [Weaver and Tillman, 2005].

The rate of soil evaporation is affected by atmospheric conditions (e.g., humidity, temperature, thermal radiation, and wind velocity), and thermal, and hydraulic properties of soil (thermal and hydraulic conductivity, porosity), all of which are strongly coupled. This strong coupling between processes leads to highly dynamic interactions between the atmosphere and soil resulting in dynamic evaporative behaviors [Sakai et al., 2011]. However, the atmospheric coupling to the soil at the land-atmospheric interface is rarely considered in most current models or practical application. This is due to the complexity of the problem in field settings and the scarcity of field or laboratory data capable of testing and refining energy and mass transfer theories. In most efforts to compute evaporation from soil, only indirect coupling is provided to characterize the interaction between multiphase flow in soil under realistic atmospheric conditions even though heat and mass flux are controlled by the coupled dynamics of the land and the atmospheric boundary layer. It is recognized that the most important process that determines the coupling between the soil water and heat is the transport of latent heat (the result of phase change) by vapor flux in the unsaturated soil pores and at the interface between the soil and the atmosphere [Bittelli et al., 2008]. Models that incorporate these processes have been developed [e.g., Jassal et al., 2003]; however, as Bittelli et al. [2008] note, a detailed experimental verification of vapor movement above the soil surface (i.e., atmospheric boundary layer) has not been conducted. Bittelli et al. [2008] suggest that the errors introduced in the vapor flow calculations are due to a number of factors that include lack of proper coupling of the thermal and mass flux processes, deficiencies in the constitutive relationships (e.g., thermal and hydraulic conductivities and soil water content) and difficulty in determining the resistance parameters at the land (soil)-atmospheric interface. For example, a prevalent modeling approach is to derive the aerodynamic and soil surface resistance terms based on semiempirical or empirical approaches and to adjust the predicted evaporation based on true conditions that depend on ambient conditions such as soil moisture, roughness, wind speed, etc. Traditionally, the influences of atmospheric conditions are applied at the soil surface and aerodynamic resistance is applied on the border between the air flow and permeable media [e.g., Bittelli et al., 2008; Camillo and Gurney, 1986; Novak, 2010; van de Griend and Owe, 1994]. In these cases, evaporation rate (E) can be determined by:

where rs is the soil surface resistance for water vapor transport, rv is the aerodynamic resistance for water vapor, (ρ) is the vapor density immediately below the soil surface (in porous medium), and (ρ) is the vapor density immediately above the soil surface (in free medium). The vapor density above the soil surface is calculated based on the measurement of relative humidity on the surface of the porous medium in the free-flow medium.

In equation (1), the aerodynamic resistance for vapor transport depends on surface roughness properties and wind speed [Bittelli et al., 2008; Campbell, 1977]. The soil surface resistance depends on soil surface water content. The relationship between aerodynamic resistance to vapor transport and soil water content is typically an exponential form; there are many exponential empirical functions used to describe this relationship [Camillo and Gurney, 1986; van de Griend and Owe, 1994]. Although this approach is widely used, modeling comparison studies have shown significant variation between model parameterizations and evaporative fluxes [Desborough et al., 1996; Schmid, 1997; Smits et al., 2012; Villagarcia et al., 2007]. Recently, with the goal of addressing the issue of coupling the land to the atmosphere, Smits et al. [2012] evaluated three different modeling approaches of bare soil evaporation formulated with different land surface boundary conditions and compared modeling results to laboratory generated experimental data. Results demonstrated that no one approach could be deemed most appropriate for every situation, demonstrating that further work focusing on the land/atmospheric interface, properly incorporating the complex interactions between the land and the atmospheric boundary layer is needed to increase the understanding of the processes that control shallow subsurface soil moisture flow that controls bare soil evaporation.

The modeling of nonisothermal single-phase (two-component) transfer in the atmosphere and two-phase (two-component) transfer in porous media have been separately investigated by many authors [e.g., Chao-Yang and Beckermann, 1993; Niessner and Hassanizadeh, 2009; Wang and Cheng, 1997]. Recently, numerical advances have been made in the coupling of free flow (Navier-Stokes) with porous media flow (Darcy flow) [Baber et al., 2012; Chidyagwai and Riviere, 2011; Mosthaf et al., 2011; Nield, 2009; Shavit, 2009], however, these models were not adequately validated with experimental data. Mosthaf et al. [2011] extended the classical single-phase coupling to two-phase flow in porous media and one phase in the free flow. Their model is based on the continuity of fluxes at the porous medium-free medium interface and use of the Beavers-Joseph boundary condition [Mosthaf et al., 2011]. Baber et al. [2012] focused on the numerical concept and its implementation into a local modeling toolbox. The numerical parametric study showed that the proposed model can predict the evaporation phenomenology correctly. They concluded that the variation of permeability influences the duration of the capillary-driven evaporation regime whereas the variation of temperature affects the magnitude of the evaporation rate. They also showed that the choice of the Beavers-Joseph coefficient has a negligible influence on the evaporation rate across the interface [Baber et al., 2012]. However, the aforementioned models did not investigate nonequilibrium multiphase flow under nonisothermal conditions.

The equilibrium assumption (i.e., instantaneous phase change) is called into question in natural evaporation from soil, in which atmospheric conditions induce mass transfer at the surface and tend to shift to equilibrium [e.g., Chammari et al., 2008]. However, the common approach when modeling the movement of liquid water, water vapor, and heat in the soil immediately below the land-atmosphere interface is to assume that water vapor in air is always in equilibrium with liquid water in the pores [e.g., Prat, 2002; Sakai et al., 2011; Shokri et al., 2009; Whitaker, 1977]. In an equilibrium formulation hypothesis, a sufficiently rapid mass exchange between the liquid water and its vapor is implied and it is assumed that the vapor pressure is always equal to the saturated vapor pressure. However, experimental studies [Armstrong et al., 1994; Chammari et al., ,; Lozano et al., 2009; Lozano et al., 2008] have observed a volatilization or condensation time. A recent comparison between equilibrium and nonequilibrium phase change modeling concepts using experimental data by Smits et al. [2011] revealed that the nonequilibrium model is able to better predict evaporation in porous media under certain conditions. Another limitation of the previously mentioned modeling efforts is that they do not take into account dispersion in the free-flow domain. However, free-flow dispersion, which depends on the Péclet number, may be the dominant regime in the free-flow region.

It is well known that a no-slip condition at the free-flow and porous domains surface is not satisfactory and indeed a slip boundary condition occurs. The slip boundary condition was first obtained experimentally by Beavers and Joseph [1967]. Beavers and Joseph proposed that the tangential component of the normal stress of the flow at the free-flow and porous medium interface is proportional to the jump of the tangential velocity across the interface [Beavers and Joseph, 1967]. The coupling condition was further studied by Saffman [1971] who concluded that the filtration velocity in porous media was much smaller than the free-flow velocity and can be neglected [Saffman, 1971]. There exist several other formulations for a slip boundary condition to include (1) using a shear stress jump condition by means of the nonlocal form of the volume averaging technique with an experimentally determined fitting parameter [Ochoa-Tapia and Whitaker, 1997], or (2) using the inertia and boundary effects [Vafai and Kim, 1990]. Alazmi and Vafai [2001] compared five different interface conditions between the porous medium and adjacent fluid layer. They concluded that the velocity field is more sensitive to variation in the boundary condition than the temperature field [Alazmi and Vafai, 2001]. They showed similar results for all five interface conditions. Therefore, the Beaver's Joseph formulation was chosen for this work based on simplicity of the implementation to the numerical model.

In this study, in order to better understand the coupling between free flow and porous media flow, we developed a theory for coupling single-phase (gas), two-component (air and water vapor) transfer in the atmosphere and two-phase (gas, liquid), two-component (air and water vapor) flow in porous media at the REV scale under nonisothermal, nonequilibrium conditions and taking into account the dispersion in free-flow medium (Figure 1). Mean wind velocity in the free fluid system, the soil thermo-physical properties and the initial conditions in the soil and atmospheric systems are the only input parameters needed for this model. In order to test the numerical formulations and codes, we performed a series of laboratory experiments under varying wind speeds using bench scale physical models and a unique low velocity porous media/wind tunnel, allowing for better control and gathering of accurate data at scales of interest not feasible in the field.

Figure 1

Schematic of the problem configuration for a single phase free fluid that interacts with two fluid phases in porous media.

In addition to bare soil evaporation, the theory developed in this work can also be used for various applications, such as environmental process involving vapor transport across interfaces of porous media and free flow (evaporation and condensation from vadose zone with a surface flows), fuel cells (gas diffusion layers of fuel cells are typically made of fibrous materials is in contact with free gas distributor), biological systems (interaction between blood flow and surrounding tissues, the dynamics of growing biofilm), and development of food drying process (interaction between air conditions and hygroscopic porous food).

2. Numerical Model Formulation

In this section, we present the formulations of macroscopic flow and, mass and energy balance equations for water, air and water vapor in the porous media as well as the gas (air and water vapor) above the porous media (Figure 2). The free flow and porous medium flow and transport equations are given in sections 2.1 and 2.2, respectively. The different boundary conditions applied on the interface between porous medium and free medium are listed in section 2.3.

Figure 2

Two dimensional problem configuration (subdomains, boundary and initial conditions) where T is temperature, w is water vapor concentration, U is wind speed, Sw is the water phase saturation, Jc and JT are mass and heat fluxes, Pc is capillary pressure, and K is intrinsic permeability.

2.1. Governing Free-Flow Medium Transport Equations

We assume nonisothermal, single-phase gas flow in the free-flow domain. The gas phase, identified here with the subscript “g”, is composed of two components, water vapor and air. The water vapor component is identified with the subscript “ν.” We consider in this domain that the density and viscosity of the moist air depends on the temperature and mass fraction of water vapor. The Navier-Stokes equation describes the fluid flow in the free domain. The continuity and Navier-Stokes equations assuming incompressible fluid, with no thermal and solutal expansion, are given, respectively, by [Bird et al., 2002]

The duct flow Reynolds numbers calculated in the range of the wind experimental velocities applied in this study show a laminar flow and transition flow regime. However, the boundary layer Reynolds numbers stay less than the critical value. In this study, we use Stokes equation for the flow motion in free-flow region, assuming a laminar, Poiseuille flow between two parallel plates and consequently, negligible vapor transfer by turbulent diffusion is valid. Nevertheless, the surface turbulence can significantly affect the transport just beneath the soil surface [Ishihara et al., 1992] and should be addressed in the future.

The energy transfer in the free-flow region is described by the following equation for the gas phase [Kaviany, 2001]

where T is the moist air temperature, and is the thermal dispersion tensor of moist air.

The component mass conservation is described by the following equation for the free fluid domain [Bird et al., 2002]

where w is the mass fraction of water vapor component in the gas-phase in free-flow domain and D is the Taylor dispersion tensor of water vapor in air.

We define the thermal and solutal Péclet numbers as

respectively, where H is the height of free-flow region, λ is the thermal conductivity of air-vapor binary mixture (or moist air thermal conductivity), u is the velocity field in x-direction and is the binary diffusion coefficient of water vapor in dry air which is a function of temperature [Campbell, 1985]. The longitudinal Taylor dispersion coefficients of equations (4) and (5) in a stratified media can be calculated as [Wooding, 1960]

2.2. Governing Porous Medium Transport Equations

Natural soil in the unsaturated zone is often represented as a three-phase system, consisting of the solid, liquid, and gas phases. In this study, we consider that two components are present in the gas-phase: dry air and water vapor. The solid phase is assumed to be inert and nondeformable. The soil pores are occupied by water as a liquid phase, and a gaseous binary mixture consisting of dry air and water vapor. We consider that the amount of air dissolved in liquid water is negligible; therefore, the liquid phase can be treated as a single component. In this section, we use the same notation as in the free medium. When necessary, the superscripts (ff) and (pm) are used to denote the free-flow and porous medium subdomains, respectively.

2.2.1. Two-Phase Flow and Transport

The flow of two noncompressible immiscible fluids in a rigid porous medium can be described by two coupled pressure-pressure partial differential equations [Bear, 1972]

where the water and gas velocity can be defined on the basis of the extended Darcy's law for multiphase flow, with definition of relative permeability as

Here the subscript ℓ and g denote the liquid phase (wetting phase) and gas phase (nonwetting phase), respectively. v and v are the velocities within the wetting and nonwetting fluids, S denotes fluid saturation, p pressure, ρ density, K intrinsic permeability tensor and k relative permeability of the phases, μ fluid dynamic viscosity, g gravitational acceleration, ϕ total porosity of the porous medium, and phase change rate between water and its vapor due to evaporation or condensation. The nonequilibrium effect is calculated and implemented to the model through this phase change rate. We will show in section 2.2.3 how it can be calculated from fluid and porous medium properties. The comparison between the equilibrium and nonequilibrium model for a porous medium-free flow coupled system will be investigated in a future study. Nonetheless, the reader can see Smits et al. [2011] for a comparison between equilibrium and nonequilibrium formulations without coupling between free flow and porous medium flow.

Equations (9) and (10) can be solved simultaneously for the main unknown pressures p and p. Since the medium is incompressible and the two fluids jointly fill the void space,

The normalized wetting saturation (or effective wetting saturation) is defined as

where S is the residual saturation related to pore scale trapping. The residual saturation changes linearly as a function of temperature as shown in She and Sleep [1998] and was implemented in this work. The difference between the nonwetting and wetting pressure is known as the capillary pressure and can defined as a function of the normalized wetting fluid saturation

To describe the relations between p, k, k, and S the analytical models of van Genuchten-Mualem can be used as [Mualem, 1976; van Genuchten, 1980]

where α and n are the Van Genuchten soil parameter . Here because p is large and the changes in other parameters are small, the equivalent height of water or capillary pressure head is used instead of capillary pressure as . The relative permeability of the wetting phase are specified according to the van Genuchten–Mualem model [Mualem, 1976; van Genuchten, 1980]

The van Genuchten–Mualem model for relative nonwetting phase permeability can be written as [Parker et al., 1987]

The mass balance for water vapor in the gas phase can be expressed as [Bear, 1972]

Because of slow flow velocities in porous media and high diffusion coefficients, dispersion can be neglected in the gas phase. Therefore, the effective vapor diffusion in porous medium can be predicted as

The tortuosity τ can be estimated as [Millington and Quirk, 1961]

In order to take into account the water vapor flux enhancement by thermal gradients, the empirical vapor enhancement factor from [Cass et al., 1984] was multiplied in this study by the vapor diffusion coefficient.

2.2.2. Energy Balance Equation

When the principle of local-scale thermal equilibrium is valid, a one-equation equilibrium model, which consists of a single transfer equation can be written as .

The characteristic time associated with thermal equilibrium is much lower than the characteristic time associated with mass transfer. Therefore, using the local thermal equilibrium assumption is acceptable for this work.

The macroscopic form of the governing equation for T can be then written as [Whitaker, 1977]

where is the effective thermal conductivity of the combined three phases and depends on structure, porosity, the ratio of the thermal properties of the solid phase on the fluid phase, and dispersion in porous media. Q is the heat loss term and L is the latent heat of water vaporization which is function of temperature [Monteith and Unsworth, 1990]. The model proposed by Campbell et al. [1994] was used to determine the effective soil thermal conductivity [Campbell et al., 1994]. It is based on the assumption that thermal conductivity of any mixture can be expressed as the weighted sum of the thermal conductivities of the individual components of the mixture. The effective thermal conductivity is given by

where λ, λ, λ are the thermal conductivities of water, gas, and soil matrix and ω, ω, and ω are the respective weighting factors for each phase calculated according to [Campbell et al., 1994].

2.2.3. Nonequilibrium Phase Change

As mentioned in the introduction, the assumption of equilibrium mass exchange is called into question for many natural drying applications. Bénet et al. [2009] evaluated several experimental studies to understand the thermodynamics of nonequilibrium phase change during soil drying. In their work, a nonequilibrium situation is created experimentally by extracting the gas phase of a small soil sample, and replacing it with dry air under isothermal conditions. Their proposed model involves three coefficients that must be determined experimentally. The volumetric phase change, , can be calculated based on the framework of linear thermodynamics of irreversible processes where fluxes are expressed as a linear function of forces

where P is the partial pressure of vapor, and P is the partial pressure of vapor in equilibrium [Bénet et al., 2009], M is the molar mass of water, and L is a phenomenological coefficient that depends on water saturation and the temperature. The phenomenological coefficient is characterized by three coefficients through the following expressions:

where L, k, and r are determined experimentally. Unfortunately, the experimental data for these three coefficients are limited to specific soil, temperatures and total gas pressures, and therefore, cannot be used in this study. In addition, hygroscopic effects have a significant influence by lowering the phase change velocity [Bénet et al., 2009]. In the experimental part of this study, we have used sandy soil; therefore, phase change will be faster compared to the soils used in [Bénet et al., 2009].

An alternative method of representing the term can also be used [Bixler, 1985; Le et al., 1995; Scarpa and Milano, 2002; Zhang and Datta, 2004]. In this study, we assumed that the phase change rate is proportional to the difference between the equilibrium density of vapor and its actual density through an appropriate time delay coefficient t (relaxation time to obtain the equilibrium within the pore space). This method also takes into account the water availability for evaporation by including soil water content.

where ρ = ρ is the vapor density and ρ is the equilibrium vapor density. The equilibrium vapor density is calculated by Kelvin's equation which assume equilibrium at a curved air-water interface as

where ρ is the saturated vapor density in the gas phase that varies with temperature as [Campbell, 1985]

High precision in equilibrium time is not expected to significantly improve the model prediction [Halder et al., 2011]. Therefore, we estimate teq using the characteristic time of pure binary diffusion for a simple cylindrical pore with an averaged pore characteristic length ℓ that the vapor has to diffuse

Equation (28) can therefore be rewritten as

2.3. Interface Boundary Conditions

To couple the two domains, suitable boundary conditions at the free and porous medium interface are needed. The different boundary conditions applied at the interface between porous medium and free medium are listed below [see also Mosthaf et al., 2011]:

2.3.1. Hydrodynamic Interface Boundary Condition

The continuity of the normal stresses (pressure jump boundary condition in the gas-phase), which is the sum of pressure term and viscous term in free fluid can be introduce as [Whitaker, 1999]

Here the term p in porous medium is defined as the sum of the water pressure and the capillary pressure. In practice, the viscous term may be small compared with the pressure, and in this case the continuity of total normal stress reduces to the continuity of pressure [Nield and Bejan, 2006].The continuity of the normal mass fluxes for the gas-phase can be expressed as

where the flux of the liquid phase has been vanished at the interface and

The Beavers-Joseph-Saffman slip boundary condition for the tangential component of the free flow velocity reads

where tj (j = 1, …, n−1) are linear independent unit tangential vectors to the boundary Γ, and α is a dimensionless slip coefficient. In this study, the surface roughness of the porous medium appears as a coefficient in the surface slip boundary condition. Therefore, the aerodynamic resistance of water vapor and surface resistance are not necessary in this model. As mentioned previously, the related humidity and temperature of the soil surface are outputs of this model. The parameter α is empirical and depends on flow conditions, interface location, surface microstructure and porosity. The slip coefficient can be estimated using an equation proposed by Neale and Nader [1974]. They imposed conditions of continuity of velocity and its gradient at the fluid porous medium boundary. By solving the Brinkman equation and comparing the resulting mass flow rate, they proposed that , where μ is the effective viscosity in Brinkman's model and μ is the fluid viscosity [Neale and Nader, 1974]. Following Whitaker [1999], the effective viscosity is explicitly given by , where ϕ the porous medium porosity. [Kim et al., 1994]The Brinkman equation contains Laplacian terms and is of the same order as Stokes equation. Therefore, through the use of Brinkman equation at the free/porous interface and in the bulk of porous matrix, the boundary conditions in both Stokes and Darcy equations become compatible. An averaging method used by Shavit et al. [2002] to study the free-flow problem at the interface of porous medium and free-flow region [Shavit et al., 2002]. They used a Cantor configuration representing the porous media. Their results show that the apparent viscosity approach does not produce a satisfactory agreement with experimental. They proposed therefore a modified Brinkman equation which fit better with flow problem. However, theoretical studies by Lundgren [1972] and Kim and Russel [1985] pointed out that the applicability of Brinkman equation is restricted only to the high porosity domain (ϕ > 0.6) where the curvatures of the streamlines in the bulk porous medium adjacent to the interface are generally greater than the pore diameter. This requirement is highly restrictive since most natural porous media have porosity <0.6. It is therefore preferable in most practical situations to use Darcy's law together with the Beavers and Joseph boundary condition [Nield, 2009].

In this study, we used a constant slip coefficient of 0.01 based on the only measurement data available in Kim et al. [1994] which has similar porous medium properties to our experimental study.

2.3.2. Continuity Boundary Condition for Temperature at the Porous Medium-Free Medium Interface

In this work, we assume a local thermal equilibrium at the interface, providing a continuity boundary condition for temperature as

2.3.3. Continuity Boundary Condition for Concentration at the Porous Medium-Free Medium Interface

We also assume the continuity of the component mass fluxes across the interface. This equilibrium is considered as

After the implementation of the initial and boundary conditions, the system of partial differential equations for the nonequilibrium model in a two-dimensional domain was simultaneously solved using the COMSOL Multiphysics software that is based on the Finite Element Method. A summary of the equations to be solved and their corresponding primary variables is listed in Table 1.

Table 1

List of Model Equations and Their Primary Variables

	Equation Type	Equation Numbers	Primary Variables
Free-flow subdomain	Continuity and Navier-Stokes	(2), (3)	vg, pg
	Component mass balance	(5)	wv
	Energy balance	(4)	T
Porous-medium subdomain	Mass balances	(9), (10)	pg, pℓ
	Component mass balance	(19)	wv
	Energy balance	(22)	T

3. Experimental Setup

To validate the proposed theoretical model, we developed a test system consisting of a wind tunnel placed above a soil tank equipped with a network of sensors to measure different soil-water variables. A series of experiments under varying boundary conditions were performed, using test sand for which the hydraulic and thermal properties were well characterized. Precision data for soil moisture, soil and air temperature and relative humidity, as well as wind velocity under well-controlled transient heat and wind boundary conditions were generated. In this section, the experimental material, methods and protocols are discussed.

3.1. Sand Material

Uniform specialty silica sand, Accusand 30/40 (effective sieve number) was used for this series of experiments. Based on the technical sheet provided by the manufacturer (Unimin Corp., Ottawa, MN), this sand is 99.8% quartz, its grain shape is classified as rounded, the uniformity coefficient is approximately 1.2, and the grain density is 2.66 g/cm3. Additional important properties of the test sand are summarized in Table 2. The capillary pressure (Pc)—water content (θ) relationship and thermal conductivity (λ)—water content (θ) relationship as a function of temperature for the sand under tight packing conditions were measured using a small Tempe cell apparatus with a network of sensors that continuously monitored soil moisture, temperature, capillary pressure, and soil thermal properties. A detailed description of the measurements of the Pc-θ and λ-θ relationships is given in Smits et al. [2010, 2012].

Table 2

Properties of Accusand 30/40 Under Tight Packing Conditions

					Van Genuchten Model Parameters (m = 1 −/n)d
d50 (mm)a	Dry Bulk Density (g cm−3)	Porosity	Residual wWater Contentb	Saturated Hydraulic Conductivity, Ks (cm s−1)c	α (m−1)	n
0.52	1.77	0.334	0.028	0.106	5.7	17.8

Estimated from sieve data provided by the manufacturer.

Measured in a separate column test.

Measured in a separate hydraulic conductivity test.

van Genuchten model parameters estimated using the computer code RETC.

3.2. Development of Experimental Apparatus

Laboratory testing was performed using a two-dimensional bench scale tank constructed of acrylic plexiglass (height = 55.0 cm, length = 25.0 cm, width = 9.0 cm, wall thickness = 1.25 cm, specific heat = 1464 J kg−1 K−1, thermal conductivity = 0.2 Wm−1 K−1, and density = 1150 kg m−3) as seen in Figure 3a (the inside volume of the tank, excluding the sensors, is 12,375 cm3). Water content and temperature distributions within the tank were continuously monitored using dielectric soil moisture sensors (Decagon Devices, Inc. ECH2O EC-5, sensor length = 5.5 cm, measurement frequency = 70 MHz, accuracy ±3%) and temperature sensors (Decagon Devices Inc. RT-1, accuracy ±0.5°C from 5 to 40°C, better than ±1.0°C from 40 to 50°C), respectively. In addition to the sensors within the tank, relative humidity and temperature was monitored at two locations on the soil surface, three locations 7.5 cm above the soil surface, and ambient conditions outside the tank using relative humidity/temperature sensors (Decagon Devices Inc. EHT RH/Temperature, accuracy ±2% from 5 to 90% RH, ±3% from 90 to 100% RH, temperature accuracy ±0.25°C). The relative humidity/temperature sensors at the soil surface were placed in direct contact with the soil grains in order to get an accurate reading directly on the soil surface.

Figure 3

Schematic view of the experimental setup: (a) tank and sensors and (b) complete setup including tank and ductwork (all dimensions are in centimeters).

A total of 25 soil moisture sensors and 19 temperature sensors were installed horizontally through the plexiglass walls of the tank. Five relative humidity/temperature sensors were installed on and above the soil surface (Figure 3a). Five channel, continuous data loggers (Decagon Devices, Inc. Em50) were used to collect data from all the sensors. During the experiment, the tank was placed on a scale (Sartorius Model 11209-95, Range = 65 kg, Resolution = ±1 g) in order to continuously monitor the cumulative weight loss of water from the tank.

To better understand the physical processes associated with evaporation in dry, bare soil conditions, heat flux, and air flow were induced over the soil surface to simulate an arid/semiarid environment. As shown in Figure 3b, heat and wind was generated by infrared heaters and a duct fan placed within galvanized steel ductwork (26.0 cm by 8.0 cm rectangular and 15.2 cm round diameter, respectively). The ductwork was aligned along the centerline of the tank to channel heated air across the soil surface. Five ceramic infrared heaters (Mor Electric Heating Assoc., Inc. Infrared Salamander Model FTE 500-240) placed in parallel within a reflector, were used to induce heat over the soil surface. An in-line duct fan (Suncourt Pro Model DB6GTP, 15.2 cm diameter) with variable speed controller (Suncourt Model VS200) was used to induce wind over the soil surface. As shown in Figure 3b, the infrared heaters were positioned upstream of the tank while the duct fan was positioned downstream in order to draw heated air from the upstream end over the soil surface and vent the air at the downstream end.

To maintain a constant air temperature over the soil surface, the infrared heaters were connected to a temperature control system (Chromalox Model 2104), which was regulated by an infrared temperature sensor (Exergen Corporation Model IRt/c.03) placed on the ductwork. The infrared sensor provides input back to the temperature control system to ensure that temperature output from the heaters remained constant. The ambient air temperature gradient at the outer boundaries of the tank was monitored using a relative humidity/temperature sensor (Decagon Devices Inc. EHT RH/Temperature sensor, accuracy as given above) upstream of the tank and an air temperature sensor (Decagon Devices Inc. ECT, sensor length = 3.0 cm, temperature accuracy as given above) downstream of the tank as seen in Figure 3b.

To control wind velocity over the soil surface, a variable speed controller, connected to the in-line duct fan, was used in combination with a galvanized steel duct damper (15.2 cm diameter) downstream of the fan and/or electrostatic register vent filters upstream of the fan. Wind velocity was monitored using a stainless steel pitot tube (Dwyer Instruments, Inc. Model 167-12, 0.32 cm diameter, 30.48 cm insertion length, accuracy ±5%) connected to a differential pressure transmitter (Omega Engineering, Inc. Model PX653-0.1D5V, Range = 0–0.1 inches of H2O, 0–25 Pascals) and an anemometer. The Pitot tube and anemometer were suspended in the center of the tank, 10 cm above the soil surface as seen in Figure 3a for the entire duration of the experiment. The pressure data from the Pitot tube were converted into wind velocity in meters per second by solving Bernoulli's equation and compiled using LabVIEW software (National Instruments Corp.). The wind velocity profile was used to obtain the average velocity over the soil surface.

3.3. Procedure

The tank was first wet-packed with Accusand 30/40 in incremental layers of approximately 1 cm in an effort to achieve uniform bulk density in accordance with the methods outlined in Sakaki and Illangasekare [2007]. The water table was initially established at the top surface of the tank. Prior to starting the experiment, the soil surface was covered with plastic wrap in order to prevent evaporation. No flow conditions were maintained along the bottom and side boundaries of the tank.

At the start of the experiment, the plastic wrap was removed from the upper boundary of the tank, to allow for evaporation, while heat and wind were induced at the soil surface using the infrared heaters and duct fan, respectively. For each experiment, air temperature was kept constant at approximately 35–40°C using the temperature control system (Chromalox Model 2104) regulated by the infrared temperature sensor (Exergen Corporation Model IRt/c.03) placed on the ductwork. Wind velocity varied for each experiment between 0.55 and 3.65 m/s. Water content, temperature, relative humidity, wind velocity, and soil tank weight were continuously monitored at 10 min intervals. A total of four experiments were conducted. Each experiment was run for approximately 12 days.

4. Results

In this section, we present a demonstration of experimental results (section 4.1.2) for a specific wind velocity. Then, the experimental results are compared in section 4.2.2 with those obtained from a numerical simulation for the horizontal temperature gradient setup shown in Figure 2.

4.1. Experimental Results Demonstration and Discussion

Results from experiment run 2 (average wind velocity of 1.22 m/s) are presented below including graphs and summary tables of measured relative humidity and temperature, as well as calculated saturation from measured soil water content. Data from this experiment are used as a base case for comparison with the remaining experimental data. Observed trends as well as differences between experiments are noted.

4.1.1. Wind Velocity

Figure 4 shows the wind velocity versus time for experiment run 2. There is a sinusoidal diurnal fluctuation in wind velocity, which may be due to variability in atmospheric conditions, e.g., barometric pressure changes to constant air density assumption. This diurnal trend was also observed in the other experiments and considered to be minor. Studying the impact of diurnal variation of temperature on evaporation is not the objective of this study; therefore, the wind velocity was averaged over the entire experiment and the mean wind velocity was used in the numerical model. The average measured maximum wind velocity for each experiment is summarized in Table 3.

Figure 4

Wind velocity over the soil surface for experiment run 2.

Table 3

Experimental Wind Velocities

Experiment Run	Average Maximum Wind Velocity (m/s)
1	0.55
2	1.22
3	3.00
4	3.65

4.1.2. Relative Humidity

Figure 5 displays the relative humidity and temperature measured on the soil surface versus time for experiment 2 (wind velocity = 1.2 m/s). Based on these data, initial relative humidity (RH) on the soil surface remains relatively constant (RH = 0.80) for a period of approximately 1.8 days, followed by a steep decrease over approximately 4 days at which point relative humidity stabilizes (RH = 0.35). In contrast, initial temperature (T) on the soil surface (T = 26°C) increases for a period of approximately 3 days at which point temperature stabilizes (T = 33°C). This general trend was observed in all the other experiments. Table 4 provides a summary of initial and final relative humidity and temperature measurements on the soil surface for each experiment. For the first 3 days of the experiment, the relative humidity of the downstream air is higher than the upstream air since the vapor flow to the free air stream moistures it (Figure 5). However, after 3 days, the trend is reversed and the RH of the upstream air is higher. The apparent reversed trend in relative humidity is due to either a difference in wall-soil contact at side of soil tank or a change in sensor contact with the soil (may have changed over time due to heating and flexibility of the sensor cable.

Figure 5

Relative humidity and temperature measured on the soil surface and relative humidity in ambient air upstream of the soil tank.

Table 4

Summary of Relative Humidity and Temperature Measured on the Soil Surface

Experiment Run	Maximum Wind Velocity (m/s)	Initial RH On Soil	Final RH On Soil	Initial Temperature On Soil (°C)	Final Temperature On Soil (°C)	Change in Temperature On Soil (°C)
1	0.55	0.75	0.35	27	31	+4
2	1.22	0.80	0.35	26	33	+7
3	3.00	0.65	0.25	29	37	+8
4	3.65	0.55	0.20	33	44.5	+11.5

Figure 6 displays relative humidity measured 7.5 cm above the soil surface versus time for experiment 2. Based on this data, relative humidity above the soil surface is approximately equal to the incoming relative humidity of the ambient air upstream of the tank, which typically varied between 0.05 and 0.10. Thus, no effect from evaporation was observed at this height above the soil surface; water evaporating from the soil surface was drawn downstream of the tank prior to reaching a height of 7.5 cm above the soil surface. This trend applies to all the experiments.

Figure 6

Relative humidity measured 7.5 cm above the soil surface at upstream, middle, and downstream of the tank.

4.1.3. Temperature

As previously discussed in sections 4 and 4.1, for each experiment, source temperature was kept constant at approximately 40°C using a temperature control system regulated by an infrared temperature sensor placed on the ductwork. Figure 7 displays the air temperature measured on the soil surface, 7.5 cm above the soil surface, upstream, in the middle and downstream of the tank.

Figure 7

Temperature measured 7.5 cm above the soil surface at upstream, middle, and downstream of the tank.

Based on this data, there is an apparent diurnal fluctuation in temperature, which may be due to changes in atmospheric conditions, variability in the temperature control system and/or temperature of the ductwork. There is also incremental heat loss from the upstream boundary to the downstream boundary of the tank. Observed heat loss from the upstream to downstream boundary may be influenced by the presence of the thermal sensors themselves; the thermal sensors have different thermal properties than the atmosphere, which could lead to perturbations of the thermal field. Overall heat loss across the tank was estimated by subtracting the upstream/incoming duct air temperature from the downstream/outgoing duct air temperature. A similar trend was observed in all experiments, as shown in Table 5.

Table 5

Summary of Ambient Air Temperature Upstream and Downstream of Tank

Experiment Run	Wind Velocity (m/s)	Air Temperature Upstream (°C)	Air Temperature Downstream (°C)	Change in Air Temperature (°C)
1	0.55	41	35	−6
2	1.22	41	35	−6
3	3.00	42	38	−4
4	3.65	46	41	−5

Figure 8a shows soil temperature versus time for experiment 2 for depths of 2.5 cm (Sensors 3), 7.5 cm (Sensors 8), 12.5 cm (Sensors 13), and 17.5 cm (Sensors 18). Sensors are numbered as labeled in Figure 3a. As evident in the figure, the deeper the sensor, the less influence the surface temperature and wind velocity have on the local temperature. We observed that soil temperature remained quite unchangeable with time below a depth of 15 cm (Sensors 16–25). Figure 8b shows temperature versus time for this experiment for a depth of 2.5 cm (Sensors 1, 3, and 5). We can observe here a difference in temperature for the same depth sensors. Sensor 3 is placed at the middle of the soil tank and, at this location, there is less heat loss compared to Sensor 1 and Sensor 5 located at the left and right side of the tank, respectively. Therefore, we observe a higher temperature for this sensor than two others. Finally, the temperature at the upstream is always higher than the downstream in the wind tunnel; therefore, Sensor 5 shows a higher temperature than Sensor 1. We will see in the next section that these differences create an asymmetrical saturation profile in the soil tank.

Figure 8

Evolution of soil temperature measured (a) vertically at the middle of soil tank and (b) horizontally at depths of 2.5 cm.

4.1.4. Saturation

Figure 9a shows the soil saturation measured over time for experiment 2 (U = 1.22 ms−1) for soil depths of 2.5, 7.5, 12.5, and 17.5 cm. For each experiment, the soil remained saturated below a depth of 15 cm (Sensors 16–25). Saturation values were calculated for each experiment using the empirical two-point α-mixing model from [Sakaki et al., 2008]. This method requires knowledge of the water saturated and air dry sensor analog-to digital converter (ADC) counts and soil porosity Although it is known that performance of dielectric permittivity sensors is temperature dependent due to the changes in the electrical characteristics of the soil with temperature [e.g., Assouline et al., 2010; Bogena et al., 2007; Kizito et al., 2008; Seyfried and Grant, 2007], the method developed by Sakaki et al. [2008] does not account for temperature dependency. Assouline et al. [2010] concluded that care must be taken in interpreting “subtle” changes in the apparent dielectric permittivity under conditions where temperature fluctuations are significant, for example, close to the soil surface under diurnal temperature conditions. Although the diurnal temperature variations in the experiments reported here were rather small (maximum of 4°C difference, 2.5 cm below the soil surface), we amended this method to correct for temperature sensitivity using the method outlined in [Cobos and Campbell, 2013]. The maximum temperature sensitivity of the measurements was relatively small (0.002 cm3 cm−3) but should be accounted for, as pointed out by Assouline et al. [2010], especially when analyzing soil water profiles in the shallow subsurface. Bogena et al. [2007] used a two-step calibration method to develop preliminary models for estimating dielectric constant as a function of temperature [Bogena et al., 2007]. According to Bogena et al. [2007], for the soil temperature ranges reported here within (approximately between 23 and 27°C), the maximum error in soil water content due to temperature effects on the sensor circuitry is 0.4 vol. %. Seyfried and Grant [2007] found that for temperatures ranging from 5 to 45°C, the maximum apparent water content change is ±0.028 cm3 cm−3 for saturated soils while the effect of temperature on dry samples was not detectable [Seyfried and Grant, 2007].

Figure 9

Evolution of soil saturation measured (a) vertically at the middle of soil tank and (b) horizontally at a depth of 2.5 cm.

Based on these saturation curves, (Figure 9a), it appears that the primary, capillary-driven drying front associated with Stage-1 evaporation moves relatively fast through the first row of sensors (i.e., Sensor 3 located 2.5 cm under the soil surface) and reaches the second row, at a depth of 7.5 cm, by the first day. This drying front moves slower through this second row of sensors, i.e., more gently sloping saturation curve for Sensors 6–10 compared to Sensors 1–5, as evaporation appears to transition to Stage-2 evaporation by about the second day. The primary drying front reaches the third row of sensors, at a depth of 12.5 cm, by the third day. Based on the gentle slope of the saturation curves by the third day, it appears that there is a secondary, diffusion-driven drying front, associated with Stage-2 evaporation, which is slowly transporting water to the soil surface. This general trend was observed in all the experiments; more gently sloping saturation curves were noted for experiments with wind velocity <1.22 m/s and more steeply dipping saturation curves for those experiments with wind velocity >1.22 m/s. Figure 9b shows saturation versus time for experiment 2 for a depth of 2.5 cm (Sensors 1, 3, and 5). We can observe a light asymmetrical distribution of saturation in the soil tank near the soil surface. This asymmetrical behavior may come from the difference in temperature upstream and downstream of the wind tunnel and as discussed above, slightly affect the saturation readings. The kinetic of the cumulative evaporation is plotted in Figure 10 for four different free flow maximum average wind speeds (U) of 0.55, 1.22, 3.00, and 3.65 m/s. This figure shows clearly that the wind speed has a considerable effect on the evaporation processes. Here increasing the wind speed increases the total evaporation. Increasing wind speed has a great effect on the initial evaporation rate then its influence becomes smaller. These observations are indicative of the two distinct stages of bare soil evaporation. Stage-1 evaporation is an atmosphere-controlled stage in which the evaporation rate is relatively high and relatively constant due to high atmospheric demands (i.e., high temperatures and wind velocity at the soil surface) and predominantly independent of soil water content [Lehmann et al., 2008; Shokri et al., 2010]. A transition regime, evaporation is driven by capillary transport and occurs at or near the rate of free-water evaporation [Bittelli et al., 2008; Dingman, 2002; Yiotis et al., 2003; Lehmann and Or, 2009]. However, experimental evidence suggests that Stage-1 evaporation may not always be high and constant but rather drop from the onset of the drying process. This drop, as seen in the experimental data presented here, is often associated with high wind velocities with a thin boundary layer and large soil pores [Shahraeeni et al., 2012]. Stage-2 evaporation, also known as the falling rate period, is a soil-controlled stage in which the evaporation rate is relatively low and controlled by the rate at which water can be transmitted to the soil surface in response to potential gradients induced by upward-decreasing soil water contents, rather than atmospheric demands. This evaporation is driven by diffusive transport and is less than the rate of free-water evaporation [Dingman, 2002; Lehmann et al., 2008; Yamanaka et al., 2004]. As anticipated, due to the greater atmospheric demand, i.e., faster wind velocity, experiments with wind velocity 1.22 m/s, the constant relative humidity (or Stage-1 evaporation) was approximately 2–3 days, whereas, experiments with wind velocity >1.22 m/s, Stage-1 evaporation was approximately 1–1.5 days.

Figure 10

Kinetic of cumulative evaporation for different wind speeds (U = maximum wind speed in free medium).

From the comparison of the slopes of the curves (Figure 10), we can see that the evaporation rate is faster for a higher wind speed for the first stage of the evaporation processes. However, the wind speed has less influence over time in the second stage (diffusion-dominant) evaporation. Increasing wind speed, increases the first stage evaporation rate and decreases the transition time from first stage to second stage evaporation. For high values of wind speed, evaporation becomes less dependent on changes in wind velocity. The evaporation rate in the diffusion-dominant regime seems to decrease for high values of wind speed. The origin of this diffusion flux reduction may be due to the turbulence at the soil surface. In order to prove this phenomenon, further experiments are needed. This issue will be investigated in future work.

In the next section, we compare experimental results with our numerical model to better understand any discrepancies between theory and experiments and the validity of the proposed model.

4.2. Numerical Simulation and Comparison With Experimental Results

To validate the proposed two-dimensional, nonisothermal, nonequilibrium coupled evaporation model, the numerical results for temperature, saturation, and evaporation rate are compared with experimental results. The subdomains, porous media properties, boundary conditions, and initial values are depicted in Figure 11. Here the heat source is imposed in the same direction as the experimental wind speed (right-hand side of the wind tunnel). In reality, changing wind speed in the experimental setup can change the heat loss rate in the free medium and influence the evaporation rate. Therefore, before comparing the experimental and numerical results, we first show the sensitivity of the model to some parameters, e.g., free-flow thermal and solutal dispersion, wind speed, wind temperature, and vapor concentration, Beavers-Joseph slip coefficient, soil permeability and porosity, and soil residual water content, respectively. In this parametric study, we have neglected the heat loss in the free-flow domain. This approach also allows us to use very small wind velocities that are difficult to generate experimentally.

Figure 11

Subdomains, boundary conditions, and initial values of experimental case.

4.2.1. Effect of Free-Flow Thermal and Solutal Dispersion

Because the Péclet numbers that correspond to the heat and mass transfer in the free-flow region of our experimentation are high, we cannot ignore thermal and solutal dispersion. An average wind speed of 1 m/s, with free-flow height of H = 0.21 m (see the experimental configuration of Figure 3) results in Pe = 11235 and Pe = 9975 from equation (6). From equations (6) and (7), the ratio of the effective coefficients to the molecular diffusion and conduction coefficients becomes very significant. These high values clearly show the importance of thermal and solutal dispersion in the free medium domain in this study. To show the importance of including solutal and thermal dispersion in the free-flow domain on the evaporation processes, simulated cumulative evaporation for U = 1 m/s with and without considering disperison are plotted in Figure 12. As we can see, there is a significant difference between the two cases, demonstrating the importance of including atmospheric system dispersion in evaporation modeling efforts.

Figure 12

Comparison between evaporation processes with and without taking dispersion in the free-flow domain into account (K = 1 × 10−10 m2, Swr = 0.075, U = 1 m/s, TU = 313 K, wvU = 0.006, αBJ = 0.01, φ = 0.33).

To better understand the effect of dispersion in the free-flow regime on the evaporation rate, a comparison was made between the results from two models, with and without dispersion. Surface plots of temperature and concentration at day 1 are depicted in Figure 13 for the two models. Comparing the temperature profiles of Figures 13a and 13b (with and without considering dispersion) shows that the thermal dispersion causes an increase in temperature with a more uniform heat distribution on the soil surface, and therefore, increasing evaporation. Here the temperature distribution on the soil surface is asymmetric without considering the dispersion (Figure 3a). However, temperatures measured on the upstream and downstream soil surface do not show a significant asymmetrical behavior (see Figure 5). Comparing Figures 3a and 3b also shows how the dispersion causes the replacement of the humid air on the soil surface with the dry air flowing through the wind tunnel. Therefore, the air demand for evaporation increases and the evaporation process is enhanced.

Figure 13

Temperature and vapor concentration surface plot with and without considering free-flow dispersion (K = 1 × 10−10 m2, Swr = 0.075, U = 1 m/s, TU = 313 K, wvU = 0.006, αBJ = 0.01, φ = 0.33).

A vertical profile of temperature and vapor concentration at the middle of the soil tank is plotted in Figure 14 at one and 11 days, for two cases (with and without considering the dispersion in free-flow region). We can see clearly that considering the dispersion in free flow increases the temperature and decreases the vapor concentration on the soil surface, which enhances the evaporation processes.

Figure 14

Comparisons of (a) temperature and (b) vapor concentration profile between with and without considering dispersion after t = 1 and 11 days for soil profile at tank center (K = 1 × 10−10 m2, Swr = 0.075, U = 1 m/s, TU = 313 K, wvU = 0.006, αBJ = 0.01, φ = 0.33).

4.2.2. Effect of Wind Speed

Figure 15 shows simulation results for a range of wind speeds between 0.002 and 6 m/s. This figure clearly shows that the wind speed has a considerable effect on the drying processes in nonisothermal system. Increasing the wind speed increases the first stage evaporation rate and decreases the transition time between two evaporative stages (soil water flow regime to water vapor diffusion controlled regime) at low velocity values; then, at high values of wind speed the evaporation rate becomes less dependent and finally independent on wind speed. On the contrary, the impact of wind speed on second stage evaporation (diffusion-dominant stage) is not significant.

Figure 15

Impact of wind speed on evaporation processes in a nonisothermal system (K = 1 × 10−10 m2, Swr = 0.075, TU = 313 K, wvU = 0.006, αBJ = 0.01, φ = 0.33).

Increasing the wind speed in a nonisothermal system increases the heat and mass transfer by convection in free-flow region. The thermal and solutal dispersion becomes more and more important. Therefore, the temperature is increased and the vapor concentration becomes the same as the wind vapor concentration above the soil surface.

As the experimental results also show in Figure 10, one can distinguish a critical velocity value from which the evaporation is no longer dependent on the wind speed. At the critical wind speed, the temperature rises to the wind temperature and the vapor concentration falls to the wind vapor concentration above the soil surface. Therefore, a higher wind speed cannot change the evaporation processes.

4.2.3. Effect of Wind Temperature and Vapor Concentration on Evaporation

The effect of wind temperature on evaporation processes from the soil was tested by comparing three theoretical models, T_U = 22°C (isothermal case), 30°C and 40°C (experimental reference case). Figure 16 provides a comparison of the cumulative evaporation with these different wind temperature levels. Increasing the wind temperature enhances the evaporation from the soil tank and decreases the transition time from a capillary-dominated evaporation regime to the diffusion-dominated regime. The phase change rate increases with an increase in soil surface temperature. As a consequence, the evaporation increases, the interface dries out faster and the transition happens earlier. However, according to Figure 16, changing the wind temperature does not significantly change the evaporation rate during Stage-2 evaporation.

Figure 16

Effect of wind temperature on drying process (K = 1 × 10−10 m2, Swr = 0.075, U = 1 m/s, wvU = 0.006, αBJ = 0.01, φ = 0.33).

The variations of the wind water vapor concentration are tested in Figure 17 for the vapor mass fractions of 0.000, 0.015, and 0.030, respectively. Here the drier the air blowing across the soil surface, the higher the cumulative evaporation in both stages. Similar to temperature effects, the soil surface dries faster with a lower vapor concentration, and therefore, transition between the two evaporation phases occurs earlier. The vapor concentration gradient and air capacity for evaporation is higher for the dryer air blowing on the soil surface; therefore, the evaporation process enhanced.

Figure 17

Effect of wind vapor quantity on drying process (K = 1 × 10−10 m2, Swr = 0.075, U = 1 m/s, TU = 313 K, wvU = 0.006, αBJ = 0.01, φ = 0.33).

4.2.4. Effect of Slip Coefficient (α)

The influence of the Beavers-Joseph coefficient on the drying process is examined by varying α in the range of 0.01–10. Results for cumulative evaporation versus time (Figure 18) show that the Beavers-Joseph slip coefficient, α, slightly influences the resulting cumulative evaporation but not significantly. This observation is consistent with the results in Baber et al. [2012]. Therefore, the choice of the Beavers-Joseph coefficient has a small influence on the evaporation rate across the interface.

Figure 18

Effect of Beavers-Joseph slip coefficient on drying process (K = 1 × 10−10 m2, Swr = 0.075, U = 1 m/s, TU = 313 K, wvU = 0.006, φ = 0.33).

4.2.5. Effect of Porous Medium Permeability (k)

We tested three different porous media permeabilities: k = 10−9 m2, 10−10 m2 and 10−11 m2. As depicted in Figure 19, changing the permeability does not change the first stage evaporation rate significantly. However, increasing the permeability increases the transition time between two evaporation stages, which affects the cumulative evaporation.

Figure 19

Effect of porous medium permeability on drying process (Swr = 0.075, U = 1 m/s, TU = 313 K, wvU = 0.006, αBJ = 0.01, φ = 0.33).

4.2.6. Effect of Porous Medium Porosity (ϕ)

The cumulative evaporation for different porous medium pore volume fractions (porosity), over 11 days is plotted in Figure 20. As these results show, porosity has a great influence on evaporation process in porous media.

Figure 20

Effect of porous medium porosity on drying process (K = 1 × 10−10 m2, Swr = 0.075, U = 1 m/s, TU = 313 K, wvU = 0.006, αBJ = 0.01).

Increasing the porous medium pore volume fraction increases considerably the amount of the cumulative water evaporation from the porous medium. It accelerates the evaporation rate in the capillary-driven evaporation regime and enhances the diffusive flux in the diffusion-dominated regime. For example, increasing the porosity of porous medium of 0.10–0.33 can increase the amount of the evaporation about four times.

4.2.7. Effect of Porous Medium Residual Saturation (Swr)

The sensitivity of the model is also examined by varying the water residual saturation (Swr) in the porous medium. The results of evolution of cumulative evaporation are plotted in Figure 21, for different Swr values. Here changing the residual saturation does not change the first stage evaporation rate significantly. The transition time occurs earlier for higher water residual saturation. However, increasing the residual saturation decreases the evaporation rate in the second evaporation stage considerably. Higher soil residual saturation causes a higher tortuosity factor at the diffusion-dominant regime which reduces the diffusion flux to the soil surface and decreases the evaporation rate.

Figure 21

Effect of residual saturation on drying process (K = 1 × 10−10 m2, U = 1 m/s, TU = 313 K, wvU = 0.006, αBJ = 0.01, φ = 0.33).

4.2.8. Impact of the Sensors

The sensors themselves have different thermal properties than the porous medium, which could lead to perturbations of the heat and mass transfer in the system. To study the impacts of the sensors on the evaporation process, we simulated the theoretical model with a 2-D geometry crossing on the soil temperature sensors. Each sensor was included in the model as a square, 0.075 m of length, corresponding to the sensor diameter. The soil moisture sensors are stainless steel and completely water proof with a thermal conductivity of 16 W/(m K), density of 798 kg/m3 and heat capacity of 520 J/(kg K). Figure 22 shows the results of the cumulative evaporation for the model with and without consideration of the sensors. Results show that the soil sensors increase the first stage evaporation rate by approximately 2%. This is principally due to the sensors improving the thermal conduction in the porous medium. Although the sensors enhance the heat transfer by conduction in the soil, they also favor heat loss from the soil. When considering the sensor in the theoretical model, the evaporation rate decreases about 5% for the second stage evaporation. The sensors decrease the tortuosity of the porous medium, and therefore, increase the diffusion flux of the water vapor from the soil to the soil surface. The comparison results in the next section come from the model geometry without the sensors. However, the results are modified according to the impact of the sensor percentage calculated in this section.

Figure 22

Impact of the sensors on drying process (K = 1 × 10−10 m2, Swr = 0.075, U = 1.2 m/s, TU = 313 K, wvU = 0.006, αBJ = 0.01, φ = 0.33)

In the next section, comparing these results with experimental results is done to understand the amount of discrepancy between theory and experience and the validity of the proposed model.

4.2.9. Comparison With Experimental Results

Observed and simulated evaporation rates for wind speeds of 0.55 and 1.2 m/s are plotted in Figures 23a and 23b, respectively. We note that the flow regime in the free-flow region is lamniar for these two wind speed values. The free-flow wind mean temperature and mean concentration are chosen from the experimental data and imposed as a boundary condition in the free-flow region upstream. Generally, there is good agreement between model and experimental results, although in both cases, the model underestimated the evaporation rate for the regime transition period of the experiments. Although not shown, this was seen in the 3.0 and 3.6 m/s test cases as well. In all comparisions, however, the model could capture the transition behavior between Stage-1 and Stage-2 evaporation.

Figure 23

Comparison between experimental and theoretical evaporation rate results for U = 0.55 m/s and U = 1.2 m/s.

The evolution over time of the observed saturation (shown before in Figure 9) is compared with the simulated results in Figures 24a and 24b. Comparison between experimental and theoretical saturation results for U = 1.2 m/s at x = 0.125 m for different depths is plotted in Figure 4a. The model can predict the first and second stage evaporation with a good precision at a depth of 2.5 cm below the soil surface. However, the results for the second stage evaporation at a depth of 7.5 and 12.5 cm are over estimated. The origin of this discrepancy is that the model overestimates the temperature variation in soil at this stage of evaporation.

Figure 24

Comparison between experimental and theoretical saturation results for U = 1.2 m/s at (a) x = 0.125 m and (b) depth of 0.025 m.

The simulated and experimental results of the saturation evolution are plotted in Figure 24b for U = 1.2 m/s at 2.5 cm below the soil surface. The numerical results of the saturation profiles do not show the asymmetric behavior across the soil tank horizontal axis.

Figures 5a and 5b show the comparison between predicted and measured temperature and saturation profiles after 6.5 days at different depths, located along the centerline of the test tank. The observed and modeled temperatures and saturations agreed well for all locations over the entire duration of the experiments. This is statistically confirmed with the MIA values (range from 0.923 to 0.988). The deviations between simulated and measured soil moistures and temperatures may, in part, be due the accuracy and resolution of the soil moisture sensors and thermistors compared to the model. The last experimental data point in Figure 25a corresponds to the temperature measured on the soil surface (see Figure 5).

Figure 25

Simulated and observed (a) temperature and (b) saturation at t = 6.5 days for soil profile at tank center (arc-length 0.125 m).

Using the measured soil surface relative humidity and temperature data, we calculated the water vapor concentrations on the soil surface and compared them with corresponding vapor concentrations from the numerical simulation at differents depths. As seen in Figure 26, simulation results 1 cm beneath the soil surface fit the experimental soil surface results better than than the simulation results from the soil surface. In an attempt to understand the relative humidity directly on the soil surface, the relative humidity sensors require good contact with the the gas phase in the soil pores so that the readings can better reflect that of the soil surface and not the surrounding air. It can be argued that in the process of maintaining good contact, the sensor is not reading the soil surface humidity but rather the relative humidity of some location beneath the soil surface. The model is able to correctly capture the kinetic form of the phemomena; however, there is a small dicrepency between the measured and simulated profiles. We note that, in this model, the evaporation rate and surface relative humidity are a direct outcome of the proposed model concept and not imposed as a boundary condition.

Figure 26

Simulated and observed vapor concentration with time on the soil surface (observed and simulated) and at locations 1, 2, and 3 cm below the soil surface for U = 1.22 m/s.

5. Conclusion

With the goal of improving our understanding of the land/atmospheric coupling to predict the important process of bare soil evaporation, we developed a model based on the coupling of Navier-Stokes free flow and Darcy flow in porous medium. The model consists of the coupled equations of mass conservation for the liquid phase (water) and gas phase (water vapor and air) in porous medium with gas phase (water vapor and air) in free-flow domain under nonisothermal, nonequilibrium phase change conditions. The boundary conditions at the porous medium-free flow medium interface include dynamical, thermal and solute equilibriums, and using the Beavers-Joseph slip boundary condition. What is unique about this coupled model is that the evaporation rate and soil surface temperature conditions come directly from the model output. In order to experimentally validate the numerical results, we developed and used a unique two-dimensional test system with a wind tunnel placed above a soil tank equipped with a network of different sensors. The theoretical model sensitivities to the important soil and air properties are discussed. The parameter sensitivity study shows that the drying processes can be physically interpreted and qualitatively match the behavior of the physical experiments. Results from numerical simulations were compared with experimental data.

The comparison results demonstrate that the coupling concept used in the integrated model formulation can predict the different stages of the drying process in porous media with good accuracy. Increasing the wind speed, in a nonisothermal system, increases the first stage evaporation rate and decreases the transition time at low velocity values; then, at high values of wind speed the evaporation rate becomes less dependent of flow in free fluid. In the opposite, the impact of the wind speed on the second stage evaporation rate (diffusion-dominant stage) is not significant. The thermal and solutal dispersion in free-flow domain should be taken into account when solving the coupled concept model. The new phase change rate equation leads to good correlation between theoretical and experimental results. The proposed theoretical model can be used to predict the evaporation process where a porous medium flow is coupled to a free flow for different practical applications.

Greek Symbols

ϕ Volume fraction of the pore or porosity.

μ Effective viscosity, Pa.s.

μ Dynamic viscosity for the α-phase, Pa.s.

ρ Vapor density in gas phase, kg/m3.

ρ Equilibrium vapor density in gas phase, kg/m3.

ρ Total mass density for the α-phase, kg/m3.

τ Tortuosity of porous medium.

α Beavers-Joseph coefficient.

Thermal dispersion tensor for moist air in free medium, W/m.K.

λ Thermal conductivity of the α-phase, W/m.K.

Effective thermal conductivity coefficient in porous medium, W/m.K.

Del Operator.

Γ Interface between free-flow medium and porous medium.

ω weighting factor in equation (24).

Subscripts, Superscripts, and Other Symbols

α ℓ (liquid), s (solid) and, g (gas) phases.

pm Porous medium.

ff Free-flow medium.