Simulation of Water Use Dynamics by Salix Bush in a Semiarid Shallow Groundwater Area of the Chinese Erdos Plateau

This study analyzed the water use of the Salix psammophila bush in a semi-arid area in northwest China using a Hydrus-1D model. The model incorporated the effect of thermally driven water flow coupling liquid water, water vapor and heat transport. The model was calibrated and validated using hourly field measurements of soil water content and temperature at different depths for a growing season of 154 days. Furthermore, another Hydrus-1D model was established to simulate environments with decreased heat, rainfall or temperature and an increased leaf area index using calibrated and validated parameters. Our results show that upward and downward thermally driven water vapor fluxes account for 0.11% and 0.28%, respectively, of the corresponding direction of total water flux during the bush’s growing season. Although the vapor flux is very small, simulations incorporating heat flow revealed alterations in the temperature and pressure head gradients over the root zone, especially during dry periods. Consequently, the cumulative contributions of groundwater to evapotranspiration (ETg) with heat flow and without heat flow were 26.9% and 40.6%, respectively, during the simulation period. Therefore, the cumulative contribution of groundwater to ETg is overestimated when heat flow is excluded. Thus, we recommended that heat transport be incorporated when evaluating ETg in arid and semi-arid areas.


Introduction
Arid and semi-arid regions compose over half of the global land area [1].These areas often suffer from frequent sandstorms, and farmland is threatened by moving sands, particularly in China [2].In the early 1980s, the Chinese government initiated a reforestation project called the "Three North Forest Shelterbelts" [3], and in 2000, the "Return Farmland to Forest and Grassland" project was initiated [4].These reforestation projects have prevented further deterioration of arid regions and have rehabilitated the eco-environment.However, recent research has revealed side effects stemming from the improper selection of plants [5], such as the decay of re-vegetation and a decrease in Water 2015, 7, 6999-7021; doi:10.3390/w7126671www.mdpi.com/journal/water Water 2015, 7, 6999-7021 groundwater levels [6].In arid areas, water resources are extremely scarce and the environment is very fragile.Surface water is usually limited, and groundwater sustains the co-existence of social-economic development and the natural environment.Thus, understanding the relationship between water use by vegetation and groundwater dynamics is crucial to managing and maintaining healthy ecosystems while providing water for human needs [7].
In desert regions, rainfall events can be characterized as rainfall pulses with a discontinuous, highly variable, and largely unpredictable frequency and intensity [8,9].This small amount of rainfall only recharges the shallow soil layer; in contrast, heavy rainfall can permeate into the deep soil layer and groundwater [2].Many perennial plants in arid and semiarid zones have adapted to these sporadic rainfall pulses by adopting a dimorphic root system composed of branched surface roots that acquire water from shallow soil layers and deeper roots that can access perennially available groundwater [10,11].Plants able to extract groundwater from shallow aquifers are commonly referred to as phreatophytes.Some phreatophytes appear to be largely dependent on groundwater, while others show only a slight dependency.Thus, phreatophytes exhibit an obligate and facultative groundwater dependency [12].However, the identification of obligate and facultative phreatophytes is complicated by species characteristics [13,14] and shifts between wetter and drier soil conditions [11].
Temperature-driven soil water flow is another process that can affect water balance and movement in arid and semi-arid areas [15].The dynamics of liquid water and water vapor movement have been determined based on the mathematical description of liquid water and vapor flux by Philip and de Vries [16] for unsaturated soils under a soil water pressure head and a soil temperature gradient.This study produced a mathematical model [1,17] based on the influence of air flow [18], meteorological conditions [19], and vegetation cover [15,20].However, there have been limited modeling studies evaluating the role of coupled fluxes and their temporal variation in desert bush-dominated areas under shallow groundwater level conditions.In particular, information about the movement of liquid water, water vapor, and heat in bush-dominated sandy soil as well as the contribution of water vapor flux to total water flux remains limited for areas such as the Mu Us Desert.
The objective of this study was to assess the water use of Salix psammophila (S. psammophila) in a semi-arid area of the Mu Us Desert of northwest China using a coupled liquid water, water vapor, and heat transport model.Systematic field measurements were conducted to investigate the relationships between S. psammophila transpiration and groundwater level changes.First, we calibrated the HYDRUS-1D model [21] using soil water and temperature measurements obtained in an S. psammophila bush plot.We then used the model to compute the contributions of groundwater to root water uptake to quantify the dependency of the S. psammophila bush on groundwater.Third, we investigated the movement of liquid water and water vapor in the absence of heat transport and under decreased rainfall, increased temperature and LAI.

Field Site Description
The experimental site was located in the Hailiutu River catchment (between 38 ˝06 1 and 38 ˝52 1 N, 108 ˝47 1 and 109 ˝23 1 E) in the Erdos Plateau of northwest China (Figure 1).The Hailiutu River is a branch of the Wuding River, which is a major tributary of the Yellow River.The total area of the Hailiutu River catchment is approximately 2645 km 2 .A hydrological station located at the outlet of the Hailiutu catchment near Hanjiamao village registered a mean annual discharge of 2.64 m 3 /s for the period 1957-2007.Studies show that the flow regime has changed dramatically over the last 51 years.Four major shifts were detected in 1968, 1986, 1992 and 2001, reflecting quasi-natural conditions, reservoir construction, combined river diversion and groundwater extraction, and increased crop area [22].Based on 50 years of meteorological data (from 1957-2007) obtained from a station approximately 40 km northwest of the study site, the long-term average annual air Water 2015, 7, 6999-7021 temperature is 8.1 ˝C, and the minimum and maximum monthly average air temperatures are ´8.6 ˝C in January and 23.9 ˝C in July.The average annual precipitation is 340 mm, approximately 70% of which occurs between July and September.The average annual potential evaporation is estimated to be 2180 mm [23].The area's main crop is maize, which is sown in mid-April and harvested in mid-October, with slight variation depending on seasonal weather conditions.Irrigation depends on local rainfall [23].
The experimental site is located approximately 11 km northeast of the Hailutu River.The dominant natural vegetation around the research site is S. psammophila, which covers approximately 30% of the area, along with other sporadically distributed herbs, primarily Artemisia ordosica, Korshinsk peashrub and Hedysarum laeve Maxim.
Water 2015, 7 4 Figure 1.Location of the Hailiutu River catchment and the research site.

Field Measurements
A set of instruments was installed to measure environmental variables, including the sap flow of S. psammophila, soil water, soil temperature, and groundwater level (Figure 2).The instruments used in the study are listed in Table 1.Detailed information regarding river catchment and groundwater level measurements can be found in a previous study [23].Here, the results of a root survey and measurements of soil water content and soil temperature are described.

Field Measurements
A set of instruments was installed to measure environmental variables, including the sap flow of S. psammophila, soil water, soil temperature, and groundwater level (Figure 2).The instruments used in the study are listed in Table 1.Detailed information regarding river catchment and groundwater level measurements can be found in a previous study [23].Here, the results of a root survey and measurements of soil water content and soil temperature are described.Soil water content (SWC, cm 3 /cm 3 ) was measured at eight depths using time-domain reflectometry (TDR) technology.All sensors were inserted horizontally into the soil.All sensors Water 2015, 7, 6999-7021 installed at 10, 20, 40, 60, 80, 100, 120 and 140 cm depths were 6005CL2 (Minitrase SEC Co. Ltd., Goleta, CA, USA, 2% resolution).Data were recorded every 10 s and stored every 1 hour by a data-logger (Minitrase SEC Co. Ltd., Goleta, CA, USA).

Measurement of Soil Temperature
The soil temperature was measured at 10 different depths using TCAV thermocouple sensors (Campbell Scientific, Logan, UT, USA).The thermocouple sensors were installed at 2, 5, 10, 20, 30, 40, 50, 60, 70 and 80 cm depths.All sensors were inserted horizontally into the soil.The data were recordedat 10-s intervals and were stored as 1-h averages using a CR1000 data-logger (Campbell Scientific, Logan, UT, USA).

Root Survey
Root distribution reflects the location of the water source used by the plants.Because root distribution varies with depth and radial distance, soil cores were collected around S. psammophila plants using a root auger (Eijkelkamp, Giesbeek, The Netherlands) after defoliation.Undisturbed uniform soil samples 10 cm in diameter and 15 cm in height were obtained, and the position of each sample was recorded at the horizontal radial distance from the stem and the depth to the point of each sample.Four symmetrical transects centered at the stem with a radius of 4 m were investigated.The total area of the root distribution was determined to be 50.24m 2 .Roots were sieved out of the core sample and washed on the day of sampling.The cleaned roots from each core sample were weighed using an electric balance (readout 0.01 g; resolution limit of 1%), and photos were taken on scaled paper (resolution: 1 mm).The total root length of fine roots (diameter < 2 mm) was measured using gvSIG (a geographic information system ) open source software [24].The root length density (cm/cm 3 ) of each sample was determined by dividing the total root length by the core volume.The vertical distribution of fine roots was determined from the average of the root-length density at the same sample depth and the radial distance within the root zone of a S. psammophila plant.Two maxima of root length density were observed: one between 0 and 30 cm and another between the 60 and 120 cm (with a peak at approximately 80-90 cm) soil layers.

Simulations of Soil Water and Heat Transport
HYDRUS-1D software [21] was used to simulate water flow, water vapor and heat transport under root water uptake conditions.The model simulation followed a 3-step approach similar to that previously published [25].The first step involves the calibration of soil parameters using hourly measurements from 26 April to 26 May 2012.The second step is the validation of the model using hourly measurements from 27 May to 27 September 2012.The third step is a complete simulation from April to September 2012 to investigate the groundwater dependency of the Salix bush.

Governing Equation
The governing equation for the one-dimensional vertical flow of liquid water and water vapor in variably saturated media is given by the following mass conservation equation [17]: where θ T is the total volumetric water content (cm 3 /cm 3 ), which is the sum of the volumetric liquid and vapor water content (θ T = θ L + θ v ); q L and q v are the flux densities of liquid water and water vapor (cm/h), respectively; t is time (h); z is the vertical axis position (cm); and S(h) is the sink term (cm 3 /cm 3 /h).
Water 2015, 7, 6999-7021 The flux density of liquid water, q L , has been defined by Philip, J.R., et al. [16]: where q Lh and q LT are the isothermal and thermal liquid water flux densities (cm/h); h is the matrix potential head (cm); T is the temperature ( ˝C); and K Lh (cm/h) and K LT (cm 2 / ˝C/h) are the isothermal and thermal hydraulic conductivities for liquid-phase fluxes due to gradients in h and T, respectively.Using the product rule for differentiation and assuming that the relative humidity in soil pores is constant at different temperatures, the flux density of water vapor, q v , can be written as [17]: where q vh and q vT are, respectively, the isothermal and thermal water vapor flux densities (cm/h); and K vh (cm/h) and K vT (cm 2 / ˝C/h) are the isothermal and thermal vapor hydraulic conductivities of water vapor, respectively.The total heat flux density is defined as the sum of the conduction of sensible heat, as described by Fourier's law (the first term on the right side), the convection of sensible heat by liquid water (the second term) and water vapor (the third term), and the convection of latent heat by vapor flow (the fourth term) [17]: where λ is the apparent thermal conductivity of the soil (J/cm/h/ ˝C); C p , C w , and C v are the volumetric heat capacities (J/cm 3 / ˝C) of the soil, liquid, and vapor phases, respectively; q L is the sum of the isothermal and thermal liquid water flux densities (cm/h); L 0 is the latent heat of vaporization of liquid water (J/cm 3 ); and q v is the sum of the isothermal and thermal vapor flux densities (cm/h).Details regarding the HYDRUS-1D model can be found in previously published work [17,21].The sink term in Equation (1), S(h), is defined as the volume of water removed from a unit volume of soil per unit time due to plant water uptake.Feddes, R.A., et al. [26] defined S(h) as: where the root-water uptake reduction factor α(h) is a prescribed dimensionless function of the soil water pressure head, and b(x) is a non-uniform distribution of the potential water uptake rate over a root zone normalized by the field root survey value determined by Huang, J.T. et al. [23].T p is the potential transpiration.

Soil Characteristics
The water retention curve is one of the most fundamental hydraulic characteristics in the soil water flow equation.The soil water retention equation is given by van Genuchten [27]: where θ is the volumetric water content (cm 3 /cm 3 ) at pressure head h (cm); θr and θs are the residual and saturated water contents (cm 3 /cm 3 ), respectively; α (> 0, in 1/cm) is related to the inverse of the air-entry pressure; and n (>1) is a measure of effect of the pore-size distribution on the slope of the retention function (m = 1 ´1/n).
Water 2015, 7, 6999-7021 The hydraulic conductivity of soil K(h) is described as: where K s is the saturated hydraulic conductivity.S e is the effective saturation: The thermal conductivity was defined by de Marsily [28] as: where C w is the volumetric heat capacity of the water and q is water flux; β t is the thermal dispersivity (cm); and λ 0 is the baseline thermal conductivity, which has been described by Chung, S.O., et al. [29] as: where b 1 , b 2 and b 3 are empirically determined parameters (W/cm/ ˝C).

Root Water Uptake Modeling
In equation ( 5), the root distribution b(x) was normalized according to the field root survey data.For α(h), we used the functional formula introduced by van Genuchten [30]: where h 50 is the pressure head at which transpiration is halved and p is an adjustable constant that determines the steepness of the transition from potential to reduced uptake rates as h decreases.

Boundary and Initial Conditions
An atmospheric boundary condition was implemented at the soil surface, while a variable pressure head condition was used at depth.The evaporation from the soil surface and transpiration by plants were simulated using the HYDRUS-1D model.To specify the potential transpiration (T p ) and potential evaporation (E p ), we calculated the potential evapotranspiration (ET p ) using the Penman-Monteith equation [31].Potential evaporation (E p ) and potential transpiration were then calculated according to the Beer's law method: where LAI is the leaf area index [´] and k is an extinction coefficient set to 0.463.In our study, LAI was derived from MODIS imagery (MOD15A2) at a temporal resolution of 8 days, with LAI reaching a peak of ~1.69.Thus, the computed potential evaporation was used as an input to calculate the actual evaporation fluxes based on a reduction of van Genuchten's equation for transpiration and the h Crit limit for soil evaporation.In our simulations, h Crit was determined from the equilibrium conditions between soil water and atmospheric water vapor.During the experiments, the observed pressure heads at the research site were used for model calibration and validation.Figure 3 summarizes the imposed surface and bottom conditions, showing the hourly values of precipitation, potential evapotranspiration, and depth to groundwater table.
Water 2015, 7, 6999-7021 with LAI reaching a peak of ~1.69.Thus, the computed potential evaporation was used as an input to calculate the actual evaporation fluxes based on a reduction of van Genuchten's equation for transpiration and the hCrit limit for soil evaporation.In our simulations, hCrit was determined from the equilibrium conditions between soil water and atmospheric water vapor.During the experiments, the observed pressure heads at the research site were used for model calibration and validation.Figure 3 summarizes the imposed surface and bottom conditions, showing the hourly values of precipitation, potential evapotranspiration, and depth to groundwater table.The upper and lower boundary conditions for heat transport were specified as temperature boundary conditions.The upper boundary values were given as the observed temperatures at 0.02 cm, and the lower boundary values were given as the groundwater temperature.The initial soil water content and soil temperature were determined from the measured values on April 26 at 0:00 h by interpolating the measured values between different depths.

Criterion for Model Calibration
Model calibration was evaluated using root mean square error (RMSE), systematic RMSE (RMSE s ), unsystematic RMSE (RMSE u ), and index of agreement (d) [32,33], defined as follows: where the RMSE s indicates errors due to under-or over prediction by showing how far the data fluctuate from the 1:1 line in a scatter plot of computed versus measured values.When the RMSE u is minimized and the RMSE s approaches RMSE, the model performs with maximum accuracy.The value of d ranges from 0 (no agreement) to 1 (a perfect fit between simulated and measured values).
N is the number of paired observations, and P i and O i are the computed and measured values, respectively.O is the mean of measured values, and Pi is defined as: where a and b are the intercept and slope of the least-squares linear regression between P i and O i .

Investigation of Liquid Water and Water Vapor Movement
To quantify the effect of water movement on the water use of the Salix bush, the calibrated model was used to investigate liquid water and water vapor movement under two scenarios.The first scenario quantified the contribution of soil water to S. psammophila under conditions with or without heat using the measured data.The second scenario evaluated liquid water and water vapor movement in the case of decreased rainfall (i.e., half the amount of rainfall during the experimental period, in agreement with historical values during periods of low rainfall) while increasing temperature and LAI (5 ˝C temperature and 1.5 LAI increases [34] were assumed, respectively) using the calibrated model.

Model Calibration and Validation
HYDRUS-1D software was used to simulate water movement at the experimental site.The soil hydraulic parameters and soil heat parameters were calibrated and validated using the soil water content and temperature data while varying the drought stress parameters h 50 and p.During the model calibration period (from April 26 to July 11, 2012), in general, the soil water content gradually increased from the ground surface to the deep layers, unless rainfall occurred.Moreover, the soil water content in the shallow layers was more sensitive to rainfall than that in the deep layers.For example, a daily rainfall of 30.3 mm occurred on June 27 and triggered an increase in the soil water content at depths up to 100 cm, while after smaller rainfall events, the increase in the soil water content only reached 40 cm.During the validation period (July 12 to September 28, 2012), the soil water content increased with soil depth to groundwater.However, a large rainfall event of 41.1 mm on July 20 and subsequent rainfall resulted in an increase in the groundwater level, which thereafter remained high.The soil water content at a depth of 100 cm remained saturated until the conclusion of our experiments.
Six parameters in the van Genuchten model [27] were calibrated using field measurements (i.e., θ r , θ s , α, n, l, and Ks).Calibration was performed by fitting the observed and modeled soil water contents using the Marquardt-Levenberg optimization algorithm.HYDRUS-1D ran the optimization process until it found the highest R 2 values [21] between observed and computed soil water content.
The soil column was schematized as six soil layers based on our in situ investigations.The parameter Ks was determined using the inverse auger method [35] for each layer.The remaining parameters were estimated using Rosetta [36], a pedotransfer function that predicts hydraulic parameters from soil texture data (Table 2).Running Hydrus-1D with Rosetta hydraulic parameter estimates and empirically determined values for saturated hydraulic conductivity in simulations produced poor agreement for the index.We therefore attempted to calibrate the soil hydraulic property model using the measured soil water content.In Equation ( 6), there are five unknown parameters for each layer.When using the inverse model in Hydrus-1D software, we found that fitting all five parameters in six layers at the same time tended to cause the inverse algorithm to fail.Thus, the five parameters were fitted layer by layer.
Water 2015, 7, This fitting method has been reported previously by other researchers [37].When the simulated and measured soil water contents exhibited a good agreement index, the hydraulic parameters were fixed and the drought stress parameters h 50 and p were fitted via the measured sap flow of S. psammophila.Finally, the hydraulic and drought stress parameters were fixed, and the thermal parameters were fitted using the soil temperature.
The temporal variation of the computed soil water content at the eight measurement depths during the calibration and validation periods was reasonably consistent with the field measurements at each depth (Figure 4).We found that the parameters α and n were more sensitive than the other parameters.Table 3 shows the calibrated hydraulic parameters.The statistical criteria for model calibration (Equations ( 13)-( 16)) are summarized in Table 4.The RMSE values were very low for both the calibration and validation periods at all measured depths.The index of agreement (d) was very high, ranging from 0.65 to 0.96 and from 0.86 to 1.00, respectively, for the model calibration and validation periods at all depths.Generally, the calibration and validation results were acceptable.The soil water contents computed using the model captured the sharp increase in the soil water content after heavy rainfall events.Reported parameter values for root-water uptake reduction by specific plants and soils range from approximately ´1000 cm to ´5000 cm for h 50 and from 1.5 to 3 for p [38].However, especially coarse soil, such as sand, is almost completely drained of water at a fairly modest pressure head (i.e., ´300 or ´400 cm).Using values for h 50 and p that are similar to those reported in the literature caused the uptake reduction to perform essentially as a step function [38].On this basis, we employed an h 50 value that was considerably lower than that reported in the literature.Likewise, a larger value of p was required to account for the steepness of the soil water retention curve.Similar to the approach reported by Zhu, Y., et al. [39], we performed simulations using a range of values for h 50 and p in an effort to calibrate the HYDRUS-1D model.As expected for Aeolian sand, the simulated water contents were not very sensitive to h 50 and p.This finding agrees with the results of Zhu, Y., et al. [39], which simulated Populus euphratica root uptake in coarse sand soil.A comparison between simulated and measured soil water contents is presented in Figure 4.For our research plants, the values of h 50 and p were found to be ´630 cm and 3, respectively, from the results of fitting the predicted and measured values for transpiration (Figure 5).It can be seen that the observed and simulated transpiration of S. psammophila fitted well (with R 2 = 0.99), indicating that the calibration values of h 50 and p were acceptable.
Hourly measurements of soil temperature at nine different depths reveal clear diurnal fluctuations (Figure 6).However, both the temperature and the amplitude of the diurnal fluctuations decrease with the increase of depth because the air temperature is higher than the groundwater temperature during the measurement period spanning from 26 April to 27 September 2011, indicating downward heat transport.The temperature also exhibited seasonal variation, increasing beginning in late April, reaching the highest temperature in July and August and decreasing in September.The calibrated HYDRUS-1D model is able to simulate both seasonal and diurnal variations in soil temperature at nine different depths (Figure 6).The statistical measures used for temperature calibration and validation are shown in Table 5.The RMSE ranges from 0.446 ˝C to 2.396 ˝C and Water 2015, 7, 6999-7021 from 0.784 ˝C to 2.175 ˝C during the calibration and validation periods, respectively.The index of agreement (d) is larger than 0.88, during both the calibration period and the validation period, indicating a good agreement between the computed and measured soil temperatures (Table 5).
account for the steepness of the soil water retention curve.Similar to the approach reported by Zhu, Y., et al. [39], we performed simulations using a range of values for h50 and p in an effort to calibrate the HYDRUS-1D model.As expected for Aeolian sand, the simulated water contents were not very sensitive to h50 and p.This finding agrees with the results of Zhu, Y., et al. [39], which simulated Populus euphratica root uptake in coarse sand soil.A comparison between simulated and measured soil water contents is presented in Figure 4.For our research plants, the values of h50 and p were found to be −630 cm and 3, respectively, from the results of fitting the predicted and measured values for transpiration (Figure 5).It can be seen that the observed and simulated transpiration of S. psammophila fitted well (with R 2 = 0.99), indicating that the calibration values of h50 and p were acceptable.

Pressure Head over the Root Zone
The mean values of the pressure head over the root zone with and without heat flow are shown in Figure 7.During the simulated period, the mean value of the pressure head ranged from ´9408.9 to ´31.11 cm with heat flow and from ´244.3 to ´32.5 cm without heat flow.In the simulation without heat, liquid connectivity breaks down (via Equation 7) and the water content cannot decrease further, whereas in the simulation with heat transport, very small amounts of water continue to be transported.These small changes have a large effect on the simulated matrix potential because the simulation approaches the residual water content.The lowest value of the pressure head occurred mainly during periods of no precipitation and a deep groundwater table (Figure 3).

Pressure Head over the Root Zone
The mean values of the pressure head over the root zone with and without heat flow are shown in Figure 7.During the simulated period, the mean value of the pressure head ranged from −9408.9 to −31.11 cm with heat flow and from −244.3 to −32.5 cm without heat flow.In the simulation without heat, liquid connectivity breaks down (via Equation 7) and the water content cannot decrease further, whereas in the simulation with heat transport, very small amounts of water continue to be transported.These small changes have a large effect on the simulated matrix potential because the simulation approaches the residual water content.The lowest value of the pressure head occurred mainly during periods of no precipitation and a deep groundwater table (Figure 3).

Evapotranspiration
The calculated evapotranspiration is shown in Figure 8, both with and without heat flow.During the simulated period, the total transpiration rates were 27.9 cm and 28.1 cm with and without heat flow, respectively.The total evaporation rates were 15.9 cm and 16.7 cm with and without heat flow, respectively.The differences were 0.2 cm and 0.8 cm for transpiration and evaporation, respectively (Figure 8a).In total, the simulated value of ET with heat flow is lower than that without heat flow; in other words, the value of ET is overestimated without heat flow.
To further clarify the effect of heat flow on transpiration, root water uptake values for two zones are shown in Figure 8b.It can be seen that the root water uptake is affected by heat flow mainly from 0 cm to 30 cm, and the value was 10.4 cm and 10.6 cm with and without heat flow, respectively.The value of root water uptake in the zone from 60 cm to 90 cm is the same-6.1 cm-both with and without heat flow.

Evapotranspiration
The calculated evapotranspiration is shown in Figure 8, both with and without heat flow.During the simulated period, the total transpiration rates were 27.9 cm and 28.1 cm with and without heat flow, respectively.The total evaporation rates were 15.9 cm and 16.7 cm with and without heat flow, respectively.The differences were 0.2 cm and 0.8 cm for transpiration and evaporation, respectively (Figure 8a).In total, the simulated value of ET with heat flow is lower than that without heat flow; in other words, the value of ET is overestimated without heat flow.
To further clarify the effect of heat flow on transpiration, root water uptake values for two zones are shown in Figure 8b.It can be seen that the root water uptake is affected by heat flow mainly from 0 cm to 30 cm, and the value was 10.4 cm and 10.6 cm with and without heat flow, respectively.The value of root water uptake in the zone from 60 cm to 90 cm is the same-6.1 cm-both with and without heat flow.

Water Storage and Bottom Flux
The change in soil water storage and bottom flux during the simulated period is shown in Figure 9.The cumulative changes in soil water storage were 11.7 cm and 11.4 cm with and without heat flow, respectively.The heat flow thus has little influence on the change in soil water storage.
During the simulated period, the cumulative bottom flux was positive, indicating the net groundwater inflow to the soil column.The total bottom fluxes were 14.4 cm and 15.1 cm with heat flow and without heat flow, respectively.Furthermore, groundwater inflow occurs during dry days, indicating the dependency of this process on groundwater.During heavy rain, the cumulative bottom fluxes decreased, indicating groundwater recharge.For example, after 4.87 cm rain fall occurred on the June 27 and June 28, value of cumulative bottom fluxes decreased from 11.18 cm to 9.55 cm without heat and from 9.42 cm to 7.89 cm with heat, respectively.This indicated value of recharge was 1.63 cm without heat and 1.53 cm with heat, respectively.

Water Storage and Bottom Flux
The change in soil water storage and bottom flux during the simulated period is shown in Figure 9.The cumulative changes in soil water storage were 11.7 cm and 11.4 cm with and without heat flow, respectively.The heat flow thus has little influence on the change in soil water storage.
During the simulated period, the cumulative bottom flux was positive, indicating the net groundwater inflow to the soil column.The total bottom fluxes were 14.4 cm and 15.1 cm with heat flow and without heat flow, respectively.Furthermore, groundwater inflow occurs during dry days, indicating the dependency of this process on groundwater.During heavy rain, the cumulative bottom fluxes decreased, indicating groundwater recharge.For example, after 4.87 cm rain fall occurred on the June 27 and June 28, value of cumulative bottom fluxes decreased from 11.18 cm to 9.55 cm without heat and from 9.42 cm to 7.89 cm with heat, respectively.This indicated value of recharge was 1.63 cm without heat and 1.53 cm with heat, respectively.Water 2015, 7, 6999-7021

Daily Movement of Liquid Water and Vapor Water in Summer
Figure 10a shows the liquid water movement temporal-space character before and after a heavy rainfall event (48.7 mm, occurred from 11:00 27 June to 4:00 29 June).Before the rainfall, the liquid water moved upward because of evapotranspiration, and the value ranged from 0 to 0.038 cm/h during the daytime; during the same period, transpiration ceased, a zero gradient zone formed at a depth of approximately 80 cm, and the liquid water below gradient zone moved downward with a value ranging from 0 cm/h to 0.08 cm/h.During the rainfall, the liquid water moved downward, and the maximum water flux reached 0.69 cm/h.Following the rainfall, because of ET, the liquid water moved upward at the upper soil profile and moved downward continually over a 54-h period.At the end of the infiltration process, the change in the liquid water showed the same characteristics as before the rainfall.
Figure 10b shows the water vapor characteristics in the temporal-space before and after the same heavy rainfall event.It can be seen that the rainfall event only influenced the shallow soil layer (around a depth of 0 cm to 40 cm) water vapor movement, and the direction was upward.During the other days, the movement direction of water vapor was different during the day (moved upward) and night (moved downward) at a depth of 0 to 40 cm.At a depth of 40 cm to the groundwater table, the water vapor moved downward.This phenomenon indicates that the vapor circulation-condensation process supplied additional moisture for root water uptake at the shallow soil profile.

Daily Movement of Liquid Water and Vapor Water in Summer
Figure 10a shows the liquid water movement temporal-space character before and after a heavy rainfall event (48.7 mm, occurred from 11:00 27 June to 4:00 29 June).Before the rainfall, the liquid water moved upward because of evapotranspiration, and the value ranged from 0 to 0.038 cm/h during the daytime; during the same period, transpiration ceased, a zero gradient zone formed at a depth of approximately 80 cm, and the liquid water below the gradient zone moved downward with a value ranging from 0 cm/h to 0.08 cm/h.During the rainfall, the liquid water moved downward, and the maximum water flux reached 0.69 cm/h.Following the rainfall, because of ET, the liquid water moved upward at the upper soil profile and moved downward continually over a 54-h period.At the end of the infiltration process, the change in the liquid water showed the same characteristics as before the rainfall.
Figure 10b shows the water vapor characteristics in the temporal-space before and after the same heavy rainfall event.It can be seen that the rainfall event only influenced the shallow soil layer (around a depth of 0 cm to 40 cm) water vapor movement, and the direction was upward.During the other days, the movement direction of water vapor was different during the day (moved upward) and night (moved downward) at a depth of 0 to 40 cm.At a depth of 40 cm to the groundwater table, the water vapor moved downward.This phenomenon indicates that the vapor circulation-condensation process supplied additional moisture for root water uptake at the shallow soil profile.

Movement of Liquid Water and Water Vapor in Autumn
Figure 11 shows the liquid water and water vapor movement in the beginning of autumn.The flux of the liquid water movement in the autumn is much smaller than that in the summer (Figure 11a).However, in arid and semi-arid areas, the air temperature dropped lower than the soil temperature at night during autumn.Consequently, the movement of water vapor shows slightly different

Movement of Liquid Water and Water Vapor in Autumn
Figure 11 shows the liquid water and water vapor movement in the beginning of autumn.The flux of the liquid water movement in the autumn is much smaller than that in the summer (Figure 11a).However, in arid and semi-arid areas, the air temperature dropped lower than the soil temperature at night during autumn.Consequently, the movement of water vapor shows slightly different characteristics compared with that in the summer (Figure 11b).At a depth of 0 cm to 40 cm, the movement of water vapor shows the same characteristics as that in the summer.However, when the air temperature is lower than the groundwater temperature, the water vapor zero flux interface moves to a depth of approximately 80 cm and the water vapor moves upward at a depth of 40 cm to 80 cm, indicating that the water vapor movement is upward.In contrast, when the air temperature is higher than the groundwater temperature, the water vapor moves downward at a depth of 40 cm to groundwater table.
Water 2015, 7, 6999-7021 characteristics compared with that in the summer (Figure 11b).At a depth of 0 cm to 40 cm, the movement of water vapor shows the same characteristics as that in the summer.However, when the air temperature is lower than the groundwater temperature, the water vapor zero flux interface moves to a depth of approximately 80 cm and the water vapor moves upward at a depth of 40 cm to 80 cm, indicating that the water vapor movement is upward.In contrast, when the air temperature is higher than the groundwater temperature, the water vapor moves downward at depth of 40 cm to groundwater table.

Contributions of Liquid Water and Water Vapor to Water Movement
The cumulative change of the average node flux with time is shown in Figure 12.The cumulative upward liquid water and water vapor flux are 26.49cm and 0.03 cm, accounting for 99.89% and 0.11% of the total upward water flux, respectively.Additionally, the cumulative downward liquid water and water vapor flux are 26.75 cm and 0.28 cm, accounting for 99.72% and 0.28% of the total downward water flux, respectively.The liquid water flux is three orders of magnitude higher than the water vapor flux.The water vapor flux is negligible.
Figure 12b shows the change in the water flux with depth.It can be seen that the zero flux interface of liquid water is formed at depths of 73 cm and 70 cm with heat flow and without heat flow.This indicates that the liquid water moves farther downward with the effect of heat flow during the simulated period.Moreover, the cumulative upward water vapor is lower than the downward water vapor (Figure 12a), so the cumulative water vapor was directed to the deep soil layer at a depth of approximately 158 cm.

Contributions of Liquid Water and Water Vapor to Water Movement
The cumulative change of the average node flux with time is shown in Figure 12.The cumulative upward liquid water and water vapor flux are 26.49cm and 0.03 cm, accounting for 99.89% and 0.11% of the total upward water flux, respectively.Additionally, the cumulative downward liquid water and water vapor flux are 26.75 cm and 0.28 cm, accounting for 99.72% and 0.28% of the total downward water flux, respectively.The liquid water flux is three orders of magnitude higher than the water vapor flux.The water vapor flux is negligible.
Figure 12b shows the change in the water flux with depth.It can be seen that the zero flux interface of liquid water is formed at depths of 73 cm and 70 cm with heat flow and without heat flow.This indicates that the liquid water moves farther downward with the effect of heat flow during the simulated period.Moreover, the cumulative upward water vapor is lower than the downward water vapor (Figure 12a), so the cumulative water vapor was directed to the deep soil layer at a depth of approximately 158 cm.

Soil Water Contributions to Evapotranspiration
According to the previous calculations under the influence of heat flow, the soil water storage is increased by 11.7 cm.The precipitation is 41.1 cm, and the ET is 43.8 cm, including 27.9 cm of transpiration and 15.9 cm of evaporation.The total bottom water flux is 14.4 cm (upward), indicating the contribution of groundwater to the soil water balance.Water 2015, 7, 6999-7021

Soil Water Contributions to Evapotranspiration
According to the previous calculations under the influence of heat flow, the soil water storage is increased by 11.7 cm.The precipitation is 41.1 cm, and the ET is 43.8 cm, including 27.9 cm of transpiration and 15.9 cm of evaporation.The total bottom water flux is 14.4 cm (upward), indicating the contribution of groundwater to the soil water balance.
To calculate the soil water contributions to evapotranspiration, the difference in the soil distribution between hydrostatic and actual conditions was used as described by Shah, N., et al. [40].When distribution of the pressure and the soil moisture reached the hydrostatic equilibrium condition, the soil acted as a vessel and the ET was supplied entirely by the groundwater without a vadose zone contribution (VZC).As the depth of the groundwater table (DWT) increased, the hydraulic connections weakened.The hydraulic connections between groundwater are lost at a rate that exceeds the upward replenishment from the saturated zone.Hence, the VZC to ET occurs in a time step of ∆t = t i ´ti´1 of the water content profile from hydrostatic equilibrium.Mathematically (modified from Shah N., et al. [39]), VZC " rpTSM eq ´TSM model q|t i´1 ´pTSM eq ´TSM model q|t i s{pt i ´ti´1 q (19) where ET vzc is the contribution of soil water to ET, P is precipitation, q bot is the calculated bottom flux, and VZC is the contribution of soil water.TSM eq is the soil water content in the column corresponding to DWT under hydrostatic equilibrium conditions, and TSM model is the soil water content computed from simulated by Hydrus-1D for the corresponding time.
The calculated soil water contributions are shown in Figure 13.It can be seen that during the small rainfall period, the soil is dry and the ET is supplied by the groundwater.In contrast, after the rainfall, the soil becomes wet and the ET mainly comes from the soil water.This result confirmed the previous research results using the soil water balance approach at the same site [23].However, during the simulation with heat flow, the cumulative amount of soil water and groundwater contribution to ET are 32.0 cm and 11.8 cm, which account for 73.1% and 26.9%, respectively, of the total ET (43.8 cm).In contrast, the cumulative soil water and groundwater contributions to ET without heat flow are 26.6 cm and 18.2 cm, which account for 59.4% and 40.6%, respectively, of the total ET (44.

Water Movement during Climate and LAI Change
Table 6 shows the calculated flux of liquid water and water vapor during periods of climate and LAI change compared to the results of the calibrated model.It can be seen that the maximum cumulative upward liquid water flux of LAI increased by 0.5, indicating that more water will be consumed if the vegetation becomes denser.However, the minimum cumulative downward liquid water was observed Water 2015, 7, 6999-7021

Water Movement during Climate and LAI Change
Table 6 shows the calculated flux of liquid water and water vapor during periods of climate and LAI change compared to the results of the calibrated model.It can be seen that the maximum cumulative upward liquid water flux of LAI increased by 0.5, indicating that more water will be consumed if the vegetation becomes denser.However, the minimum cumulative downward liquid water was observed when the rainfall decreased by half.This result indicated that water use by the investigated bush will increase and that the percolation amount decreased as rainfall decreased.The cumulative water vapor flux did not change much in comparison with the liquid water flux with decreases in rainfall, increases in temperature, or increases in LAI.This result is consistent with the results presented in Figure 12, suggesting that the amount of water vapor was only slightly affected by a decrease in rainfall and by increases in temperature and LAI.Therefore, the vapor from groundwater will not allow a resilience of the investigated bush ecosystem in case of rainfall decrease, temperature increase and LAI increase.

Discussion
In arid and semi-arid areas, thermally driven soil water movement is an important component of the soil-plant-atmosphere interaction.The temperature gradient changes due to diurnal temperature variations, resulting in the evaporation and condensation of water vapor in the soil.Zeng, et al. [41] found that the movement of water vapor can be described by three stages in shallow soil layers over the course of a day.Furthermore, the evaporation and condensation of water vapor are determined by the pressure head and temperature gradient.The results of water vapor movement simulation in this paper show similar characteristics as previously reported [41] during summer, when the air temperature was higher than the temperature of the groundwater (Figure 10b).However, the zero flux interface is formed in deeper soil in the autumn (Figure 11b).This phenomenon indicates that water vapor condensation is more frequently caused by both diurnal and seasonal air temperature changes under shallow groundwater table conditions in semi-arid desert regions.During the summer, both the air temperature and the groundwater temperature increase, so the alternation and condensation of vapor occur in shallow soil layers determined by diurnal soil temperature fluctuations.In contrast, the air temperature decreased and the groundwater temperature increased due to the soil temperature penetration time lag, resulting in the alternation of evaporation and condensation in shallow and deep soil layers.
Anthony, et al. [42] examined the magnitude of the water vapor flux reported in experimental studies by various authors and found that the maximum magnitude of moisture change due to vapor flow ranges from 7.2 ˆ10 ´3 to 2.5 ˆ10 ´2 cm/h.The simulated results in this study are in agreement with these values (Figures 10b and 11b).Regarding the water vapor contribution to total water flux, Parlange, et al. [43] found that water vapor contributed between 10% and 30% of the total flux based on one week of field observations in bare sandy loam soil.Deb, et al. [44] found that the total upward water vapor flux in the unsaturated soil layer of furrow-irrigated sandy loam soil was approximately 10.4% of the upward total water flux over a 50-day period.Compared to those results, the water vapor contribution to the total water flux in our simulation was much smaller: the upward and downward contributions were 0.11% and 0.28%, respectively.These smaller percentages are caused by the transpiration of the investigated bush from liquid water in semi-arid areas with shallow groundwater table conditions.In other words, the liquid water was dominant due to the Water 2015, 7, 6999-7021 large amount of transpiration, even though the rainfall decreased by half (Table 6).This phenomenon agrees with the results reported by Garcia, et al. [15].Although the amount of water vapor flux is small, the mean pressure over the root zone is very low during dry periods with heat flow (Figure 7).Water vapor could provide a small but noteworthy source of water for plants during the driest period of the year [45].Thus, simulations without heat flow did not capture the dynamics of pressure head changes that affect the rate and direction of water flow [15].
In semi-arid and arid regions, rainfall variations can lead to changes in water use strategy by plant species [11].This study confirmed that an increase in rainfall could lead to an increased use of rain water and decreased groundwater use (Figure 13) during the growing season by the Salix bush.Nevertheless, quantifying groundwater use is another issue.The contribution of groundwater to vegetation water use varies by crop [46], grass [47], and trees [39] according to the Richards equation.However, these studies did not consider the effect of thermal driving flow.Based on our simulation results, the groundwater contribution to ET is different depending on heat flow condition (Figure 10).
This study is limited to thermally driven water movement on water use of S. psammophila.However, water use of ecological system includes competition and coordination mechanism.To evaluate vegetation competition for water and light limitation, Brolsma, et al. [48,49] created a coupled bio-physical-vegetation growth and variably saturated three-dimensional hydrological model to simulate both short-term and long-term vegetation dynamics along a hill slope in temperate forests.It found that once a tree is established under slightly unfavorable soil moisture conditions it cannot be outcompeted by smaller trees with better soil moisture uptake capabilities, both in dry as in wet conditions.To investigate water flow within and around a root network, Schneider, et al. [50] proposed a standalone root water uptake model to investigate the role of root architecture on the spatial distribution of root water uptake.Model simulations show a redistribution of water uptake from more densely to less densely rooted layers with time.Water use of the investigated bush vegetation system may involve both the competition and coordination mechanism referring to the bush ecological consistent as described in section 2.1.Moreover, Specific ecological features could also impact water-use strategies as point by Bertrand, et al. [51].Therefore, these vegetation water use strategies should be taken into account in further modeling efforts.

Conclusions
A Hydrus-1D model was calibrated and validated using hourly field measurements of soil water content; the temperature over a period of 154 days; and liquid water, water vapor and heat transport.The soil water content and temperatures computed by the model fit well with the empirical values at all depths.The model was used to investigate the effects of heat flow on the soil water flux and the groundwater contribution to evapotranspiration by the Salix bush.The main conclusions of the study can be summarized as follows: 1.
The mean pressure head over the root zone when heat flow is included is smaller than when heat flow is not included on dry days.This means that simulated heat transport varied with temperature and pressure head gradients, suggesting a mechanism for moisture redistribution within the root zone.

2.
The zero flux interface of thermally driving water vapor flux varies daily and is affected seasonally by temperature gradients.During the summer, water vapor is condensed in shallow soil layers.However, the zero flux interfaces are formed in both shallow and deeper soil layers in the autumn.This will help us understand the pattern of water vapor movement over space and time.

3.
In semi-arid areas, the water use of the Salix bush depends on rainfall infiltration.During the driest period, more groundwater is used for transpiration.

4.
The percentages of groundwater contribution to ET were 26.9% and 40.6% with heat and without heat flow over the course of the 154-day simulation, respectively.Therefore, the Water 2015, 7, 6999-7021 groundwater contribution will be overestimated when thermally driven water vapor flow is not taken into account.

Figure 1 .
Figure 1.Location of the Hailiutu River catchment and the research site.

Figure 3 .
Figure 3. Summary of the modeled soil column boundary conditions (P = precipitation; ETp = potential evapotranspiration; Tp = transpiration; GWT = depth to groundwater table).The upper and lower boundary conditions for heat transport were specified as temperature boundary conditions.The upper boundary values were given as the observed temperatures at 0.02 cm, and the lower boundary values were given as the groundwater temperature.The initial soil water content and soil temperature were determined from the measured values on April 26 at 0:00 h by interpolating the measured values between different depths.

Figure 4 .Figure 4 .
Figure 4.The fit of measured soil water content at eight different depths based on soil profile (points) using model-computed values (solid lines) for the calibration period (26 April to 12 July) and the validation period (13 July to 27 September).

Figure 5 .
Figure 5.The fit of the measured (SFm) and computed (SFs) cumulative transpiration rates of the S. psammophila bush.

Figure 6 .
Figure 6.The fit of computed soil temperatures (solid lines) to measured temperatures (points) at nine depths in the soil profile used for calibration (26 April to 12 July) and validation (13 July to 27 September).

Figure 7 .
Figure 7.Comparison of the mean pressure head values over the root zone with heat flow (hRoot with heat) and without heat flow (hRoot without heat).

Figure 7 .
Figure 7.Comparison of the mean pressure head values over the root zone with heat flow (hRoot with heat) and without heat flow (hRoot without heat).

Water 2015, 7 , 18 Figure 8 .
Figure 8.(a) Cumulative evapotranspiration; and (b) root water uptake for selected depth intervals with heat flow and without heat flow.

Figure 9 .
Figure 9.A comparison of the cumulative change in soil water storage and cumulative bottom water flux, both with and without heat flow.

Figure 8 .
Figure 8.(a) Cumulative evapotranspiration; and (b) root water uptake for selected depth intervals with heat flow and without heat flow.

3. 4 . 18 Figure 8 .
Figure 8.(a) Cumulative evapotranspiration; and (b) root water uptake for selected depth intervals with heat flow and without heat flow.

Figure 9 .
Figure 9.A comparison of the cumulative change in soil water storage and cumulative bottom water flux, both with and without heat flow.

Figure 9 .
Figure 9.A comparison of the cumulative change in soil water storage and cumulative bottom water flux, both with and without heat flow.

Figure 10 .
Figure 10.Simulated vertical distribution of the isothermal and thermal liquid water and water vapor fluxes during a typical period from the experimental period before and after a rainfall event (0:00 26 June to 23:00 1 July); (a) liquid water fluxes; (b) thermal vapor water fluxes.

Figure 10 .
Figure 10.Simulated vertical distribution of the isothermal and thermal liquid water and water vapor fluxes during a typical period from the experimental period before and after a rainfall event (0:00 26 June to 23:00 1 July); (a) liquid water fluxes; (b) thermal vapor water fluxes.

Figure 11 .
Figure 11.The movement of (a) liquid water; and (b) water vapor in autumn.

Figure 11 .
Figure 11.The movement of (a) liquid water; and (b) water vapor in autumn.

Water 2015, 7 21 Figure 12 .
Figure 12.Changes in the cumulative liquid water and water vapor flux with time (a) and with depth; (b) during the investigation period.

Figure 12 .
Figure 12.Changes in the cumulative liquid water and water vapor flux with time (a) and with depth; (b) during the investigation period.
8 cm) without heat flow.The simulated results indicate that the contribution of groundwater to ET is overestimated without considering heat flow.Water 2015, 7 22 simulation with heat flow, the cumulative amount of soil water and groundwater contribution to ET are 32.0 cm and 11.8 cm, which account for 73.1% and 26.9%, respectively, of the total ET (43.8 cm).In contrast, the cumulative soil water and groundwater contributions to ET without heat flow are 26.6 cm and 18.2 cm, which account for 59.4% and 40.6%, respectively, of the total ET (44.8 cm) without heat flow.The simulated results indicate that the contribution of groundwater to ET is overestimated without considering heat flow.

Figure 13 .
Figure 13.The contribution of soil water to ET during the simulation period (ETvzc is the soil water contribution, and ETcvzc is the cumulative contribution): (a) soil water movement with heat flow; and (b) soil water movement without heat flow.

Figure 13 .
Figure 13.The contribution of soil water to ET during the simulation period (ET vzc is the soil water contribution, and ET cvzc is the cumulative contribution): (a) soil water movement with heat flow; and (b) soil water movement without heat flow.

Table 1 .
Instruments used at the research site.

Table 1 .
Instruments used at the research site.

Height (m) Horizontal Distance to the Bush Stem (m) Instrument or Sensor Type
2.2.1.Measurement of Soil Water Content

Table 2 .
Soil texture at the observation site.

Table 3 .
The calibrated parameters for the HYDRUS-1D model.

Table 4 .
Statistical measures of HYDRUS-1D model performance for simulations of volumetric water content.

Table 5 .
Statistical measures of HYDRUS-1D model performance for simulations of soil temperature at different soil depths.

Table 6 .
The cumulative flux of liquid water and water vapor during climate and LAI changes.