Published online 26 May 2006
Published in Vadose Zone J 5:784-800 (2006)
DOI: 10.2136/vzj2006.0007
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Numerical Analysis of Coupled Water, Vapor, and Heat Transport in the Vadose Zone
Hirotaka Saitoa,*,
Jiri 
im
neka and
Binayak P. Mohantyb
a Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92521
b Biological and Agricultural Engineering, Texas A&M Univ., College Station, TX 77843-2117
* Corresponding author (hirotaka.saito{at}ucr.edu)
Received 12 January 2006.
 |
ABSTRACT
|
|---|
Vapor movement is often an important part in the total water flux in the vadose zone of arid or semiarid regions because the soil moisture is relatively low. The two major objectives of this study were to develop a numerical model in the HYDRUS-1D code that (i) solves the coupled equations governing liquid water, water vapor, and heat transport, together with the surface water and energy balance, and (ii) provides flexibility in accommodating various types of meteorological information to solve the surface energy balance. The code considers the movement of liquid water and water vapor in the subsurface to be driven by both pressure head and temperature gradients. The heat transport module considers movement of soil heat by conduction, convection of sensible heat by liquid water flow, transfer of latent heat by diffusion of water vapor, and transfer of sensible heat by diffusion of water vapor. The modifications allow a very flexible way of using various types of meteorological information at the soil–atmosphere interface for evaluating the surface water and energy balance. The coupled model was evaluated using field soil temperature and water content data collected at a field site. We demonstrate the use of standard daily meteorological variables in generating diurnal changes in these variables and their subsequent use for calculating continuous changes in water contents and temperatures in the soil profile. Simulated temperatures and water contents were in good agreement with measured values. Analyses of the distributions of the liquid and vapor fluxes vs. depth showed that soil water dynamics are strongly associated with the soil temperature regime.
Abbreviations: DOY, day of the year • TDR, time domain reflectometry
|
INTRODUCTION
|
|---|
THE simultaneous movement of liquid water, water vapor, and heat in the vadose zone plays a critical role in the overall water and energy balance of the near-surface environment of arid or semiarid regions in many agricultural and engineering applications. Moisture near the soil surface is influenced by evaporation, precipitation, liquid water flow, and water vapor flow, most of which are strongly coupled. In arid regions, vapor movement is often an important part of the total water flux since soil moisture contents near the soil surface usually are very low (Milly, 1984).
In agricultural applications, spatial and temporal changes in surface soil moisture need to be well understood to achieve efficient and optimum water management. Vapor transport is also important since the actual contact area between liquid water and seeds is often very small such that seeds need to imbibe water from vapor to germinate (Wuest et al., 1999).
In engineering applications, accurate assessment of the volume of leachate from landfills, which is essential for long-term landfill management, requires estimates of water and heat movement through engineered landfill covers (e.g., Khire et al., 1997). Understanding surface energy and water balances as well as liquid water, water vapor, and heat movement in soils is critical for the performance evaluation of engineered surface covers for waste containment in landfills in arid or semiarid regions (Scanlon et al., 2005). Traditional engineered covers typically consist of multilayered resistive barriers with low hydraulic conductivities to minimize water movement into the underlying waste. Recently, evapotranspiration (ET) covers have been developed (e.g., Hauser and Gimon, 2004) to increase the water storage capacity of landfill covers and their protective function by removing water through evapotranspiration. These ET cover designs thus need to be assessed mainly in terms of the near surface water and energy balance.
Another engineering application requiring an assessment of coupled liquid water, water vapor, and heat transport involves nuclear waste repositories. Radioactive decay in repositories may create steep temperature gradients in adjacent soils and rocks, which can lead to unwanted consequences. Generated heat may cause evaporation of soil water, and subsequent migration and condensation of water vapor in cooler areas. This process may significantly change the physical characteristics of surrounding materials due to the precipitation and dissolution of various minerals (e.g., Spycher et al., 2003). Finally, in military applications, an accurate prediction of water contents and temperatures near the soil surface may allow improved detection of buried land mines since their presence can significantly change these two variables (imnek et al., 2001).
Early pioneering studies on interactions between liquid water, water vapor, and heat movement were reported by Philip and de Vries (1957), who provided a mathematical description of liquid water and water vapor fluxes in soils driven by both pressure head (isothermal) and soil temperature (thermal) gradients. They derived the governing flow equation for nonisothermal flow as an extension of the Richards equation, which originally considered only the pressure head gradient. The theory of Philip and de Vries (1957) was later extended by Nassar and Horton (1989), who additionally considered the effect of an osmotic potential gradient on the simultaneous movement of water, solute, and heat in soils.
Heat transport and water flow are coupled by the movement of water vapor, which can account for significant transfer of latent energy of vaporization. Soil temperatures may be significantly underestimated when the movement of energy associated with vapor transport is not considered. For example, Cahill and Parlange (1998) reported that 40 to 60% of the heat flux in the top 2 cm of a bare field soil of Yolo silt loam was due to water vapor flow. Fourier's law describing heat transport due to conduction (e.g., Campbell, 1985) thus needs to be extended to include heat transport by liquid water and water vapor flow. The general heat transport model then considers movement of soil heat by (i) conduction, (ii) convection of sensible heat by liquid water flow, and transfer of (iii) latent and (iv) sensible heat by diffusion of water vapor (e.g., Nassar and Horton, 1992).
The soil–atmosphere interface is an important boundary condition affecting subsurface movement of liquid water, water vapor, and heat under field conditions. Direct measurement of all components needed to fully evaluate the water and energy balance at the soil surface, such as precipitation, evapotranspiration, heat flux, and radiation, are rarely available at short time intervals. In many situations, only standard daily meteorological data from weather stations are available. If detailed simulations of water and heat fluxes are required, the components of the surface mass and energy balances need to be calculated from standard daily information, such as daily maximum and minimum air temperatures, daily maximum and minimum relative humidities, the shortwave incoming solar radiation, and daily average wind speed. One common approach to address this issue is to first calculate approximate diurnal changes in the meteorological variables using relatively simple models (Ephrath et al., 1996) and to then use estimated values for further investigation. Very few studies have demonstrated and evaluated such approaches for simulating soil moisture and temperature in field soils.
Although it is widely recognized that the movement of liquid water, water vapor, and heat are closely coupled and strongly affected by each other, their mutual interactions are rarely considered in practical applications. The effect of heat transport on water flow is often neglected, arguably because of model complexity and a lack of data to fully parameterize the model, but partly also because of a lack of user-friendly simulation codes. The two major objectives of this study thus were to develop a numerical model that (i) solves the coupled equations governing liquid water, water vapor, and heat transport, together with the surface water and energy balance, and (ii) provides considerable flexibility in accommodating various types of meteorological information that can be collected at selected time intervals (daily, hourly, or other time intervals). Complete evaluation of the movement of liquid water, water vapor, and heat in the subsurface can be accomplished by simultaneously solving the system of equations for the surface water and energy balances, subsurface heat transport, and variably saturated flow. These equations need to be solved numerically since they are highly nonlinear and strongly coupled. This coupled movement of liquid water, water vapor, and heat in the subsurface, as well as interactions of these subsurface processes with the energy and water balances at the soil surface, was implemented in the HYDRUS-1D code (imnek et al., 1998a). The HYDRUS-1D is a widely used, well-documented, and tested public domain code for simulating water and solute transport in soils (e.g., Scanlon et al., 2002).
Model performance of the modified HYDRUS-1D code was tested in this study by simulating continuous changes in water contents, temperatures, and fluxes of a bare field soil. Results were compared against field soil temperature and water content data collected at different depths at a field site near the University of California Agricultural Experimental Station in Riverside, CA, during the fall of 1995. We demonstrate how standard daily meteorological variables can be used for simulating diurnal changes in soil water contents, temperatures, and liquid water, water vapor, and heat fluxes. In addition, the impact of vapor flow on predictions of soil water contents and soil temperatures was investigated.
|
MATERIALS AND METHODS
|
|---|
Models
Liquid Water and Water Vapor Flow
The governing equation for one-dimensional flow of liquid water and water vapor in a variably saturated rigid porous medium is given by the following mass conservation equation:
 | [1] |
where 
is the total volumetric water content (m3 m–3), qL and qv are the flux densities of liquid water and water vapor (m s–1), respectively, t is time (s), z is the spatial coordinate positive upward (m), and S is a sink term usually accounting for root water uptake (s–1). The total volumetric water content is defined as
 | [2] |
where l and v are volumetric liquid water and water vapor contents (expressed as an equivalent water content, m3 m–3).
The flux density of liquid water, qL, is described using a modified Darcy law as given by Philip and de Vries (1957):
 | [3] |
where qLh and qLT are isothermal and thermal liquid water flux densities (m s–1), respectively, h is the pressure head (m), T is the temperature (K), and KLh (m s–1) and KLT (m2 K–1 s–1) are the isothermal and thermal hydraulic conductivities for liquid-phase fluxes due to gradients in h and T, respectively.
The flux density of water vapor, qv, can also be separated into isothermal, qvh, and thermal, qvT, vapor flux densities (m s–1) as follows (Philip and de Vries, 1957):
 | [4] |
where Kvh (m s–1) and KvT (m2 K–1 s–1) are the isothermal and thermal vapor hydraulic conductivities, respectively. Combining Eq. [1], [3], and [4], we obtain the governing liquid water and water vapor flow equation:
 | [5] |
where KTh (m s–1) and KTT (m2 K–1 s–1) are the isothermal and thermal total hydraulic conductivities, respectively, and where:
 | [6] |
 | [7] |
Soil Hydraulic Properties
The pore-size distribution model of Mualem (1976) was used to predict the isothermal unsaturated hydraulic conductivity function, KLh(h), from the saturated hydraulic conductivity and van Genuchten's (1980) model of the soil water retention curve, Eq. [9]:
 | [8] |
where Ks is the saturated hydraulic conductivity (m s–1), Se is the effective saturation (unitless), and l and m are empirical parameters (unitless). The parameter l was given a value of 0.5 as suggested by Mualem (1976). The parameter m was determined by fitting van Genuchten's analytical model (van Genuchten, 1980)
 | [9] |
to water retention data using the RETC code (van Genuchten et al., 1991). In Eq. [9], s and r are the saturated and residual water contents (m3 m–3), respectively, and 
(m–1), n (unitless), and m (=1 – 1/n) are empirical shape parameters.
The thermal hydraulic conductivity function, KLT, in Eq. [5] is defined as follows (e.g., Noborio et al., 1996b):
 | [10] |
where GwT is the gain factor (unitless), which quantifies the temperature dependence of the soil water retention curve (Nimmo and Miller, 1986), 
is the surface tension of soil water (J m–2), and 0 is the surface tension at 25°C (=71.89 g s–2). The temperature dependence of is given by
 | [11] |
where is in g s–2 and T in °C.
The isothermal, Kvh, and thermal, KvT, vapor hydraulic conductivities are described as (e.g., Nassar and Horton, 1989; Noborio et al., 1996b; Fayer, 2000)
 | [12] |
 | [13] |
where D is the vapor diffusivity in soil (m2 s–1), 
w is the density of liquid water (kg m–3), sv is the saturated vapor density (kg m–3), M is the molecular weight of water (M mol–1, =0.018015 kg mol–1), g is the gravitational acceleration (m s–2, =9.81 m s–2), R is the universal gas constant (J mol–1 K–1, =8.314 J mol–1K–1), 
is the enhancement factor (unitless, Cass et al., 1984), and Hr is the relative humidity (unitless). The vapor diffusivity, D, of the soil is defined as
 | [14] |
where a is the air-filled porosity (m3 m–3), 
is the tortuosity factor as defined by Millington and Quirk (1961):
 | [15] |
and Da is the diffusivity of water vapor in air (m2 s–1) at temperature T (K):
 | [16] |
The saturated vapor density, sv (kg m–3), as a function of temperature is expressed as
 | [17] |
while the relative humidity, Hr, can be calculated from the pressure head, h, using a thermodynamic relationship between liquid water and water vapor in soil pores as (Philip and de Vries, 1957):
 | [18] |
When the liquid and vapor phases of water in soil pores are in equilibrium, the vapor density of the soil can be expressed as the product of the saturated vapor density and the relative humidity. The HYDRUS-1D code uses an enhancement factor, , to describe the increase in the thermal vapor flux as a result of liquid islands and increased temperature gradients in the air phase (Philip and de Vries, 1957). An equation for the enhancement factor was derived by Cass et al. (1984) and is expressed as
 | [19] |
where fc is the mass fraction of clay in the soil (unitless).
Heat Transport
The governing equation for the movement of energy in a variably saturated rigid porous medium is given by the following energy conservation equation:
 | [20] |
where Sh is the storage of heat in the soil (J m–3), qh is the total heat flux density (J m–2 s–1), and Q accounts for sources and sinks of energy (J m–3 s–1). The storage of heat is given by
 | [21] |
where T is a given temperature (K), n is the volumetric fraction of solid phase (m3 m–3), Cn, Cw, Cv, and Cp are volumetric heat capacities (J m–3 K–1) of the solid phase, liquid water, water vapor, and moist soil (de Vries, 1963), respectively, and L0 is the volumetric latent heat of vaporization of liquid water (J m–3) given by
 | [22] |
where Lw is the latent heat of vaporization of liquid water (J kg–1; Lw [J kg–1] = 2.501 x 106 – 2369.2T [°C]).
The total heat flux density, qh, is defined as the sum of the conduction of sensible heat as described by Fourier's law, sensible heat by convection of liquid water and water vapor, and latent heat by vapor flow (de Vries, 1958):
 | [23] |
where 
() is the apparent thermal conductivity of soil (J m–1s–1K–1), and qL and qv are the flux densities of liquid water and water vapor (m s–1), respectively.
Combining Eq. [20], [21], and [23] results in the governing equation for heat movement (e.g., Nassar and Horton, 1992; Fayer, 2000):
 | [24] |
where the last term on the right side represents a sink of energy associated with root water uptake.
Soil Thermal Properties
The apparent thermal conductivity of soil, (), in Eq. [23] combines the thermal conductivity of the porous medium in the absence of flow and the macrodispersivity, which is assumed to be a linear function of velocity (de Marsily, 1986; Hopmans et al., 2002). The apparent thermal conductivity, (), may then be expressed as (e.g., imnek and Suarez, 1993):
 | [25] |
where ß is the thermal dispersivity (m). Since the thermal dispersivity plays an important role only when the liquid water flux is very large, only a few studies have actually determined its values (e.g., Hopmans et al., 2002). The thermal conductivity, 0(), accounts for the tortuosity of the porous medium, and can be described with a simple equation given by Chung and Horton (1987):
 | [26] |
where b1, b2, and b3 are empirical regression parameters (W m–1 K–1). Chung and Horton (1987) provided average values for the b coefficients for three textural classes (i.e., clay, loam, and sand), which are implemented in the HYDRUS-1D code. Soil thermal properties are much more influenced by water content than by textural differences since thermal properties of the different mineral components of the solid phase are all approximately of the same order of magnitude (Jury and Horton, 2003).
Initial inspection of the fully coupled set of equations for liquid water, water vapor, and heat transport may suggest that significantly more parameters are needed compared with more approximate decoupled water flow and heat transport models. This actually is not the case. The soil hydraulic properties are fully described by six soil-specific hydraulic parameters that appear in the van Genuchten–Mualem model (i.e., r, s, , n, Ks, and l); these parameters are needed to simulate variably saturated liquid water flow based on the Richards equation. The soil thermal properties are described by means of three b coefficients that appear in the thermal conductivity function. Considering additionally thermal liquid flow, and both thermal and isothermal vapor flow, does not increase the demand for input parameters since KLT, Kvh, KvT, and Sh (Eq. [10], [12], [13], and [21]) can be fully described using literature values of various physical properties, such as surface tension, the diffusivity of water vapor, and the heat capacities of the different soil components.
Surface Energy Balance
Surface precipitation, irrigation, evaporation, and heat fluxes are used as boundary conditions for liquid water and water vapor flow and heat transport in field soils. Evaporation and heat fluxes can be calculated from the surface energy balance (e.g., van Bavel and Hillel, 1976; Noborio et al., 1996a; Boulet et al., 1997):
 | [27] |
where Rn is net radiation (W m–2), H is the sensible heat flux density (W m–2), LE is the latent heat flux density (W m–2), L is the latent heat (J kg–1), E is the evaporation rate (kg m–2s–1), and G is the surface heat flux density (W m–2). While Rn and G are positive downward, H and LE are positive upward.
Net radiation is defined as (e.g., Brutsaert, 1982)
 | [28] |
where Rns is the net shortwave radiation (W m–2), Rnl is the net longwave radiation (W m–2), a is the surface albedo (unitless), St is the incoming (global) shortwave solar radiation (W m–2), 
s is the soil surface emissivity (unitless) representing the reflection of longwave radiation at the soil surface, Rld
is the incoming (thermal) longwave radiation at the soil surface (downward flux) (W m–2) as emitted by the atmosphere and cloud cover, and Rlu
is the sum of the outgoing (thermal) longwave radiation emitted from the surface (vegetation and soil) into the atmosphere (W m–2).
The surface albedo can be described with a simple linear expression relating the albedo with the surface water content since it depends, especially for bare soils, on the soil surface wetness. Van Bavel and Hillel (1976) proposed the following simple formulae to calculate the surface albedo from the surface wetness:
 | [29] |
where top is the water content at the surface.
Using the Stefan–Boltzmann law (e.g., Monteith and Unsworth, 1990), the net longwave radiation can be rewritten as (e.g., Brutsaert, 1982):
 | [30] |
where the subscripts a and s are used for variables of the atmosphere and the soil, respectively. The emissivity of the atmosphere depends both on the air temperature and humidity, while the emissivity of the soil surface depends on the water content and the vegetation at the soil surface. The emissivity of bare soil can be expressed as a function of the volumetric water content (van Bavel and Hillel, 1976) as
 | [31] |
leading to an emissivity of 0.9 for a dry surface and 
0.98 when the soil is saturated ( 0.45). Following Idso (1981), the atmospheric emissivity can be calculated as follows:
 | [32] |
where ea is the atmospheric vapor pressure (kPa), which can be expressed as a function of Ta as
 | [33] |
For a sky partially covered with clouds, the longwave radiation is calculated with (Monteith and Unsworth, 1990)
 | [34] |
where c (unitless) is the fraction of cloud cover (or the cloudiness factor) of the sky. This equation was derived from the average temperature difference between the atmosphere and the clouds in England (Monteith and Unsworth, 1990). The cloudiness factor is readily estimated from the atmospheric transmission coefficient for solar radiation, Tt, as (Campbell, 1985)
 | [35] |
A value of the incoming shortwave solar radiation St (W m–2) at any given time and location can be calculated by taking into account the position of the sun (e.g., van Bavel and Hillel, 1976) using the following equation (Campbell, 1985):
 | [36] |
where Gsc is the solar constant (1360 W m–2), and Tt (unitless) is defined as the ratio of the measured daily global solar radiation Stm (W m–2) and the daily potential global (extraterrestrial) radiation Ra (Wm–2):
 | [37] |
The last term of Eq. [36], e (rad), is the solar elevation angle given by (Monteith and Unsworth, 1990)
 | [38] |
where t is the local time within a day and t0 is the time of solar noon.
Among a number of available accepted models (e.g., Wang and Bras, 1998), the sensible heat flux, H, can be simply defined as (e.g., van Bavel and Hillel, 1976)
 | [39] |
where Ta is the air temperature (K), Ts is the soil surface temperature (K), Ca is the volumetric heat capacity of air (J m–3 K–1), and rH is the aerodynamic resistance to heat transfer (s m–1).
Since surface evaporation is controlled by atmospheric conditions, surface moisture, and moisture transport in the soil, a model is needed that accounts for all of these factors. Evaporation, E, is often calculated as
 | [40] |
where vs is the water vapor density at the soil surface (kg m–3), va is the atmospheric vapor density (kg m–3), rv is the aerodynamic resistance to water vapor flow (s m–1), and rs is the soil surface resistance to water vapor flow that acts as an additional resistance along with the aerodynamic resistance (s m–1) (Camillo and Gurney, 1986). Evaporation models commonly use only the aerodynamic resistance in the denominator of Eq. [40] (e.g., Milly [1984] and some existing codes, such as SHAW [Flerchinger et al., 1996] and UNSATH [Fayer, 2000]). This approach may not be correct since it overestimates the evaporation rate for dry soils because of the assumption of Eq. [18] that equilibrium exists in the soil (Kondo et al., 1990). When the soil surface is dry, water vapor in larger pores is dynamically transported toward the atmosphere and its density in pores is not in equilibrium with the average water content of a given depth. The soil surface resistance is thus used to account for an extra resistance to water vapor flow in soil pores. This resistance depends strongly on both soil structure and texture. The relationship between the soil surface resistance, rs, and the surface water content has been empirically formulated in many studies and typically has an exponential form (e.g., van de Griend and Owe, 1994), indicating that the surface resistance increases dramatically as the soil dries out. The soil surface resistance in this study is determined using the following formula (Camillo and Gurney, 1986):
 | [41] |
Calculation of the sensible heat flux, H, or evaporation rate, E, using Eq. [39] or [40] requires knowledge of the aerodynamic resistances to heat transfer, rH, and to water vapor flow, rv, respectively. These may be calculated from the surface wind speed as follows (Campbell, 1985):
 | [42] |
where U* is the friction velocity (m s–1), k is the von Karman constant, which has a value of 0.41, zref is the reference height of temperature measurements (m), zH is the surface roughness for the heat flux (m), d is the zero-plane displacement (m), and 
H is the atmospheric stability correction factor for the heat flux (unitless). The assumption that the resistance to heat flux is equal to that to vapor flow has long been used (e.g., van Bavel and Hillel, 1976; Flerchinger et al., 1996). The friction velocity in Eq. [42] is defined as (Campbell, 1985)
 | [43] |
where u is the mean wind speed (m s–1) at height zref, zm is the surface roughness for the momentum flux (m), and m is the atmospheric stability correction factor for the momentum flux (unitless). For bare soils, d is equal to zero, while typical surface roughness values of 0.001 m are used for both zH and zm (Oke, 1978). In general, calculations of rH and rv require an iterative procedure since stability correction factors must be evaluated from variables that themselves depend on the stability correction factors. Camillo and Gurney (1986) proposed an approximate approach to calculate the stability correction factor, which is used also in this study so no need exists for the iterative procedure.
Meteorological Variables Generation
Solving the energy balance equation (Eq. [27]) at a time interval of interest requires values of meteorological variables at the same or similar time intervals. Weather stations, however, do not always provide standard meteorological data at time intervals of interest. Although several models exist for calculating diurnal changes from daily average values (e.g., Ephrath et al., 1996), finding the best model was not an objective of this study. We used relatively simple approaches to generate continuous values of the meteorological variables from available daily information.
When the daily maximum and minimum air temperatures are the only information provided by a weather station, continuous estimates of the air temperature, Ta, can be obtained using a trigonometric function with a period of 24 h as follows (Kirkham and Powers, 1972):
 | [44] |
where 
is an average daily temperature (°C), AT is the amplitude of the cosine wave (°C) that can be calculated from the difference between the daily maximum and minimum temperatures, and t is a local time within the day. The argument of the cosine function shows that the highest temperature is assumed to occur at 1300 h and the lowest at 0100 h.
Since daily temperature variations typically show a cyclic behavior throughout the day, it is reasonable to assume that the relative humidity also shows such a cyclic pattern. Similarly as for air temperature (Eq. [44]), a trigonometric function can be used to calculate continuous values of the relative humidity from daily information (e.g., Gregory et al., 1994). When both the daily maximum and minimum relative humidity values are available, continuous values of the relative humidity can be generated as well using a cosine function with a period of 24 h as follows:
 | [45] |
where 
is the average daily relative humidity (unitless), AH is the amplitude of the cosine wave that can be calculated from the difference between the maximum and minimum relative humidity values, t is a local time within the day, and tmax is the hour of the day when the relative humidity is at its maximum. The relative humidity is generally at its minimum during daytime when the air temperature and the wind speed are at their maxima. Gregory et al. (1994) showed that the relative humidity was highest at around 0500 and 0600 h in Lubbock, TX. In this study, we used 0500 h for the peak time of the relative humidity. The diurnal change in the saturation vapor density may be calculated from generated temperature values using Eq. [17], while the corresponding atmosphere vapor density, va, can be calculated by multiplying the generated relative humidity (Eq. [45]) and the saturation vapor density.
It is well known that wind transports heat and effectively mixes the soil–atmospheric boundary layer. Wind is generally highly variable in space and time since it involves mostly turbulent flow and as such is characterized by random fluctuations in its speed and direction (Campbell, 1977). Because many meteorological models require continuous inputs of the wind speed, continuous values of the wind speed must somehow be calculated from daily information. In this study, we used the following simple approach based on the maximum-to-minimum wind speed ratio Ur, which is defined as
 | [46] |
where Umin and Umax (m s–1) are the unknown minimum and maximum wind speeds of the day, respectively. This ratio may be determined from prior knowledge or calibrated on available data. The maximum and minimum wind speeds can then be calculated from the daily average wind speed as follows:
 | [47] |
 | [48] |
where 
(m s–1) is a daily average wind speed. Cyclic behavior of the wind speed U during the day is then modeled using the daily average wind speed as follows (Gregory, 1989):
 | [49] |
where t is the time of the day, and tmax is the time when the maximum wind speed occurs. In Gregory's (1989) study, the maximum wind speed occurred at 1530 h.
Calculation of continuous values of the incoming shortwave radiation using Eq. [36] requires continuous values of the transmission coefficient, Tt (Eq. [37]). When only the daily incoming shortwave solar radiation, Stm, is available, a daily average transmission coefficient must be used since no reliable generation is available. Calculation of the daily average transmission coefficient using Eq. [37] requires values of the daily potential global solar radiation, which depends on the latitude of the location of interest and the day of year as follows (Campbell, 1985):
 | [50] |
where dr is the relative distance from the earth to the sun, 
s is the sunset hour angle (rad), 
is the angular distance of the sun (i.e., solar declination, rad), and 
is the latitude of the location (rad).
Field Data
The coupled numerical model and other meteorological modules implemented in HYDRUS-1D were tested using a soil temperature and moisture data set collected at the field site near the University of California Agricultural Experimental Station in Riverside, CA. Temporal soil temperature and water content variations were measured using horizontally inserted thermocouples and time domain reflectometry (TDR) probes, respectively, near the soil–air interface at depths of 2, 7, and 12 cm during the fall of 1995 (24 November (DOY 328)– 5 December (DOY 339); see Mohanty et al., 1998). The TDR probes consisted of three 20-cm-long rods spaced 2.5 cm apart (total width 5 cm). The experimental field was irrigated twice during the measurement period on DOY 334 and 335 with 0.55 and 0.20 cm of water, respectively, using a sprinkler system containing several laterals and outlet ports. Irrigation rates were measured using catch cans near the measuring location. Soil temperatures and water contents at each depth were measured at 20- and 40-min intervals, respectively.
Field Measurements
Figure 1
shows measured water contents and temperatures at three depths at the measurement location. Soil water contents were not recorded from noon of DOY 330 to noon of DOY 331 due experimental problems. The soil temperature data show a typical sinusoidal diurnal behavior at all three depths, while the soil water content data do not reveal a distinct diurnal pattern, only noisy fluctuations. When TDR is used to measure the soil water content, noisy fluctuations are often inherent in the interpretation of TDR waveforms (e.g., Cahill and Parlange, 1998). Irrigations on DOY 334 and 335 increased the water contents at depths of 2 and 7 cm, while the water content remained almost constant at the depth of 12 cm during the experiments. As a result of irrigation, average water contents before and after irrigation were significantly different (P < 0.01), even at the depth of 7 cm. More details on temperature and water content variations at the site can be found in Mohanty et al. (1998).
View larger version (50K):
[in this window]
[in a new window]
|
Fig. 1. Soil temperatures and water contents measured at three depths (2, 7, and 12 cm) from 24 November (Day of the Year [DOY] 328) to 5 December (DOY 339) at the experimental site in Riverside, CA.
|
|
Soil Hydraulic Data
The soil at the experimental site was characterized as an Arlington fine sandy loam (coarse-loamy, mixed, thermic Haplic Durixeralf). Soil water retention curves were measured in the laboratory using the Tempe Cell (Soil Moisture Equipment Corp.) method for pressure heads down to –8 m and a pressure chamber (Soil Moisture Equipment Corp.) for lower pressure heads (down to –150 m). Figure 2
shows a plot of the water retention curve for a soil core sample collected 4 m east of the measurement location. The figure includes both measured drainage (drying) and imbibition (wetting) data. The van Genuchten (1980) analytical model (Eq. [9]) was fitted to the retention data using the RETC code (van Genuchten et al., 1991), leading to r = 0.011 m3 m–3, s = 0.445 m3 m–3, = 0.0277 cm–1, and n = 1.38.
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 2. Water retention curve measured on a soil core sample collected at the experimental site in Riverside, CA. Both drainage and imbibition curves are plotted (open symbols: drainage, closed symbols: imbibition). The van Genuchten (1980) model (VG) is fit to the measured data.
|
|
Several studies previously measured the saturated hydraulic conductivity, Ks, of the Arlington fine sandy loam soil using different methodologies (e.g., imnek et al., 1998b; Wang et al., 1998a, 1998b). Table 1 lists soil hydraulic conductivities obtained by the different researchers. Measured values vary depending on the methodology used. In this study, we used the value of 34.2 cm d–1 obtained by Wang et al. (1998a) using a tension infiltration experiment for two reasons. First, their measurements were conducted at the same site as used in our study. Also, their Ks value was close to the median value obtained with the three different methods listed in Table 1. The clay fraction of the Arlington fine sandy loam at Riverside needed for the enhancement factor (Eq. [19]) was 8.8% (imnek et al., 1998b).
View this table:
[in this window]
[in a new window]
|
Table 1. Soil hydraulic conductivities of Arlington fine sandy loam estimated by different researchers using different methodologies. Listed values are all geometric means of measured values.
|
|
Meteorological Data
Standard daily meteorological data in California can be obtained from approximately 400 weather stations that are a part of the California Irrigation Management Information System (CIMIS) Network operated by the California Department of Water Resources (DWR). Data collected daily by DWR were taken from the University of California Integrated Management Program website (http://www.ipm.ucdavis.edu/WEATHER/about_weather.html, accessed November 2005, verified 26 Mar. 2006). Data were obtained for the automatic CIMIS weather station located at UC Riverside, Agricultural Operations Department Area, CA (33°58'N, 117°20'W, elevation 310.9 m). Daily meteorological variables used in this study are summarized in Table 2. Data were also available from the CIMIS website (http://wwwcimis.water.ca.gov/cimis/welcome.jsp, accessed November 2005, verified 26 Mar. 2006), from which we downloaded daily and hourly values of standard meteorological variables. We used hourly values from this website to validate the various generations of continuous weather parameters (see also Table 2 for variables used in this study).
View this table:
[in this window]
[in a new window]
|
Table 2. Weather variables recorded at the University of California-Riverside Agricultural Operations Department Area.
|
|
Numerical Implementation
The total variably saturated water flow and heat transport equations (Eq. [5] and Eq. [24]) were both solved numerically using the finite element method for spatial discretization and finite differences for the temporal discretization. The HYDRUS-1D software package (imnek et al., 1998a) was modified to consider the simultaneous movement of liquid water, water vapor, and heat, as well as the surface energy and water balances. The modified version of HYDRUS-1D either can consider continuous values of net radiation and other meteorological variables measured at any time interval, which are then used directly in the energy balance equation (Eq. [27]), or can calculate continuous values from daily meteorological data as described above.
At each time step, the net radiation (Rn), the sensible heat flux (H), and the latent heat flux (LE) were calculated to obtain the surface heat flux (G) by solving the surface energy balance equation (Eq. [27]), which is subsequently used as a known heat flux boundary condition. Calculation of Rn, H and LE requires knowledge of the temperature and pressure head of the soil surface. These were obtained by sequentially solving the governing equations of the surface energy balance, variably saturated water flow, and heat transport twice, while updating temperatures and water contents between the two runs.
Numerical solutions of transport equations often exhibit undesired numerical oscillations. An appropriate combination of space and time discretization can often prevent such oscillations. In the original HYDRUS-1D, undesired oscillations were minimized or eliminated by using two dimensionless numbers that characterize the space and time discretization, i.e., the grid Peclet and Courant numbers (imnek et al., 1998a). A dimensionless Courant number is used also in the modified HYDRUS-1D code to stabilize the numerical solution of the heat transport equation. Neglecting terms associated with the transfer of latent heat by water vapor in Eq. [24], the Courant number, Cu, may be defined as
 | [51] |
where 
t and z are temporal and spatial steps, respectively. The maximum permitted time step, tmax, is then calculated assuming a value of 1 for the Courant number.
The numerical solution can also be stabilized by specifying the maximum nodal temperature difference, Tmax, allowed during a given time step (e.g., 0.25°C). The maximum allowed time step is calculated as
 | [52] |
where Ti is the maximum nodal temperature change at the previous time step t. The maximum permitted time step for the heat transport calculations is then the smaller value of Eq. [51] or [52].
Initial and Boundary Conditions
The soil profile was considered to be 50 cm deep, with observation nodes located at depths of 2, 7, and 12 cm for comparing calculated temperatures and volumetric water contents with measured values. A constant nodal spacing of 2 mm was used, leading to 251 discretization nodes across the problem domain. Calculations were performed for a period of 12 d from 24 Nov. (DOY 328) to 5 Dec. (DOY 339) 1995. No natural precipitation was recorded at the weather station during this time period.
Initial water contents (a constant value of 0.13 m3 m–3) and soil temperatures (variable with depth) were determined from measured values on DOY 328. Boundary conditions at the soil surface for liquid water, water vapor, and heat transport were determined from the surface water and energy balance equations. When the soil was irrigated on DOY 334 and 335 (0.55 and 0.20 cm of water, respectively), irrigation rates were used as the flux boundary condition. For heat transport, a heat flux at the soil surface, G, was calculated from the surface energy balance equation (Eq. [27]). The temperature of irrigation water was assumed to be equal to the air temperature. Zero pressure and temperature gradients were used as the bottom boundary conditions at the depth of 50 cm. These conditions assume that the water table is located far below the domain of interest and that heat transfer across the lower boundary occurs only by convection of liquid water and water vapor.
|
RESULTS AND DISCUSSION
|
|---|
Meteorological Model Verification
Generated continuous changes in the air temperature, relative humidity, and wind speed are compared with hourly measured values obtained from the nearby weather station in Fig. 3
for the entire simulation period. A conventional value of 3 was used for the maximum-to-minimum wind speed ratio, Ur (Eq. [46]). Most calculated values compared well with hourly measured values, even for the wind speed, which fluctuated more than the other variables due to inherent randomness. Although the general trends of the relative humidity were well described, some maximum values were slightly overestimated. Using the generated continuous air temperatures and relative humidities, diurnal variations in the atmosphere vapor pressure were calculated using Eq. [17] and compared with hourly vapor pressures obtained from the weather station (Fig. 4
). As shown in Fig. 4, diurnal changes in the vapor pressure calculated from the generated meteorological variables closely fitted the hourly vapor pressure values at the weather station. Results presented in Fig. 3 and 4 demonstrate the validity of using the generated continuous diurnal variations in meteorological variables.
View larger version (31K):
[in this window]
[in a new window]
|
Fig. 3. Diurnal changes of three meteorological variables—air temperature, relative humidity, and wind speed—generated from daily information during the simulation period along with hourly values measured at the weather station near the study site from 23 November (Day of the Year [DOY] 327) to 6 December (DOY 340).
|
|
View larger version (19K):
[in this window]
[in a new window]
|
Fig. 4. Vapor densities calculated from generated air temperatures (Eq. [44]) and relative humidities (Eq. [45]) along with vapor pressures calculated hourly at the weather station from 23 November (Day of the Year [DOY] 327) to 6 December (DOY 340).
|
|
Meteorological variables, including net radiation, air temperature, and relative air humidity, were measured at the experimental site every 5 min from 21 Oct. (DOY 294) to 1 Nov. (DOY 305) and every 10 min from 2 Nov. (DOY 306) to 25 Nov. (DOY 329) 1995 before the soil temperature and water content were measured. To investigate whether or not generated diurnal changes in the meteorological variables can be used to reproduce continuous changes in the net radiation, calculations were performed from DOY 294 to 329 with the same domain conditions as described above. Daily meteorological data during this period were again obtained from the weather station in Riverside (Table 2). Diurnal changes in the meteorological variables needed to calculate the net radiation at any given time were calculated from the daily data.
Figure 5
shows temporal changes in both the measured and simulated net radiations for a 10-d period. Positive values indicate an incoming flux to the land surface (downward), while negative radiation values imply an outgoing flux from the land surface (upward). Calculated and measured daytime net radiations matched reasonably well, whereas nighttime agreement was less good, especially around 0000 h of DOY 310, 311, 315, and 319. This may have been caused by clouds or fog on these days, which reduced the net radiation toward zero. The current model takes cloudiness into account only during the daytime (Eq. [35]); modeling such effects during nights is not easy. Figure 5 also depicts simulated net shortwave and longwave radiations, which were used to calculate net radiation values. Simulated net radiation values before and after the displayed period in Fig. 5 matched measured values similarly well. Results presented in Fig. 5 demonstrate that our choice of meteorological models to calculate the net radiation is reasonable.
View larger version (33K):
[in this window]
[in a new window]
|
Fig. 5. Net radiation calculated from 6 Nov. 1995 (Day of the Year [DOY] 310) to 16 Nov. 1995 (DOY 320) using the standard daily meteorological data obtained from the weather station at the UC-Riverside Agricultural Operations Department Area in Riverside, CA. Measured hourly net radiation values were obtained at the experimental site near the University of California Agricultural Experimental Station in Riverside. Net radiation was calculated as the sum of the net incoming shortwave solar radiation and net outgoing longwave radiation.
|
|
Simulated Soil Temperatures and Water Contents
Soil water contents and soil temperatures at the experimental site were numerically simulated from DOY 328 to 339 using the coupled liquid water, water vapor, and heat transport module implemented in the HYDRUS-1D code. Figure 6
shows measured and simulated water contents at depths of 2, 7, and 12 cm. Predicted water contents follow fairly well the measured values at all three depths during the entire simulation period, including the gradual decrease in the water content at the 2-cm depth before irrigation. The observed noise in measured values is inherent to the TDR measurements (e.g., Cahill and Parlange, 1998). Rapid increases in the water content at the 2-cm depth after 0.55 and 0.20 cm of irrigation water was applied on DOY 334 and 335, respectively, were predicted reasonably well, with a slight overestimation after the first irrigation and some underestimation after the second irrigation. The small increase in at the 7-cm depth after the two irrigations was also well predicted. The simulated and measured water contents at the 12-cm depth were almost constant during the simulation period.
View larger version (25K):
[in this window]
[in a new window]
|
Fig. 6. Simulated and measured volumetric water contents during the simulation time period (Day of the Year [DOY] 328–DOY 339) at three depths (2, 7, and 12 cm) at the experimental site.
|
|
Figure 7
depicts simulated and measured soil temperatures at three depths. The simulated and measured temperatures both show typical sinusoidal diurnal behavior. The amplitude of both the simulated and measured daily temperature variations decreased with depth due to attenuation of the transported heat energy. Before the irrigation on DOY 334, the simulated temperature amplitude was larger than the observed one. Jury and Horton (2003) presented an analytical solution for the temperature amplitude at a particular depth z:
 | [53] |
where A0 and A are temperature amplitudes at the soil surface and the specified depth (°C), respectively, is the angular frequency (s–1), and KT is the thermal diffusivity (m2 s–1), being the ratio of the thermal conductivity and the heat capacity. The temperature amplitude can be reduced by either increasing the soil heat capacity (Eq. [21]) or decreasing the soil thermal conductivity, both of which would lead to a lower thermal diffusivity. For example, reducing the thermal diffusivity by 50% would lead to a much better fit of the measured temperatures at all three depths. After the irrigation, simulated values almost perfectly reproduced the measured ones. Apart from the temperature amplitude, simulated and measured temperatures generally agreed at all three depths, with most simulated values being within a few degrees of the measured values. A large number of formulas can be used to calculate various meteorological parameters, such as air emissivity or boundary layer resistances. For example, Mahfouf and Noilhan (1991) showed that, depending on the choice of an evaporation model, simulated soil surface temperatures can be easily altered by 3°C or more. Differences between simulated and measured values in our study decreased with depth and were <2°C at the 12-cm depth.
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 7. Simulated and measured soil temperatures during the simulation time period (Day of the Year [DOY] 328–DOY 339) at three depths (2, 7, and 12 cm) at the experimental site.
|
|
Calculated surface energy fluxes during the simulation period are depicted in Fig. 8
. Although the net radiation values were fairly constant during the simulation period, variations in daily surface heat and sensible heat fluxes were more pronounced. This confirms that the dynamics of the surface energy components and the partitioning of the energy are influenced by many different land and atmospheric attributes and as such cannot be predicted solely from the net radiation. While surface heat fluxes are needed as boundary conditions for subsurface heat transport, their variations are directly related to simulated soil temperature variations (Fig. 7). Latent heat fluxes are in general very small during dry periods because of a lack of soil moisture for evaporation (i.e., before irrigation); however, after the soil was irrigated on DOY 334 and 335, the latent heat fluxes increased significantly during the next few days as soil moisture near the surface increased. The peak of the latent heat flux was predicted to be approximately at 1200 h every day as expected. The sensible heat flux decreases after irrigation since a large part of the incoming energy is then used for evaporation (Fig. 8). The direction of the sensible heat flux varies depending on the time of day; the soil surface is heated during daytime, leading to upward sensible heat fluxes (positive values), while the sensible heat flux in general is downward during night (negative values) when the soil surface cools down.
View larger version (23K):
[in this window]
[in a new window]
|
Fig. 8. Simulated surface energy balance components at the experimental site: net radiation, sensible heat flux, latent heat flux, and surface heat flux from 24 November (Day of the Year [DOY] 328) to 5 December (DOY 339). All fluxes are positive upward, except for the net radiation where positive value indicates incoming radiation to the soil surface (downward) and negative value outgoing radiation from the soil surface (upward).
|
|
The temporal variation in the latent heat flux may be investigated also by analyzing temporal variations in the pressure head and relative humidity at the soil surface (Fig. 9
). As soil water evaporates into the atmosphere during the daytime, the pressure head of the soil surface is reduced to as low as –10 000 cm (pF = 4). The pressure head (or relative humidity) then increases sharply during the nighttime. The mechanisms of increases in the pressure head and relative humidity during the night are discussed below.
View larger version (22K):
[in this window]
[in a new window]
|
Fig. 9. Temporal changes of the pressure head and relative humidity at the soil surface during the simulation period from 24 November (Day of the Year [DOY] 328) to 5 December (DOY 339). The corresponding surface relative humidity was calculated from the surface pressure head using Philip and de Vries' equilibrium model (Eq. [18]).
|
|
Figure 10
shows temporal variations in the calculated transmission coefficient, the surface albedo, the surface emissivity, and the air emissivity during the simulation period. All of these variables are necessary to calculate the net radiation. The transmission coefficient and air emissivity show large temporal variations. Variations in the air emissivity were due to changes in the vapor pressure (Eq. [32]), temporal variations of which are also distinct (Fig. 4). The soil surface albedo and soil emissivity, on the other hand, were fairly constant during the simulation period, except for the albedo immediately after irrigation on DOY 334 and 335. Wet surfaces do not reflect radiation as much as dry surfaces, resulting in higher net radiation values during and after the irrigations (Fig. 7). Figure 10 demonstrates that temporal variations in most variables are not simple and need to be well estimated before they should be used in calculations of the net radiation.
View larger version (16K):
[in this window]
[in a new window]
|
Fig. 10. Temporal changes of transmission coefficient, surface albedo, surface emissivity, and atmosphere emissivity required for calculations of the net radiation at a given time of day for the simulation time period (Day of the Year [DOY] 328–DOY 339). These variables are calculated directly by the program.
|
|
Liquid Water and Water Vapor Fluxes
Calculated distributions of the liquid water and water vapor fluxes vs. depth for both isothermal and thermal components before and after irrigation are depicted in Fig. 11
and 12
, respectively. While the plots in Fig. 11 show fluxes on a typical dry day (from 1200 h of DOY 329 to 1200 h of DOY 330), the impact of irrigation on DOY 335 on these fluxes is shown in Fig. 12. The thermal liquid flux was almost always negligible compared with the other three fluxes on the dry day.
View larger version (16K):
[in this window]
[in a new window]
|
Fig. 11. Calculated vertical distributions of the thermal and isothermal fluxes of liquid water and water vapor at 1200 and 0000 h of a typical dry 24-h period during Days of the year [DOY] 329 and 330. Positive values indicate upward and negative values downward fluxes.
|
|
TOP
ABSTRACT
INTRODUCTION
