Numerical Simulation of Effect of Sand Thickness on Soil Evaporation

Evaporation from the soil is an important component of evapotranspiration, and mulching greatly affects soil evaporation. We conducted numerical simulations to study the effect of the thickness of sand mulch on soil evaporation. We tested nine treatments: mulching with sand thicknesses of 1, 3, 5, 6, 8, 10, 15 and 20 cm plus an unmulched control (CK). Accumulated evaporation was significantly lower, and the resistance to evaporation was significantly higher, for the mulched treatments than CK. The volumetric soil water content (SWC) was significantly higher for the mulched treatments than CK, but SWC varied little for thicknesses >5 cm. Heating was slower and more uniform for the mulched treatments than for CK. With the increase of the thickness of sand, the temperature transmission was slowed down. The change of soil temperature was not obvious at thicknesses >15 cm. A thickness of 5 cm was the most effective for storing water and preserving heat. Our results provide a theoretical basis and technical guidance for the effective use and management of soil water in farmland in arid regions.


INTRODUCTION
Soil evaporation is the diffusion of soil water into the air in the form of vapour and is a special component of the hydrological cycle and also an important component of evapotranspiration (Morillas et al. 2013). Soil evaporation represents an important loss of soil water on farms in semi-arid regions, a key factor restricting agricultural planting (Mellouli et al. 2000, Danierhan et al. 2013, Lei et al. 2014. Mulching can effectively restrain the evaporation of farmland water, thus preserving soil water (Li 2003, Yuan et al. 2009). Mulching with sand is a common and effective practice with a long history for decreasing soil evaporation (Fairbourn et al. 1973, Nachtergaele et al. 1998, Zhao et al. 2017a. Studying soil evaporation under sand mulching is therefore important for the management of farmland water.
Many studies have reported that adding sand mulching to the soil surface can effectively restrict the evaporation of soil water. Yamanaka et al. (2004) found that resistance to the flow of water vapour increased exponentially with the thickness of a layer of gravel mulch on the soil surface. Gravel-sand mulches are more effective than bare soil in conserving soil water, and the water content increases with mulch thickness (Ma et al. 2011). Kaseke et al. (2012) reported that net and total evaporation were 2.03 and 3.42-fold higher from bare soil than from soil mulched with gravel. Gravel mulch substantially increases resistance to evaporation, and mulched soils have much lower accumulated evaporation (Diaz et al. 2005, Qiu et al. 2014. A numerical simulation is an efficient approach for studying changes to soil water and temperature during soil evaporation (Assouline et al. 2014, Balugani et al. 2018. The variables can be effectively controlled to ensure the accuracy of the test. Simulations can also address more complex problems than experimentation. VADOSE/W is finite-element software based on a soil-atmosphere model that mainly models infiltration, soil evaporation, temperature and water balance. It has been used with increasing popularity during the last decade for predicting the infiltration and evaporation of soil water (Weeks et al. 2006). Li et al. (2016) found that VADOSE/W could simulate the behaviour of the flow of water in unsaturated loessial soil, and Zhang et al. (2016) found that VADOSE/W could accurately simulate and predict unsaturated flow associated with capillary barrier covers. Zhao et al. (2017b) reported that VADOSE/W could reliably simulate soil evaporation under sand mulching and sand inclusion.
At present, sand mulching only depends on experience and lacks corresponding theoretical guidance in actual production. It is necessary to determine the effect of cover thickness on soil evaporation and temperature to optimize the thickness of the sand layer. The objectives of this study were to determine (i) the effect of sand thickness on soil-wa-Vol. 20, No. 1, 2021 • Nature Environment and Pollution Technology ter content (SWC) and soil temperature during evaporation, and (ii) the optimum mulch thickness. This study provides a theoretical reference for the rational use of land and the transfer of soil water and heat for agricultural production in sand-mulched fields.

Soil-Atmosphere Boundary Conditions
The key to numerical simulation is to determine reasonable boundary conditions for the model. We chose the flux boundary for the soil surface and the flow-boundary conditions for quantifying surface infiltration and actual evaporation. Wilson et al. (1994) considered atmospheric humidity, wind speed, solar radiation and other conditions at the soil surface and derived an equation for calculating soil-surface evaporation under the conditions of soil-atmosphere coupling: were to determine (i) the effect of sand thickness on soil-water content (SWC) and soil temperature during evaporation, and (ii) the optimum mulch thickness. This study provides a theoretical reference for the rational use of land and the transfer of soil water and heat for agricultural production in sandmulched fields.

Soil-Atmosphere Boundary Conditions
The key to numerical simulation is to determine reasonable boundary conditions for the model.
We chose the flux boundary for the soil surface and the flow-boundary conditions for quantifying surface infiltration and actual evaporation. Wilson et al. (1994) considered atmospheric humidity, wind speed, solar radiation and other conditions at the soil surface and derived an equation for calculating soil-surface evaporation under the conditions of soil-atmosphere coupling: Where, AE is actual vertical evaporative flux (mm/day), Γ is the slope of the curve for saturated vapour pressure versus temperature at mean air temperature (kPa/°C), Qn is the total net radiation at the soil surface (mm/day), η is a psychrometric constant and Ea = f(u)ea(B-A), where f(u) = 0.35(1+0.15 Wa), Wa is wind speed (km/h), ea is the water-vapour pressure of the air above the soil surface (kPa), B is the inverse of atmospheric relative humidity and A is the inverse of the relative humidity at the soil surface.
Soil heat moves through the soil to the atmosphere, and atmospheric heat moves through the surface into the soil, so the temperature of the soil surface represents the temperature boundary condition and can be estimated (for conditions where no snowpack is present) as Wilson (1990): Where, Ts is the temperature of the soil surface (°C) and Ta is the temperature of the air above the soil surface (°C).

General Flow Law
VADOSE/W is formulated on the basis that the flow of water, heat, vapour and gas through both saturated and unsaturated soil follows an appropriate form of a Darcy-type flow law: Where, q is specific flux (m/s), k is conductivity (m/s) and i is a gradient of the potential.

…(1)
Where, AE is actual vertical evaporative flux (mm/day), Γ is the slope of the curve for saturated vapour pressure versus temperature at mean air temperature (kPa/°C), Q n is the total net radiation at the soil surface (mm/day), h is a psychrometric constant and E a = f(u)e a (B-A), where f(u) = 0.35(1+0.15 W a ), W a is wind speed (km/h), e a is the watervapour pressure of the air above the soil surface (kPa), B is the inverse of atmospheric relative humidity and A is the inverse of the relative humidity at the soil surface.
Soil heat moves through the soil to the atmosphere, and atmospheric heat moves through the surface into the soil, so the temperature of the soil surface represents the temperature boundary condition and can be estimated (for conditions where no snowpack is present) as Wilson (1990): were to determine (i) the effect of sand thickness on soil-water content (SWC) and soil temperature during evaporation, and (ii) the optimum mulch thickness. This study provides a theoretical reference for the rational use of land and the transfer of soil water and heat for agricultural production in sandmulched fields.

Soil-Atmosphere Boundary Conditions
The key to numerical simulation is to determine reasonable boundary conditions for the model.
We chose the flux boundary for the soil surface and the flow-boundary conditions for quantifying surface infiltration and actual evaporation. Wilson et al. (1994) considered atmospheric humidity, wind speed, solar radiation and other conditions at the soil surface and derived an equation for calculating soil-surface evaporation under the conditions of soil-atmosphere coupling: Where, AE is actual vertical evaporative flux (mm/day), Γ is the slope of the curve for saturated vapour pressure versus temperature at mean air temperature (kPa/°C), Qn is the total net radiation at the soil surface (mm/day), η is a psychrometric constant and Ea = f(u)ea(B-A), where f(u) = 0.35(1+0.15 Wa), Wa is wind speed (km/h), ea is the water-vapour pressure of the air above the soil surface (kPa), B is the inverse of atmospheric relative humidity and A is the inverse of the relative humidity at the soil surface.
Soil heat moves through the soil to the atmosphere, and atmospheric heat moves through the surface into the soil, so the temperature of the soil surface represents the temperature boundary condition and can be estimated (for conditions where no snowpack is present) as Wilson (1990): Where, Ts is the temperature of the soil surface (°C) and Ta is the temperature of the air above the soil surface (°C).

General Flow Law
VADOSE/W is formulated on the basis that the flow of water, heat, vapour and gas through both saturated and unsaturated soil follows an appropriate form of a Darcy-type flow law: Where, q is specific flux (m/s), k is conductivity (m/s) and i is a gradient of the potential.

…(2)
Where, T s is the temperature of the soil surface (°C) and T a is the temperature of the air above the soil surface (°C).

General Flow Law
VADOSE/W is formulated on the basis that the flow of water, heat, vapour and gas through both saturated and unsaturated soil follows an appropriate form of a Darcy-type flow law: Where, q is specific flux (m/s), k is conductivity (m/s) and i is a gradient of the potential.

Partial Differential Equations for Water and Heat Flow
The general governing differential equation for two-dimensional seepage can be expressed as:

Partial Differential Equations for Water and Heat Flow
The general governing differential equation for two-dimensi Where, P is pressure (kPa), Pv is the vapour pressure of conductivity in the x-direction (m/s), ky is hydraulic conductivity applied boundary flux (m/s), Dv is a coefficient for the diffusion m/(kN•s)), y is the elevation head (m), ρ is the density of wa acceleration (m/s 2 ) and t is time (s).
For the transfer of heat: , kty is thermal conductivity in the y-direction and is a is the Darcy water velocity in the x-direction (m/s), Vy is the Darc (m/s), Qt is the applied thermal-boundary flux (J/s) and Lv is the la

Coupling Heat and Mass Equations
Examination of the equations governing the transfer of heat parameters: pressure (P), temperature (T) and vapour pressure (P these parameters was necessary to solve the equations. This relat widely accepted thermodynamic relationship proposed by Edlefse Where, Pvs is the saturation vapour pressure (kPa) of so atmospheric relative humidity, Ψ is the total potential of the li equivalent matric potential (m), Wv is the molecular weight of wate gas constant (8.314 J/(mol•K)) and T is the temperature (K).

Geometric Model Partial Differential Equations for Water and Heat Flow
The general governing differential equation for two-dimensional seepage can be expressed as: Where, P is pressure (kPa), Pv is the vapour pressure of soil water (kPa), kx is hydraulic conductivity in the x-direction (m/s), ky is hydraulic conductivity in the y-direction (m/s), Q is the applied boundary flux (m/s), Dv is a coefficient for the diffusion of water vapour through soil (kg m/(kN•s)), y is the elevation head (m), ρ is the density of water (kg/m 3 ), g is the gravitational acceleration (m/s 2 ) and t is time (s).
For the transfer of heat: Where, ρc is the volumetric specific heat (J/(m 3 •°C), ktx is thermal conductivity in the x-direction (W/(m•°C)), kty is thermal conductivity in the y-direction and is assumed equal to ktx (W/(m•°C)), Vx is the Darcy water velocity in the x-direction (m/s), Vy is the Darcy water velocity in the y-direction (m/s), Qt is the applied thermal-boundary flux (J/s) and Lv is the latent heat of vaporization (J/kg).

Coupling Heat and Mass Equations
Examination of the equations governing the transfer of heat and mass identified three unknown parameters: pressure (P), temperature (T) and vapour pressure (Pv). The third relationship between these parameters was necessary to solve the equations. This relationship can be described using the widely accepted thermodynamic relationship proposed by Edlefsen & Anderson (1943): Where, Pvs is the saturation vapour pressure (kPa) of soil water at temperature T, hr is atmospheric relative humidity, Ψ is the total potential of the liquid-water phase expressed as an equivalent matric potential (m), Wv is the molecular weight of water (0.018 kg/mol), R is the universal gas constant (8.314 J/(mol•K)) and T is the temperature (K).

…(4)
Where, P is pressure (kPa), P v is the vapour pressure of soil water (kPa), k x is hydraulic conductivity in the x-direction (m/s), k y is hydraulic conductivity in the y-direction (m/s), Q is the applied boundary flux (m/s), D v is a coefficient for the diffusion of water vapour through soil (kg m/(kN•s)), y is the elevation head (m), r is the density of water (kg/m 3 ), g is the gravitational acceleration (m/s 2 ) and t is time (s).
For the transfer of heat:

Partial Differential Equations for Water and Heat Flow
The general governing differential equation for two-dimensional se Where, P is pressure (kPa), Pv is the vapour pressure of soil w conductivity in the x-direction (m/s), ky is hydraulic conductivity in the applied boundary flux (m/s), Dv is a coefficient for the diffusion of wa m/(kN•s)), y is the elevation head (m), ρ is the density of water (kg acceleration (m/s 2 ) and t is time (s).
For the transfer of heat: Where, Pvs is the saturation vapour pressure (kPa) of soil wat atmospheric relative humidity, Ψ is the total potential of the liquid-w equivalent matric potential (m), Wv is the molecular weight of water (0.01 gas constant (8.314 J/(mol•K)) and T is the temperature (K).

Partial Differential Equations for Water and Heat Flow
The general governing differential equation for two-dimensional seepage can be expressed as: Where, P is pressure (kPa), Pv is the vapour pressure of soil water (kPa), kx is hydraulic conductivity in the x-direction (m/s), ky is hydraulic conductivity in the y-direction (m/s), Q is the applied boundary flux (m/s), Dv is a coefficient for the diffusion of water vapour through soil (kg m/(kN•s)), y is the elevation head (m), ρ is the density of water (kg/m 3 ), g is the gravitational acceleration (m/s 2 ) and t is time (s).
For the transfer of heat: Where, ρc is the volumetric specific heat (J/(m 3 •°C), ktx is thermal conductivity in the x-direction (W/(m•°C)), kty is thermal conductivity in the y-direction and is assumed equal to ktx (W/(m•°C)), Vx is the Darcy water velocity in the x-direction (m/s), Vy is the Darcy water velocity in the y-direction (m/s), Qt is the applied thermal-boundary flux (J/s) and Lv is the latent heat of vaporization (J/kg).

Coupling Heat and Mass Equations
Examination of the equations governing the transfer of heat and mass identified three unknown parameters: pressure (P), temperature (T) and vapour pressure (Pv). The third relationship between these parameters was necessary to solve the equations. This relationship can be described using the widely accepted thermodynamic relationship proposed by Edlefsen & Anderson (1943): Where, Pvs is the saturation vapour pressure (kPa) of soil water at temperature T, hr is atmospheric relative humidity, Ψ is the total potential of the liquid-water phase expressed as an equivalent matric potential (m), Wv is the molecular weight of water (0.018 kg/mol), R is the universal gas constant (8.314 J/(mol•K)) and T is the temperature (K).

…(5)
Where, rc is the volumetric specific heat (J/(m 3 •°C), k tx is thermal conductivity in the x-direction (W/(m•°C)), k ty is thermal conductivity in the y-direction and is assumed equal to k tx (W/(m•°C)), V x is the Darcy water velocity in the x-direction (m/s), V y is the Darcy water velocity in the y-direction (m/s), Q t is the applied thermal-boundary flux (J/s) and L v is the latent heat of vaporization (J/kg).

Coupling Heat and Mass Equations
Examination of the equations governing the transfer of heat and mass identified three unknown parameters: pressure (P), temperature (T) and vapour pressure (P v ). The third relationship between these parameters was necessary to solve the equations. This relationship can be described using the widely accepted thermodynamic relationship proposed by Edlefsen & Anderson (1943):

Partial Differential Equations for Water and He
The general governing differential equation fo Where, P is pressure (kPa), Pv is the vapou conductivity in the x-direction (m/s), ky is hydraul applied boundary flux (m/s), Dv is a coefficient fo m/(kN•s)), y is the elevation head (m), ρ is the acceleration (m/s 2 ) and t is time (s).
For the transfer of heat:

Coupling Heat and Mass Equations
Examination of the equations governing the tr parameters: pressure (P), temperature (T) and vapo these parameters was necessary to solve the equati widely accepted thermodynamic relationship propo Where, Pvs is the saturation vapour pressur atmospheric relative humidity, Ψ is the total pote equivalent matric potential (m), Wv is the molecular gas constant (8.314 J/(mol•K)) and T is the tempera

…(6)
Where, P vs is the saturation vapour pressure (kPa) of soil water at temperature T, h r is atmospheric relative humidity, Ψ is the total potential of the liquid-water phase expressed as an equivalent matric potential (m), W v is the molecular weight of water (0.018 kg/mol), R is the universal gas constant (8.314 J/(mol•K)) and T is the temperature (K).

Geometric Model
The experiment simulating evaporation used an indoor test soil column with bare soil as the control group (CK) and columns with the soil surface mulched with sand 1.0, 3.0, 5.0, 6.0, 8.0, 10.0, 15.0 and 20.0 cm thick, simplified as a two-dimensional numerical model (Fig. 1). The diameter of the soil column was 10 cm, the height of the bare soil was 35.0 cm. The heights of the sand-mulched treatments were thus 36.0, 38.0, 40.0, 41.0, 43.0, 45.0, 50.0 and 55.0 cm. VADOSE/W was used to model the soil columns.

Model Parameters
Clay loam and sand were used in the experiments. The initial state of the soil column and the physical properties of the soil and sand are presented in Table 1.

Simulated Experimental Conditions
The soil-atmosphere boundary was the upper boundary of the model, which can be used for the exchange of water and heat, and the lower boundary was impermeable. The recharge of groundwater was not included in the model. The relative humidity was 10% during soil-column evaporation, the wind speed was 0 m/s and potential evaporation was 10 mm/d. Each mulching treatment was also tested with two temperature treatments, an endothermic process (ENP) and an exothermic process (EXP), for studying the evaporation The experiment simulating evaporation used an indoor test soil column with bare soil as the control group (CK) and columns with the soil surface mulched with sand 1.0, 3.0, 5.0, 6.0, 8.0, 10.0, 15.0 and 20.0 cm thick, simplified as a two-dimensional numerical model (Fig. 1)

Model Parameters
Clay loam and sand were used in the experiments. The initial state of the soil column and the physical properties of the soil and sand are presented in Table 1.

Simulated Experimental Conditions
The soil-atmosphere boundary was the upper boundary of the model, which can be used for the exchange of water and heat, and the lower boundary was impermeable. The recharge of groundwater was not included in the model. The relative humidity was 10% during soil-column evaporation, the wind speed was 0 m/s and potential evaporation was 10 mm/d. Each mulching treatment was also tested with two temperature treatments, an endothermic process (ENP) and an exothermic process (EXP), for studying the evaporation of sand-mulched soil at different temperatures. ENP had sand    (i = 0,1,3,5,6,8,10,15,20) represent the endothermic process and exothermic process of CK and the soil of sand thickness with 1,3,5,6,8,10,15,20 cm, respectively. of sand-mulched soil at different temperatures. ENP had sand surface temperatures of 30°C and soil internal temperature of 15°C, and EXP had sand surface temperatures of 15°C and soil internal temperature of 30°C (Table 2).

Effect of Sand Thickness on Cumulative Evaporation
The effect of sand thickness on cumulative evaporation for ENP and EXP are shown in Fig. 2. Cumulative evaporation for ENP and EXP was significantly higher in CK than the mulched treatments, regardless of sand thickness. This result was consistent with the findings of Wesemael et al. (1996), who also noted that sand mulching reduced evaporation. The reduction of cumulative evaporation decreased as sand thickness increased for all thickness treatments. The reduction was lowest for a thickness of 1 cm, but this thickness still considerably reduced evaporation compared to bare soil, because soil evaporation depends mainly on the capillary water movement (Xing et al. 2019). Mulching with sand can prevent capillary action in soil, so water can only diffuse as vapour, which greatly weakens the evaporation of soil water (Yamanaka et al. 2004).
Cumulative evaporation for CK and mulching with 1, 3, and 5 cm of sand were higher for ENP than EXP by 5.89, 2.1, 0.38 and 0.19 mm, respectively. Cumulative evaporation, however, was similar for ENP and EXP for sand thicknesses >5 cm. Cumulative evaporation decreased as sand thickness increased, consistent with the findings of Diaz et al. (2005), but did not vary significantly in our study for thicknesses >5 cm. Yamanaka et al. (2004) and Zhao et al. (2017b) also reported that resistance to evaporation did not increase with mulch thicknesses >5 cm. Wang et al. (2014), however, reported that evaporation was low when mulch thickness was >7 cm. The difference between these two conclusions may be due to different soil textures, particle-size distributions and study methods.

Effect of Sand Thickness on SWC
The effect of sand thickness on SWC for ENP and EXP is shown in Fig. 3. Mean SWC for ENP and EXP was lower in the surface layer (0-6 cm) than in other layers, especially for CK. Mean SWC for CK after evaporation was 0.13 m 3 /m 3 lower in the 0-6 cm than the 18-36 cm layer, and mean SWC in the sand-mulched treatments was lower in the 0-6 cm than the other layers, but the difference was small. Mean SWC was similar among the layers for sand thicknesses >5 cm.
Mean SWC at the end of evaporation was 35.5% lower than initial SWC for CK and 16.5, 4.9 and 2.2% lower for sand thicknesses of 1, 5 and 15 cm, respectively. Mean SWC for each layer was higher for the sand-mulched treatments than CK and was lower for ENP than EXP, especially for CK. These results indicate that sand mulching resists the transfer of water during evaporation, and the thicker the sand, the more the resistance, consistent with the conclusions by Li (2003) and Wang et al. (2004). Govers et al. (2010) and Modaihsh et al. (1985) also found that the degree of inhibition of evaporation and thus the reduction of water loss depended on the thickness of gravel mulch, with the most effective thicknesses of 5 cm (Govers et al. 2010) and 6 cm (Modaihsh et al. 1985), although Modaihsh et al. (1985) did not test a thickness of 5 cm.

Effect of Sand Thickness on Soil Temperature
The effect of sand thickness on soil temperature for ENP and EXP is shown in Fig. 4. Soil temperature for the treatments gradually increased with evaporation for ENP. Temperature after 15 days was highest for mulches 1 and 3 cm thick and then remained unchanged, while the other groups did not reach the highest temperature after evaporation. Heating was slower and more uniform for the mulched treatments than for CK. With the increase of the thickness of sand, the temperature transmission was slowed down. The change of soil temperature was not obvious at thicknesses >15 cm. Mean temperature after evaporation for ENP was higher for mulches 1, 3 and 5 cm thick than for CK. EXP and ENP were symmetrical, indicating that sand mulching preserved heat, consistent with the conclusions by Xie et al. (2010) and Lü et al. (2013).   Mean SWC at the end of evaporation was 35.5% lower than initial SWC for CK and 16.5, 4.9 and 2.2% lower for sand thicknesses of 1, 5 and 15 cm, respectively. Mean SWC for each layer was higher for the sand-mulched treatments than CK and was lower for ENP than EXP, especially for CK. These results indicate that sand mulching resists the transfer of water during evaporation, and the thicker the sand, the more the resistance, consistent with the conclusions by Li (2003)   Sand mulching could slow down the temperature transfer of ENP and EXP. Wang et al. (2014) also suggested that the thermal conduction effect of the mulch would lessen when the thickness of mulch is increased beyond a certain limit, which was 15 cm in our study.

CONCLUSION
Numerical simulation is to study the effects of sand thickness on cumulative evaporation, SWC and soil temperature under endothermic and exothermic soil conditions during evaporation.
(i) Sand mulching inhibited the movement of soil water during evaporation; the thicker the sand, the stronger the resistance. Cumulative evaporation was lower, and mean SWC was higher, for the mulched treatments than CK. Cumulative evaporation did not vary significantly for sand thicknesses >5 cm, and mean SWC for each layer varied little for ENP and EXP. The effect of sand thickness on soil temperature for ENP and EXP is shown in Fig. 4. Soil temperature for the treatments gradually increased with evaporation for ENP. Temperature after 15 days was highest for mulches 1 and 3 cm thick and then remained unchanged, while the other groups (ii) Sand mulching could slow down the temperature transfer. The changes in soil temperature for both ENP and EXP were slower and more uniform with sand mulching than for CK. The change of soil temperature was not obvious at thicknesses >15 cm. A sand thickness of 5 cm could thus be more effective for storing water and preserving heat.