Numerical Simulation of Soil Evaporation with Sand Mulching and Inclusion

Abstract

A model of unsaturated soil-water movement using a prediction model of basic physical soil properties for calculating correlation functions was developed using VADOSE/W. The reliability of the model was assessed by comparing the results with those of a soil-column test. Coefficients of determination, R², between the simulated and the measured daily evaporation for sand-mulch thicknesses of 0 (control, CK), 1.7, 3.6 and 5.7 cm were 0.8270, 0.8214, 0.8589 and 0.9851, respectively. R², between the simulated and measured cumulative evaporation for mulch thicknesses of 0, 1.7, 3.6 and 5.7 cm were 0.9755, 0.9994, 0.9997 and 0.9983, respectively. The fits were, thus, good, verifying the reliability of the model. The program accurately predicted the distribution of cumulative evaporation and volumetric water content during evaporation from a soil column with mulch thicknesses of 1, 1.3, 1.5, 1.7, 2, 3, 5 cm and depths of sand inclusion thick of 0, 5, 10 and 15 cm for 20 days. Cumulative evaporation of sand inclusion was lower than in CK. Cumulative evaporation was independent of the mulch thickness and depended only on the depth of the inclusion: the deeper the inclusion, the higher the evaporation. The best mulch thickness was 5 cm, and the best inclusion depth was 5 cm. This study offers a new method to study the evaporation process with sand mulching and inclusion, which can provide guidance for improving the utilization efficiency of soil water.

1. Introduction

Soil evaporation in arid and semiarid regions is an invalid loss of rainwater [1,2] and is an important component of evapotranspiration [3], which is the main method of soil water depletion in the field [4]. Soil evaporation is also a key factor that restricts agricultural planting [5,6]. Farmers have, thus, adopted sand mulching farming patterns on farmland in the arid and semi-arid regions of Northwestern China [7]. This ancient practice has been used in China for more than 300 years [8]. Evaporation involves many aspects, such as air temperature, solar radiation, relative humidity, wind speed, and so on [9]. Many studies have shown that adding sand mulch to the soil surface can effectively restrict the evaporation of soil water, storage, and preservation of soil moisture, and retain heat [10,11]. Wang [12] suggested that a method should be developed to quantify the dry soil layer. Soil texture in nature is affected by hydrology, meteorology, crustal movement, and human and biological activities. Most soil profiles have a layered structure rather than a single homogeneous layer. Some scholars have begun to study the effects of soil properties with adding sand inclusions at various depths [13,14].

Numerical simulation has recently been widely used in various fields. It can effectively control the variables, and ensure the accuracy of the test. In addition, it can deal with more complex problems and make up for the lack of experimentation. VADOSE/W (GEO-SLOPE International Ltd., Calgary, AB, Canada) is a finite-element software, based on a soil-atmosphere model focused mainly on the deformation of expansive soil, soil-water balance, and slope stability. Vanapalli et al. [15] found that expansive soil deformation simulated by VADOSE/W agreed well with the results of field studies. Benson et al. [16] found that water-balance data from a test section simulating a monolithic alternative cover were better for VADOSE/W than another numerical model, UNSAT-H. At present, the study of soil evaporation under sand mulching is focused on indoor experiments, and the numerical simulation test is lacking, so the software is applied in this field.

We used VADOSE/W software to simulate the indoor soil column evaporation experiment from bare soil, soil mulched with various thicknesses of sand, and sand inclusions at various depths. Cumulative soil evaporation and the distribution of volumetric water content were measured and simulated for sand-mulch thicknesses of 1, 1.3, 1.5, 1.7, 2, 3 and 5 cm and sand inclusions at depths of 0, 5, 10 and 15 cm for 0–20 days. The specific objectives are to (i) provide a theoretical basis for the movement of water through the calibrated appropriate thickness of the sand mulching and depth of the inclusion by evaporation with sand mulching and inclusion; (ii) provide a new research method for the evaporation of soil under the condition of sand mulching and inclusion.

2. Materials and Methods

2.1. Experimental Materials

Here, only a brief overview of the field experiment is provided. Further details of these experiment could be found in [17]. Clay loam and sand were used in the experiments. Soil samples were collected in the field from the top 50 cm, air dried (residual moisture content of 0.028 m³/m³), rolled, cleaned of debris, and sieved through a 2-mm mesh. Sand was air dried (residual moisture content of 0.017 m³/m³), cleaned of debris, and sieved through 0.2–2 mm meshes. The physical properties of the soil and sand are shown in

Table 1. Physical properties of the soil and sand in the experiment.

Soil Texture	Particle Size Distribution (%)						Field Capacity (m3/m3)	Residual Volume Moisture Content (m3/m3)	Hydraulic Conductivity (m/day)
	2.0–1.0 mm	1.0–0.5 mm	0.5–0.2 mm	0.2–0.02 mm	0.02–0.002 mm	<0.002 mm
Clay loam	0.07	0.38	2.54	36.99	39.16	20.86	0.3636	0.028	0.04
Sand	10.17	16.58	73.25	-	-	-	0.1532	0.017	17.856

.

Soil properties, such as the volumetric water-content function, the permeability coefficient function, can be predicted through indicators of basic physical soil properties by using VADOSE/W. The parameter-prediction method can be applied to the model based on the functions in the table, such as the moisture and the permeability coefficient of the product function.

2.2. Experimental Design

The experimental system included columns, a heating system, and a weighing system. The columns were made of Plexiglas and were 10 cm in diameter, 50 cm in height, and had seven drainage holes 1 cm in diameter at the bottom. The columns had pairs of 15-mm circular sampling ports at intervals of 5 cm vertically staggered at 90° for collecting samples from various depths for the analysis of soil moisture. Two layers of gauze were placed at the bottom of the columns to prevent the loss of soil particles through the drainage holes before the soil was loaded. The dried soil materials were packed into the columns to a height of 37.3 cm. The prepared sand was then loaded on the top at thicknesses of 0, 1.7, 3.6 and 5.7 cm, the initial volumetric water content was saturated, and the column was weighed before evaporation was initiated. All treatments had three replicates. The columns were heated by infrared light bulbs 30 cm above the soil surfaces at a constant intensity of illumination that maintained the temperature at 30 °C throughout the experiment. The columns were weighed during the evaporation by an electronic scale at 10:00 for 20 days for calculating daily evaporation. The potential evaporation of water from the surface was calculated using the cross-sectional areas of the columns. The potential evaporative capacity was 10 mm/day.

3. Model Formulation

3.1. Formulation of a Geometric Model

We used VADOSE/W for numerically simulating homogeneous bare soil as a column (control, CK) and for soil surfaces covered by a homogeneous horizon of sand, which simplified the modeling as a two-dimensional problem (Figure 1). The model control column was 10 cm in diameter partitioned into ten equal parts, and 37.3 cm in height partitioned into 38 equal parts, for a total of 380 units (Figure 1a). Three model sand-mulched columns with sand thicknesses of 1.7, 3.6 and 5.7 cm were 10 cm in diameter partitioned into ten equal parts and 39, 40.9 and 43 cm in height partitioned into 39, 42 and 44 equal parts, for totals of 390, 420 and 440 units, respectively (Figure 1b–d). The upper boundary of the model was the soil-atmosphere boundary, where water and heat can exchange, and the bottom boundary was impermeable. Recharge of groundwater is not considered in the model.

The above boundary conditions for water flow and temperature required that meteorological conditions be provided when the instantaneous solution, including temperature, humidity, and wind speed. The temperature was maintained at 30 °C, relative humidity was 10%, wind speed was 0 m/s and potential evaporation was 10 mm/day during soil-column evaporation. The soil column was initially saturated and was not supplemented with water during the test period of continuous evaporation for 20 days (these values of parameters are from [17]). The amount of evaporation was represented by the degree of the reduction in the quality of the soil column, and the reliability of the software was assessed by comparing simulated and measured daily and cumulative evaporation.

3.2. Geometric Models of Soil Columns with Different Sand-Mulch Thicknesses and Depths of Sand Inclusions

The experiment contained seven groups, each group has four treatments, for a total of 28 treatments. The seven groups were sand-mulch thicknesses of 1, 1.3, 1.5, 1.7, 2, 3 and 5 cm, and the four treatments were sand inclusions at depths of 0, 5, 10 and 15 cm from the upper surface of the soil column.

3.3. Boundary Conditions

The key for the numerical simulation is to determine reasonable boundary conditions for the model. We chose the flux boundary for the soil-surface element and the flow-boundary conditions for the quantification of surface infiltration and actual evaporation. The equation is as follows [18]:

$$AE=\frac{\mathsf{\Gamma }{Q}_{\mathrm{n}}+\eta E_{\mathrm{a}}}{\mathsf{\Gamma }+\eta A}$$

(1)

where AE is actual vertical evaporative flux (mm/day), Γ is the slope of the curve for saturation vapor pressure versus temperature at the mean air temperature (kPa/°C), Q_n is the total net radiation at the soil surface (mm/day), η is psychrometric constant and E_a = f(u)e_a(B−A), where f(u) = 0.35(1 + 0.15 W_a), W_a is the wind speed (km/h), e_a is the water-vapor pressure of the air above the soil surface (kPa), B is the inverse of the atmospheric relative humidity, and A is the inverse of the relative humidity at the soil surface.

Soil heat moves through the soil surface into the atmosphere, and atmospheric heat moves through the surface into the soil, thus, the temperature of the soil surface represents the boundary condition of the temperature and can be estimated (for conditions where no snow pack is present) with the following relationship [19]:

$$T_{\mathrm{s}}=T_{\mathrm{a}}+\frac{1}{\eta f(u)}(Q_{\mathrm{n}}-AE)$$

(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).

3.4. General Flow Law

VADOSE/W is formulated on the basis that the flow of water, heat, vapor, and gas through both saturated and unsaturated soil follows an appropriate form of a Darcy-type flow law which states that:

$$q=ki$$

(3)

where q is the specific flux (m/s), k is the conductivity (m/s), and i is gradient of potential.

3.5. Partial Differential Water and Heat Flow Equations

The general governing differential equation for two-dimensional seepage can be expressed as:

$${\lambda }_t\frac{\partial P}{\partial t}=\frac{1}{\rho }\frac{\partial }{\partial x}\left(D_v\frac{\partial P_v}{\partial x}\right)+\frac{1}{\rho }\frac{\partial }{\partial y}\left(D_v\frac{\partial P_v}{\partial y}\right)+\frac{\partial }{\partial x}\left(k_x\frac{\partial \left(P/\rho g+y\right)}{\partial x}\right)+\frac{\partial }{\partial y}\left(k_y\frac{\partial \left(P/\rho g+y\right)}{\partial y}\right)+Q$$

(4)

where P is the pressure (kPa), P_v is the vapor pressure of soil moisture(kPa), k_x is the hydraulic conductivity in the x-direction (m/s), k_y is the hydraulic conductivity in the y-direction (m/s), Q is the applied boundary flux (m/s), D_v is the diffusion coefficient of the water vapor through the soil (kg m/(kN s)), y is the elevation head (m), ρ is the density of water (kg/m³), g is the acceleration due to gravity (m/s²), and t is the time (s).

For heat transfer:

$${\lambda }_t\frac{\partial T}{\partial t}=L_v\frac{\partial }{\partial x}\left(D_v\frac{\partial P_v}{\partial x}\right)+L_v\frac{\partial }{\partial y}\left(D_v\frac{\partial P_v}{\partial y}\right)+\frac{\partial }{\partial x}\left(k_{tx}\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(k_{ty}\frac{\partial T}{\partial y}\right)+Q_t+\rho c{V}_x\frac{\partial T}{\partial x}+\rho c{V}_y\frac{\partial T}{\partial y}$$

(5)

where ρc is the volumetric specific heat value (J/(m³·°C), k_tx is the thermal conductivity in the x-direction (W/(m·°C)), k_ty is the thermal conductivity in the y-direction and 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).

3.6. Coupling Heat and Mass Equations

Examination of the governing heat and mass transfer equations reveals that there are three unknown parameters, namely pressure (P), temperature (T), and vapor pressure (P_v). In order to solve the equations, a third relationship between these parameter is necessary.

Their relationship can be described by using the widely accepted thermodynamic relationship given by [20]:

$$P_v=P_{vs}{h}_r=P_{vs}{e}^{\mathsf{\Psi }g{W}_v/RT}$$

(6)

where P_vs is the saturation vapor pressure (kPa) of the soil water at soil 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 temperature (K).

4. Results and Discussion

4.1. Model Validation

Daily and cumulative evaporation from the soil columns for the three thicknesses of sand mulch simulated by VADOSE/W agreed well with the measured results from the indoor experiment of the soil columns (Figure 2a,b). The correlation coefficients for the simulated and measured values of daily evaporation for CK and sand thicknesses of 1.7, 3.6 and 5.7 were 0.9094, 0.9063, 0.9268 and 0.9925, respectively. Simulated daily evaporation in CK was significantly lower than the measured value on day 3. The reason is that the dry soil layer is formed on the surface of the soil column, which significantly inhibits the evaporation of soil. This is consistent with [6], which studied the effects of different gravel-sand mulching degrees and size on soil moisture evaporation. It was also found that the evaporation of the soil surface was more intense and easy to form a dry soil layer. The coefficients of determination, R², for the simulated and measured daily evaporation for CK and sand thicknesses of 1.7, 3.6 and 5.7 were 0.8270, 0.8214, 0.8589 and 0.9851, respectively. The correlation coefficients for the simulated and measured cumulative evaporation for CK and sand thicknesses of 1.7, 3.6 and 5.7 were 0.9877, 0.9997, 0.9999 and 0.9991, respectively. R² for the simulated and measured cumulative evaporation for CK and 1.7, 3.6 and 5.7 were 0.9755, 0.9994, 0.9997 and 0.9983, respectively. The simulated and measured values of daily and cumulative evaporation, thus, agreed well for all soil columns, and the accuracy was high, indicating that VADOSE/W could be used to simulate soil conditions for improving the test accuracy.

The relationship between cumulative evaporation and time for the three thicknesses of sand mulch (

Table 2. Fitting equations of the relationship between the simulated and measured values of the cumulative evaporation and the time under sand mulching differing in thickness.

Sand Mulching (cm)	Fitting Equations		Determination Oefficient (R2)
CK	measured values	Y = 17.336ln(X) + 9.3423 (1 ≤ X ≤ 20)	0.9968
	simulated values	Y = 17.612ln(X) + 5.9751 (1 ≤ X ≤ 20)	0.9861
1.7	measured values	Y = 0.9182 X + 0.2321 (1 ≤ X ≤ 20)	0.9994
	simulated values	Y = 0.9429 X + 0.3334 (1 ≤ X ≤ 20)	1
3.6	measured values	Y = 0.2818 X + 0.6786 (1 ≤ X ≤ 20)	0.9999
	simulated values	Y = 0.2814 X + 0.8184 (1 ≤ X ≤ 20)	0.9999
5.7	measured values	Y = 0.185 X + 0.5465 (1 ≤ X ≤ 20)	1
	simulated values	Y = 0.1724 X + 1.1366 (1 ≤ X ≤ 20)	0.9953

) can be described by:

Y = Aln(X) + B

(7)

Y = KX + C

(8)

Equation (7) is a logarithmic model, where A and B are fitting parameters. Equation (8) is a linear model, where K is the slope, and C is a constant. X and Y in the two equations are the evaporation time (d) and the cumulative evaporation (mm), respectively. The relationship between cumulative evaporation and time in CK was logarithmic (fitted by Equation (7)). The relationships between the three thicknesses of sand mulch and time of cumulative evaporation were linear (fitted by Equation (8)). This regularity was consistent with Yuan et al. [21]. This phenomenon is mainly due to the inhibition of soil water evaporation with sand mulching, making the evaporation rate is greatly reduced. The slope (K) of the fitting equations for the sand-mulched treatments decreased as the thickness of the sand increased: K_5.7 < K_3.6 < K_1.7. R² was high (0.9861–1) for the fitting equations, and the fitting precision was high.

4.2. Model Application

4.2.1. Effects of Sand Mulching and Inclusion on Cumulative Evaporation

The effects of sand mulching and inclusion on soil cumulative evaporation for the various thicknesses and depths are shown in Figure 3. The cumulative evaporation under the simulated conditions of steady-state potential evaporation within the 20 days of evaporation for the sand-inclusion treatments was lower than CK. For example, cumulative evaporation for the sand-mulch thickness is 1 cm, and sand inclusions at depths of 5, 10, and 15 cm were 23.0, 32.4 and 40.51 mm lower, respectively, than the cumulative evaporation for CK of 60.80 mm. Cumulative evaporation for the inclusion depths followed the order 5 cm < 10 cm < 15 cm < CK, and the cumulative evaporation for each depth was the same in each group, indicating that the evaporation was mainly from the soil above the inclusion layer. The cumulative evaporation was independent of the thickness of sand-inclusion layer and depended only on the depth of the inclusions: the deeper the inclusion, the higher the cumulative evaporation. Cumulative evaporation for the sand-mulch thickness is 1, 1.3, 1.5, 1.7, 2, 3 and 5 cm, and sand inclusion at a depth of 0 cm (surface sand) was 55.15, 27.94, 22.91, 19.19, 17.07, 8.16 and 4.76 mm, indicating that the thicker the surface sand, the lower the cumulative evaporation. This result was consistent with the conclusion obtained by Qiu et al. [22], who found that the inhibition of sand gravel on the evaporation is closely related to the thickness and particle size. This paper only considered the thickness and depth of the sand, and sand mulching particle size factors will be considered in a future study.

The fitting equations for the relationship between cumulative evaporation and time for the various sand-mulch thicknesses and inclusion depths are given in

Table 3. Fitting equations of the relationship between the cumulative evaporation and time under sand layers of different thicknesses and at different layers.

Sand Inclusion Treatment		Fitting Equations	Determination Coefficient (R2)	Sand Inclusion Treatment		Fitting Equations	Determination Coefficient (R2)
Sand Thickness	Sand Layer			Sand Thickness	Sand Layer
1 cm	0 cm	Y = 2.7503X + 0.1684	1	2 cm	0 cm	Y = 0.8339X + 0.3918	1
	5 cm	Y = 4.5307ln(X) + 8.6873	0.9899		5 cm	Y = 3.9044ln(X) + 9.3992	0.9987
	10 cm	Y = 7.2833ln(X) + 10.643	0.9986		10 cm	Y = 7.0774ln(X) + 10.672	0.9984
	15 cm	Y = 10.04ln(X) + 10.565	0.9992		15 cm	Y = 9.9055ln(X) + 10.653	0.9990
1.3 cm	0 cm	Y = 1.3842X + 0.2598	1	3 cm	0 cm	Y = 0.3772X + 0.611	1
	5 cm	Y = 4.3125ln(X) + 8.9165	0.9948		5 cm	Y = 3.6174ln(X) + 9.8299	0.9926
	10 cm	Y = 7.2352ln(X) + 10.656	0.9985		10 cm	Y = 7.197ln(X) + 10.719	0.9983
	15 cm	Y = 9.9755ln(X) + 10.597	0.9991		15 cm	Y = 9.8946ln(X) + 10.699	0.9989
1.5 cm	0 cm	Y = 1.1306X + 0.2976	1	5 cm	0 cm	Y = 0.1881X + 1.0383	0.9987
	5 cm	Y = 4.1757ln(X) + 9.0677	0.9971		5 cm	Y = 3.3917ln(X) + 10.288	0.9777
	10 cm	Y = 7.217ln(X) + 10.661	0.9985		10 cm	Y = 7.2355ln(X) + 10.802	0.9980
	15 cm	Y = 9.9481ln(X) + 10.618	0.9991		15 cm	Y = 9.9159ln(X) + 10.800	0.9987
1.7 cm	0 cm	Y = 0.9429X + 0.3334	1	Note: 1 ≤ X ≤ 20
	5 cm	Y = 4.0563ln(X) + 9.21	0.9984
	10 cm	Y = 7.2028ln(X) + 10.669	0.9985
	15 cm	Y = 9.9274ln(X) + 10.631	0.9991

. Cumulative evaporation and time was best described by Equation (8) for an inclusion depth of 0 cm and by Equation (7) for the other groups. A and B in the fitting equations for the same inclusion depth in each group were similar, further indicating that the cumulative evaporation was dependent on the soil thickness above the inclusion layer, independent of the thickness of the sand mulch.

4.2.2. Effects of Sand Mulching and Inclusion on Volumetric Water Content

The changes of volumetric water content at the end of the evaporation for the various thicknesses and depths are shown in Figure 4. The volumetric water contents were similar for each group, and the means for the sand-mulch thickness at 1, 3 and 5 cm, and sand inclusion at a depth of 0 cm were 0.243, 0.319, 0.339 and 0.352 m³/m³, respectively. Volumetric water content for the inclusion at 0 cm was consistent with the surface sand mulch. Volumetric water content increased with the thickness of the sand mulch. The volumetric water content of each group of inclusion treatment meets the same regularity, the volumetric water content was a transition in the location of the sand inclusion layer. Volumetric water content above the sand inclusion was significantly lower than below, indicating that the inclusion layer could effectively inhibit evaporation from below the inclusion. The volumetric water content below 5 cm differed little from the content below 10 and 15 cm, indicating that an inclusion depth deeper than 5 cm had little effect on the volumetric water content. The best depth of inclusion was, thus, 5 cm. The effect on the volumetric water content below the inclusion was also very small when the thickness of the sand mulch was 5 cm. The best thickness of sand mulch was, thus, 5 cm. This conclusion is the same as that of Govers et al. [23] and Modaihsh et al. [24], who found that the different thickness of gravel in the reduction of water loss, and the degree of inhibition of evaporation is different, of which 5 cm [23] or 6 cm [24] thickness sand mulching to reduce the loss of water is the most effective. The difference between the two conclusions is because Modaihsh et al. [24] just set the sand thickness of 0, 2 and 6 cm, and the lack of thickness at 5 cm to study. In addition, the deficiency of this study is mainly to verify the factors of gravel particle size, and the validation of other factors needs to be further studied.

5. Conclusions

We used an indoor soil-column test to verify the reliability of numerical simulation based on the VADOSE/W finite element. The correlation coefficients for the simulated and measured daily evaporation for sand-mulch thicknesses of 0, 1.7, 3.6 and 5.7 were 0.9094, 0.9063, 0.9268 and 0.9925, with R² of 0.8270, 0.8214, 0.8589 and 0.9851, respectively. The correlation coefficients for the simulated and measured cumulative evaporation for mulch thicknesses of 0, 1.7, 3.6 and 5.7 were 0.9877, 0.9997, 0.9999 and 0.9991, with R² of 0.9755, 0.9994, 0.9997 and 0.9983, respectively, indicating that the numerical simulation using VADOSE/W was highly reliable.

The cumulative evaporation and volumetric water content were simulated for sand mulches and inclusions of various thicknesses and depths. Cumulative evaporation was not correlated with the thickness of the mulch, only with the depth of the inclusion, increasing with depth. The sand inclusion effectively inhibited evaporation from below it. The best mulch thickness was 5 cm, and the best inclusion depth was 5 cm. The model and program used in this study, generally, accurately represented the transport of water in the tested soil for the various mulch thicknesses and inclusion depths.