Monitoring the variations of evapotranspiration due to land use/cover change in a semiarid shrubland

Evapotranspiration (ET) is an important process in the hydrological cycle, and vegetation change is a primary factor that affects ET. In this study, we analyzed the annual and inter-annual characteristics of ET using continuous observation data from eddy covariance (EC) measurement over 4 years (1 July 2011 to 30 June 2015) in a semiarid shrubland of Mu Us Sandy Land, China. The Normalized Difference Vegetation Index (NDVI) was demonstrated as the predominant factor that influences the seasonal variations in ET. Additionally, during the land degradation and vegetation rehabilitation processes, ET and normalized ET both increased due to the integrated effects of the changes in vegetation type, topography, and soil surface characteristics. This study could improve our understanding of the effects of land use/cover change on ET in the fragile ecosystem of semiarid regions and provide a scientific reference for the sustainable management of regional land and water resources.


Introduction
Arid and semiarid biomes cover approximately 40 % of the Earth's terrestrial surface (Fernández, 2002).Previous studies have shown that more than 50 % of precipitation (P) is consumed by evapotranspiration (E T ) (Yang et al., 2007;Liu et al., 2002).Moreover, a slight change in E T could have significant influences on water cycle and the ratio of E T /P could increase to even 90 % or more in these regions (Mo et al., 2004;Glenn et al., 2007).In terms of physical processes, E T is affected by net radiation (Valipour et al., 2015), water vapor pressure deficit (Zhang et al., 2014), wind speed (Falamarzi et al., 2014), and soil water stress (Allen et al., 1998).
Vegetation change mainly includes phenological change (temporal) and land use/cover change (spatial).Phenological change reflects the response of plants to climate change (vegetation greening and browning processes) (Ge et al., 2015), which actively controls E T through internal physiologies such as stomatal conductance (Pearcy et al., 1989), as well as the number and sizes of stomata (Turrell, 1947).In general, transpiration is directly proportional to stomatal conductance at the leaf scale (Leuning et al., 1995).At the canopy scale, E T is positively proportional to surface conductance, which is an integration of stomatal conductance and leaf area (Ding et al., 2014).Thus, as a good indicator of vegetation phenological change, many studies have found that E T is positively related to vegetation indexes such as the Normalized Difference Vegetation Index (NDVI) (Gu et al., 2007).Land use/cover change influences E T by modifying vegetation species with different transpiration rates, radiation transfers within the canopy (Martens et al., 2000;Panferov et al., 2001), topography (Lv et al., 2006), albedos (Zeng and Yoon, 2009), soil texture (Maayar and Chen, 2006), litter coverage (Wang, 1992), and biological soil crusts (BSCs) (Yang et al., 2015;Fu et al., 2010;Liu, 2012).These complex processes result in no consensus on the effects of land use/cover change on E T .For example, during the land degradation process, some researchers found that warming air temperature was the main cause of making E T increase (Zeng and Yang, 2008;Li et al., 2013;Feddema and Freire, 2001).By contrast, a decline in E T was found along with the deforestation process because of less transpiration (Snyman, Published by Copernicus Publications on behalf of the European Geosciences Union. 2001; Souza and Oyama, 2011) or higher albedo (Zeng et al., 2002).Moreover, no changes in E T during the land degradation process were reported either (Hoshino et al., 2009).Thus, there has been an important push to better understand how E T responds to vegetation change, especially to the land use/cover change.
Three methods were usually employed to assess the effects of vegetation change on E T : numerical models, paired comparative approaches, and in situ field observations.In these methods, numerical models are widely used (Twine et al., 2004;Kim et al., 2005;Li et al., 2009;Cornelissen et al., 2013;Mo et al., 2004).However, model parameterization of vegetation conditions is a big challenge, as the aforementioned complex underlying mechanisms may not be completely considered in the models.Therefore, the simulated effects of vegetation change on E T are highly dependent on model parameterizations, which may induce uncertainty (Cornelissen et al., 2013;Li et al., 2009).The paired comparative approach is often considered the best method; nonetheless, it is difficult to find two sites with similar meteorological conditions but different vegetation conditions (Li et al., 2009;Lorup et al., 1998).Moreover, the method of in situ field observations is widely used to investigate long-term landatmosphere exchanges.However, the land use/cover conditions at sites are generally stable, and only the response of E T to vegetation phenological change can be observed, such as the E T variations in grassland (Y.Zhang et al., 2005), mixed plantation (cork oak, black locust, and arborvitae) (Tong et al., 2017), vineyard (Li et al., 2015), and grazed steppe (Chen et al., 2009;Vetter et al., 2012).Continuous field observations under both land degradation and vegetation rehabilitation processes have rarely been documented, especially in the semiarid shrubland.
The Mu Us Sandy Land is a semiarid shrubland ecosystem on the northern margin of the Loess Plateau in China.The area covers only 40 000 km 2 (Dong and Zhang, 2001) and is ecologically fragile (Yang et al., 2007).In such an ecosystem, sand dunes and BSCs are commonly observed (Gao et al., 2014;Yang et al., 2015;Li and Li, 2000;Liu, 2012).Due to the existence of BSCs and dry sand layers (Z.Wang et al., 2006;Feng, 1994;Liu et al., 2006;Yuan et al., 2008), soil evaporation has been effectively retained; therefore, the Mu Us Sandy Land contains abundant groundwater (Li and Li, 2000).During the past decades, rapid land use/cover changes have occurred in this region due to agricultural reclamation (Wu and Ci, 2002;Ostwald and Chen, 2006;Zhang et al., 2006), leading to dramatic changes in vegetation conditions.With respect to the specific question of whether land use/cover change will lead to increases in E T or not, a continuous measurement of E T under different land use/cover conditions is required in this region.Coincidentally, two processes of land use/cover changes (land degradation and vegetation rehabilitation) have occurred at the edge of the Mu Us Sandy Land, providing us with a unique opportunity to study the effects of land use/cover change on E T .
Hence, based on the 4-year measurement of E T by eddy covariance techniques, this study analyzed the seasonal and inter-annual variations in E T , and discussed the possible reasons for the responses of E T to land use/cover change.Our results were expected to provide a scientific reference for the sustainable management of regional land and water resources in the context of intensive agricultural reclamation.
2 Case study and data

Site description
The study was conducted at the Yulin flux site (38 • 26 N, 109 • 28 E, 1233 m), which was established in June 2011.This site is located in a landform transition zone that changes from the Mu Us Sandy Land to the north Shaanxi Loess Plateau (Fig. 1).This site is a semiarid area with temperate continental monsoon climate.According to long-term climate data  from a meteorological station in Yulin (Fig. 1), the annual precipitation varied from 235 to 685 mm, with a mean of 402 mm, and more than 50 % of annual precipitation fell in the monsoon season (July-September).The mean annual air temperature was 8.4 • C over the past 61 years.The dominant soil type is sand (98 % sand) (saturated soil water content of 0.43 m 3 m −3 , field capacity of 0.16 m 3 m −3 , residual moisture content of 0.045 m 3 m −3 ).There are widely distributed fixed sand dunes and semi-fixed sand dunes around the site, and the depth of the dry sand layer is 10 cm (Z.Wang et al., 2006).The mean groundwater depth at our study site from 1 July 2011 to 30 June 2015 was 3.5 m.
Shortage of water is the critical limiting factor for vegetation growth in this site, and drought-enduring vegetation (e.g., shrubs) prevails as a result of droughts (Wang et al., 2002;Wu, 2006).The study site is mainly covered with mixed vegetation: the native drought-enduring shrubs with low water demand (e.g., Artemisia ordosica and Salix psammophila) (Fig. 2a) and the sparse grass (mainly distributed at the bottom of sand dunes because of the better soil moisture condition) (Lv et al., 2006).The maximum root depth of the shrubs was approximately 160 cm.Xiao et al. (2005) reported that the growing season of Artemisia ordosica and Salix psammophila spanned from late April to late September.Therefore, we defined the period from 1 May to 30 September as the vegetation growing season for data analysis in this study.On 15 August and 7 September 2011, we did surveys of the vegetation coverage by randomly selecting seven samples around the flux tower (5 × 500 cm × 500 cm and 2 × 1000 cm × 1000 cm).We found that the vegetation coverage was 28.2 % in August and 27.9 % in September.At the end of June 2012, the land use/cover condition around the eastern portion of the flux tower began to be changed by farmers (leaves and branches were cut, and the sand dunes were bulldozed) (Fig. 2c), converting part of the natural vegetated land to bare land, with the planning of planting potatoes in the future.As time went on, natural grass gradually grew out in the area of bare land before potatoes were planted.Thus, our study period (1 July 2011 to 30 June 2015) was divided into four periods according to the land use/cover conditions: (a) Period I (1 July 2011 to 30 June 2012), the period with the natural land use/cover condition (i.e., mixed sparsely distributed shrubs and grass) (Fig. 2a and b); (b) Period II (1 July 2012 to 30 June 2013), the transitional period when the land use/cover condition started to change (some natural vegetation removed and sand dunes bulldozed); (c) Period III (1 July 2013 to 30 June 2014), the period when the land use/cover condition constituted two parts: the natural vegetation zone and the bare soil zone (Fig. 2c); and (d) Period IV (1 July 2014 to 30 June 2015), the period when the bare soil zone was gradually covered by regrowing grass (Fig. 2d).

Eddy covariance system measurements
Net exchange of water vapor between atmosphere and canopy at this site is measured by the eddy covariance (EC) flux measurement, which assesses the fluxes of land-atmosphere (such as water and energy) (Baldocchi et al., 2001).The data are essential for the estimation of the water and energy balance (Franssen et al., 2010).At our site, the EC system is installed at a height of 7.53 m above the ground surface, using CSAT3 three-dimensional sonic anemometers (Campbell Scientific Inc., Logan, UT, USA) for wind and temperature fluctuation measurements and a LI-7500A openpath infrared gas analyzer (LI-COR, Inc., Lincoln, NE, USA) for water vapor content measurement.

Other measurements
Net radiation (R n ) is measured by a net radiometer (CNR-4; KIPP&ZONEN, Delft, the Netherlands), including four radiometers measuring the incoming and reflected short-wave radiation (R S ), and incoming and outgoing long-wave radiation (R L ).Sunshine duration (D S ) is measured by a sunshine recorder (CSD3; KIPP&ZONEN, Delft, the Netherlands).Wind speed and direction (05103, Young Co. Traverse City, MI, USA) are measured at 10 m above the ground surface.Precipitation (P, mm) is recorded with a tipping bucket rain gauge (TE525MM; Campbell Scientific Inc., Logan, UT, USA) installed at a height of 0.7 m above the ground surface.Air temperature (T a ) and relative humidity (R H ) are measured by a temperature and relative humidity probe (HMP45C; Campbell Scientific Inc., Logan, UT, USA) at a height of 2.6 m above the ground surface.Soil water content (θ ) is measured by time domain reflectometry (TDR) sensors (CS616; Campbell Scientific Inc., Logan, UT, USA), soil temperature (T s ) is measured by thermocouples (109; Campbell Scientific Inc., Logan, UT, USA), and soil heat flux (G) is measured by heat flux plates (HFP01SC; Campbell Scientific Inc., Logan, UT, USA) at a depth of 0.03 m below the ground surface.These ground variables (G, θ , T s ) are measured beneath the surface at two profiles: a plant canopy profile and a bare soil profile.θ and T s are measured at depths of 5, 10, 20, 40, 60, 80, 120, and 160 cm below the ground surface.The groundwater table is measured by an automatic sensor (CS450-L; Campbell Scientific Inc., Logan, UT, USA), which is installed in a groundwater well close to the tower.

Flux data processing
The 10 Hz three-dimensional wind speed and water vapor concentrations that were collected by the EC technique were processed to half-hourly latent heat flux (λE T ) using Eddypro processing software (v5.2.0, LI-COR, Lincoln, NE, USA).The main principle is that λE T can be expressed as ρ a w q (where w is the fluctuation of vertical wind speed, q is the fluctuation of specific humidity, and ρ a is the air density).The software also applies the quality control of data, including spike removal, tilt correction, time lag compensation, turbulent fluctuation blocking, and spectral corrections.The percentages of half-hourly λE T values removed (includ- ing missing and rejected) through the quality control procedure were 17.3 % in Period I, 20.2 % in Period II, 16.5 % in Period III, and 18.6 % in Period IV.Almost all the removed λE T values occurred during the nighttime (89.1 % in Period I, 91.3 % in Period II, 92.6 % in Period III, and 88.7 % in Period IV).During the nighttime, the change in λE T was small, and E T values were close to zero.Therefore, after removal of the nighttime λE T values, the errors of the gapfilled nighttime values based on the neighboring good data were small.Moreover, nighttime λE T values accounted for only a small proportion of the daily E T .Furthermore, the percentages of rejected and missing data in our study are similar to those reported by other scholars, and these percentages are in a range of 15-31 % (Falge et al., 2001;Wever et al., 2002;Mauder et al., 2006).Therefore, the λE T data set was considered reliable after a quality control procedure.
After quality control, missing and rejected data were gap-filled in order to create continuous data sets.Three methods were applied in the gap-filling procedure: (1) linear interpolation was used to fill gaps of less than 1 h by calculating an average of the values before and after the data gap; (2) for gaps that are larger than 1 h but smaller than 7 days, the mean diurnal variation (MDV) method (Falge et al., 2001) was used; (3) for gaps that are larger than 7 days but smaller than 15 days in daily λE T values, we fitted the relationship between daily λE T (W m −2 ) and the daily available energy flux (R n − G) (W m −2 ) in each period.We chose the function f with the highest coefficient of correlation (R) in each period (Yan et al., 2013), and the function Then, we used the fitted function f in each period to estimate the daily λE T values of large gaps.In addition, gaps that are larger than 7 days but smaller than 15 days mostly appeared in the winter, which accounted for a small proportion of annual λE T .

Footprint model
In order to determine the contributing source area of flux at our site, the scalar flux footprint model proposed by Hsieh et al. (2000) was used.The analytic model accurately describes the relationship between the footprint, observation height, surface roughness, and atmospheric stability.The fetch F f was calculated as follows: where k is the von Karman constant (= 0.40), D and Q are similarity constants (for stable conditions, D = 0.28 and Q = 0.59; for near neutral and neutral conditions, D = 0.97 and Q = 1; for unstable conditions, D = 2.44 and L is the Obukhov length, Z m is the height of the wind instrument (= 10.0 m), and Z u is defined as (Hsieh et al., 2000) where Z 0 is the height of the momentum roughness (0.05 m).

Method of analyzing controlling factors on E T
It is generally recognized that potential evapotranspiration (E TP ), vegetation condition, and soil water stress are the three main factors that control E T (Lettenmaier and Famiglietti, 2006;Chen et al., 2014).In order to decouple the effect of vegetation change from the integrated effects of these three factors on E T , we used a simple equation which was similar to the FAO single crop coefficient method (Irrigation and Drainage Paper No. 56 (FAO-56)).This equation can be expressed as follows: where f v (vegetation) represents the effect of vegetation change on E T and f s (soil water) represents the effect of soil water stress on E T .Moreover, f v (vegetation) can be regarded as the normalized E T , which eliminates the effects of atmospheric and soil water stress on E T and can be expressed by rearranging Eq. ( 3): (4)

Potential evapotranspiration
E TP (mm day −1 ) was estimated by the following equation (Maidment, 1993), which is a modification of the Penman equation: where the units of R n and G are mm day −1 ; ρ a is the air density (= 3.486 P a 275+T a , kg m −3 , where P a is the atmospheric pressure in kPa and T a is in degrees Celsius); c p is the specific heat of moist air (= 1.013 kJ kg −1 • C −1 ); is the slope of the saturation vapor-pressure-temperature curve (kPa • C −1 ); VPD is the difference between the mean saturation vapor pressure (e s , kPa) and actual vapor pressure (e a , kPa); and λ is the latent heat of vaporization of water (= 2.51 MJ kg −1 ).γ is the psychrometric constant (kPa • C −1 ), which is calculated by the following equation: where ε is the ratio of the molecular weight of water vapor to that of dry air (= 0.622).
r a is the aerodynamic resistance, which can be calculated as follows (Penman, 1948): where Z h is the height at which meteorological variables are measured (2 m), and Z ao is the aerodynamic roughness of the surface (0.00137 m) (Penman, 1963); U 2 is the daily wind speed at a height of 2.0 m (m s −1 ), and it was calculated by the wind speed at a height of 10.0 m (U 10 , m s −1 ):

Vegetation parameters
In this study, vegetation phenology was represented by Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI data when the land use/cover conditions were fixed.NDVI is sufficiently stable to reflect the seasonal changes of any vegetation (Huete et al., 2002).Higher NDVI generally reflects the greater photosynthetic capacity (greenness) of the vegetation canopy (Gu et al., 2007;Tucker, 1979).The daily NDVI was calculated by daily surface reflectance data: where NIR is the spectral response in the near-infrared band (857 nm) and VIS is the visible red radiation band (645 nm).
In this study, NDVI was calculated by using MODIS/Terra data (MOD09GQ) (NDVI Terra ) and MODIS/Aqua data (MYD09GQ) (NDVI Aqua ) (http://reverb.echo.nasa.gov),respectively.As we found that there were slight differences (|NDVI Terra − NDVI Aqua | = 0.01 ± 0.0075) between NDVI Terra and NDVI Aqua , we calculated NDVI by averaging NDVI Terra and NDVI Aqua in order to eliminate the impacts of such differences.The calculated NDVI values were then filtered to remove anomalous hikes and drops (Lunetta et al., 2006), and the smoothing spline method was used to produce a smoother profile.
Theoretically, land use/cover change can be evaluated by comparing the land use/cover maps in two different periods.However, transient land use/cover maps were unavailable at our site.Therefore, we separated the study area within the footprint into two zones: the undisturbed zone without any land use/cover change was deemed zone A and the disturbed zone with land use/cover change was deemed zone B. In zone A, vegetation change included only vegetation phenological change; however, in zone B, there were not only vegetation phenological changes, but also land use/cover changes.Based on the assumption that the phenological changes caused by climate in the two zones were the same, we defined an indicator (D lu ) as a measure of land use/cover change: where M A and M B are the monthly vegetation coverages of zone A and zone B, respectively.The monthly vegetation coverage was calculated by monthly NDVI values (Gutman and Ignatov, 1998): where NDVI max is the maximum value (0.8 in this study) and NDVI min is the minimum value (0.05 in this study) (Gutman and Ignatov, 1998).The calculated monthly M values (27.6 and 24.2 %) were consistent with the measured vegetation coverages in August 2011 (28.2 %) and September 2011 (27.9 %) at our study site.

Soil water stress
The effects of the soil water stress on E T can be described in three stages (Idso et al., 1974).Stage 1: the soil water is enough to satisfy the potential evaporation rate (f s = 1); stage 2: the soil is drying and water availability limits E T (0 < f s < 1); and stage 3: the soil is dry and evaporation can be considered negligible (f s = 0).We used daily soil water content in the root depth (θ r ) to estimate f s by the following expression (Morillas et al., 2013): where θ w is the wilting value and θ k is the stable field capacity which is considered to be equivalent to 60 % of the field capacity (Lei et al., 1988;Wang et al., 2008).θ r was calculated by measured soil water contents at different depths (θ i where i = 5, 10, 20, 40, 60, 80, 120, and 160 cm).From land surface to the depth of 5 cm, the soil water profile was assumed triangular, while at other depths, the soil water profiles were assumed trapezoidal.Therefore, the soil moisture of root zone was calculated as where θ i (i = 5, 10, 20, 40, 60, 80, 120, and 160 cm) was calculated by taking a weighted average of the measured values in the canopy and bare surface patches, where θ i,c and θ i,b refer to the measured soil water contents of the canopy patch and bare soil patch at the depth of i cm, respectively.

Statistical analysis
In this study, we chose daily data in Period I to analyze the correlations between E T and the three controlling factors (E TP , NDVI, and f s ).We used several common functions (e.g., an exponential function, a linear function, a logarithmic function, and a quadratic function) to fit these correlations.We found that the determination coefficient (R 2 ) of the linear function was generally the highest.Therefore, in this study, we chose the linear function to fit the correlations between E T and the three controlling factors.Additionally, a significant t test was performed to evaluate the degrees of these correlations.Moreover, data on rainy days were removed because E T values were gap-filled rather than measured.

Footprint and energy balance closure
Based on the footprint model, we got the half-hourly scatter data (Eq.2), and, according to the wind rose diagram (Fig. 3a), the prevailing wind directions at this site were northwesterly and southeasterly.Therefore, we chose an ellipse to enclose the scatters and simulate the footprint (Fig. 3b).Under unstable conditions, 93 % of the half-hourly flux data are plotted within the ellipse.Additionally, we measured the boundary of zone B in October 2013 when the land use/cover condition in zone B had stopped changing (Fig. 3b).There were 11 pixels (250 m × 250 m) in zone A and 19 pixels (250 m × 250 m) in zone B, and thus, when calculating the weight-averaged NDVI (NDVI w ) within the footprint, we chose the weighted coefficient as β = 11/(11 + 19).
EC system performance was assessed by the energy balance closure which was calculated by conducting the linear regression between available energy (R n −G) and the sum of surface fluxes (λE T + H), which is also used to examine the quality of flux data (Wilson et al., 2002).The linear regression yielded a slope of 0.87, an intercept of −1.42 W m −2 , and an R 2 of 0.82.These indicators suggested that the measurements at our experimental site provided reliable flux data and that the EC measurements underestimated the sum of the surface fluxes to the extent of 13 %.Many researchers have investigated energy imbalance (Barr et al., 2006;Wilson et al., 2002;Franssen et al., 2010), and there is a consensus that it is difficult to examine the exact reasons for the imbalance.

Characteristics of environmental variables
A brief summary of key environmental variables is presented in this section.Four-year andlong-term (1954-2014) average monthly values of D s , T a , R H , and P are shown in Fig. 4. Monthly D s was much higher than the long-term average monthly values, except in July and September.The highest value of D s was observed in May (299.5 h) and the low-   The inter-annual characteristics of daily T a , D s , R H , θ r , groundwater level (GWL), and total P in the growing season of each period are listed in Table 1.
The values of T a , R H , P, and θ r in the growing season of Period IV were the lowest compared to those in the other three periods.Periods I-III were all wet years, while Period IV was a dry year.The values of θ r in Periods I-III were similar; however, θ r decreased by 0.0113 m 3 m −3 in Period IV.The mean GWL in Period III was the shallowest.

Seasonal variations in E T due to climate variability and vegetation phenology
The seasonal curve of E T in each year had a single peak value (Fig. 5a), with higher E T appearing mostly in the growing season, while lower E T appeared in the non-growing season.
The daily E T ranged from 0.0 mm day −1 to 6.8 mm day −1 during the four periods; the highest E T was observed on 22 June 2013, which was the day after a continuous rainfall event that extended from 19 to 21 June 2013 (90.3 mm).
The lowest E T appeared on 28 November 2012, which was in the frozen period (late November to early March at our study site).On rainy days, E TP (Fig. 5b) was low due to low net radiation and air temperature.E TP ranged from 0.2 mm day −1 in December 2011 to 17.9 mm day −1 in September 2013.The seasonal NDVI curve for natural land use/cover conditions (in zone A during Periods I-IV and in zone B during Period I) represented the process of natural vegetation phenology, and it had a single peak value in each year (Fig. 5c).In early May, the seasonal NDVI curve began to increase as the native vegetation entered the growing season, and a maximum value (0.27 ± 0.01) was reached in July or August.In the winter, the daily NDVI remained relatively constant (0.13 ± 0.01).f s (Fig. 5d) increased rapidly in response to rainfall events of more than 5 mm a day and decreased rapidly 1 or 2 days after rainfall events.From late November to early March, there was a frozen period when the soil water content was below the wilting point.The groundwater level changed obviously in the monsoon season (July to September) and mildly in the winter (December to February).
The linear correlations between E T and the three controlling factors all passed the t test at a 95 % confidence level.The R 2 value of the correlation between E T and NDVI w (NDVI w = NDVI A × β + NDVI B × (1 − β)) was the largest, indicating that NDVI was highly correlated with the daily variations in E T .To better quantify the effects of the phenological process on E T , the correlation between daily f v and NDVI w in Period I was analyzed (Fig. 7a).
A positive linear regression was found between f v and NDVI w (Fig. 7a).The slope of the linear regression was used to evaluate the degree of the correlation between f v and the vegetation phenological process.We found that when w increased by 1 unit, f v increased by approximately 1.86 units.

Inter-annual variations in E T due to land use/cover change
During the four periods, in zone A, the NDVI values of each period were similar because the land use/cover condition did not change, while in zone B, the peak values of NDVI first declined from 0.28 to 0.15 (Period I to Period III) due to the land use/cover condition changing from mixed vegetation to bare soil.The peak NDVI values then increased to 0.22 (Period IV) due to grass recovery (Fig. 5c).An interesting phenomenon was observed accompanied by the changing process of land use/cover conditions: E T in the growing season gradually increased from Period I to Period III (Table 2), while it increased greatly in Period IV even with less precipitation, because a mass of soil water and groundwater was consumed to satisfy the E T demand (Fig. 5e).
Compared with Period I, D lu values in Period II and Period III gradually increased, while D lu in Period IV decreased.Taking August in each period as an example, in Period I, D lu was 0.2 %, while in Periods II-IV, D lu were 2.9, 12.6, and 8.6 %, respectively.In order to eliminate the influence of vegetation phenological change on E T , we chose the Hydrol.Earth Syst.Sci., 21, 863-877, 2017 www.hydrol-earth-syst-sci.net/21/863/2017/  growing season of each period to analyze the correlation between f v and D lu .
The quantitative results of the correlation between D lu and f v are shown in Fig. 7b.From Period I to Period III, as land surface characteristics changed (the natural vegetation in zone B was cleared, the fixed and semi-fixed sand dunes were bulldozed, and the BSCs and dry sand layers disappeared), f v increased, and this increase was more evident in Period III (from 78.5 to 88.1).When the land use/cover conditions in zone B gradually changed from bare soil to sparse grassland due to the self-restoring capacity of nature, f v increased significantly (from 88.1 to 111.3).Table 2. Typical values of total evapotranspiration (E T , mm), total potential evapotranspiration (E TP , mm), the indicator of land use/cover change (D lu , %), the soil water stress of the root zone (f s ), and normalized E T (f v (= E T /(E TP ×f s ))) in the growing season of each period (because there were some missing data in Period IV (from 12 September to 23 November 2014 and from 13 March to 22 April 2015), we removed the values of E T , E TP , and f s in these two time ranges of Periods I-III).The correlations between E T and its controlling factors suggest that at our experimental site, NDVI is the predominant factor that influences the seasonal variations in E T .The positive linear relationship between f v and NDVI suggests that transpiration is likely controlled by the stomatal conductance and the numbers of stomata, which are proportional to the leaf area (Pearcy et al., 1989;Turrell, 1947), rather than the atmospheric water demand represented by E TP .
Various studies have assessed the correlation between vegetation phenological change and E T , and these results generally reflected consistent and positive linear relationships (Nouri et al., 2014;Rossato et al., 2005;Duchemin et al., 2006;Glenn et al., 2007).However, for different vegetation species, phenological change has effects on E T to different degrees.Relatively strong regressions between NDVI and E T have been reported at forested sites (Loukas et al., 2005;Nouri et al., 2014;Lo Seen Chong et al., 1993) and grass-covered sites (Kondoh and Higuchi, 2001;Nouri et al., 2014), with determination coefficients higher than 0.7.These results reflect the strong control between phenological changes and E T .Thus, we speculate that for high vegetated ecosystems, phenological change may have a significant control on E T .However, in low vegetated ecosystems such as the sparse shrubland in this study, the relationship between E T and phenological change is thus positive but relatively weak.

Possible reasons for the effects of land use/cover changes
During Periods I-IV, the land use/cover conditions at our experimental site underwent changes associated with two processes: land degradation process (Periods II-III) and vegetation rehabilitation process (Period IV).Notable results were observed during these two processes: (1) E T and normalized E T values both increased and (2) normalized E T increased much faster during the vegetation rehabilitation process than it did during the land degradation process.
The effect of phenological change on E T demonstrates that E T decreases with leaf browning.Thus, we expect that E T will also decrease if leaves are cleared by human activities.However, during Periods II-III, not only were leaves cleared, but other land surface properties were also changed (all branches were cut, sand dunes (fixed and semi-fixed) Hydrol.Earth Syst.Sci., 21, 863-877, 2017 www.hydrol-earth-syst-sci.net/21/863/2017/ were bulldozed, and the dry sand layers and BSCs were destroyed), resulting in complex land use/cover conditions.These altered land surface properties might contribute to the increase in E T .Previous studies demonstrated that dry sand layers and BSCs could effectively restrict the soil evaporation rate (Z.Wang et al., 2006;Lv et al., 2006;Liu et al., 2006;Dong et al., 1999;Yang et al., 2015;Fu et al., 2010;Liu, 2012).However, the bulldozing of sand dunes at our experimental site made the elevation of the flat soil surface lower than the average elevation of the undisturbed soil surface (approximately 1.5 m lower, Fig. 2d), making the groundwater depth much shallower than the pre-disturbance depth.Thus, the formation of dry sand layers was restricted due to the shallow groundwater level.In this situation with the destroyed BSCs and the disappeared dry sand layers, the sufficient groundwater supply (Li and Li, 2000) accelerated the loss of water that was stored in shallow soil through evaporation.The enhanced soil evaporation offset the inhibiting effect of transpiration due to leaves clearing, which made E T increase.
A secondary reason for the increase in soil evaporation was that the soil layer absorbed more solar radiation during the land degradation process.In Period I, the radiation absorbed by the shadowed soil was the solar radiation transmitted into the canopy of shrubs and grass.However, when the natural vegetation was cleared, the leaves and the branches were also removed, which made the shadowed soil exposed and enhanced the radiation absorbed by the soil, thereby increasing soil evaporation (Martens et al., 2000;Panferov et al., 2001).Moreover, the removal of leaves and branches and the disappearance of sand dunes both altered the land surface albedo, which could have directly altered the solar radiation absorbed by the land surface (Dirmeyer and Shukla, 1994;Greene et al., 1999), subsequently leading to the change in E T .
Some inconsistent results regarding the E T dynamics during land degradation process were reported.A portion of studies reported that E T decreased during the land degradation process for different reasons.For example, Souza and Oyama (2011) and Snyman (2001) demonstrated that E T decreased during the land degradation process due to decreased transpiration in semiarid regions.Lu et al. (2011) considered that the low soil water content was the main reason for the decrease in E T during the land degradation process.Mao and Cherkauer (2009) also reported a decrease in E T when the land use/cover condition was converted from forest to grass or cropland in the Great Lakes region.However, contrasting results were also reported regarding the effects of land degradation on E T .Hoshino et al. (2009) found that there was no difference in E T during the land degradation process associated with overgrazing in a semiarid Mongolian grassland, and they hypothesized that the reason for this lack of change might be the short grazing time (2 years).Li et al. (2013) demonstrated that the warming air temperature was the main cause of increased E T during the land degradation process on the Qinghai-Tibet Plateau.Throughout the above studies of E T during land degradation processes, we found it difficult to accurately describe the trends in E T , even when the land degradation was only manifested by less vegetation coverage.Therefore, at our study site with complex land surface properties (sand dunes, dry sand layers, and BSCs), the effect of land degradation on E T was much more complicated.
During the vegetation rehabilitation process (Period IV), f v increased significantly due to the rehabilitation of grass in zone B, even though less precipitation was observed compared with other periods (Periods I, II, and III).The rehabilitation of grass, rather than shrubs, was due to the sufficient groundwater supply, which resulted from bulldozing the sand dunes.Previous researchers reported that sparse shrubs more commonly grew at the top of sand dunes and that grass grew at the bottom of sand dunes because the difference between groundwater level and the top of sand dunes was larger than that between groundwater level and the bottom of the sand dunes (Lv et al., 2006;Dong et al., 1999).Because transpiration increases with vegetation greening (as demonstrated in Sect.4.3), the regrowing grass would enhance plant transpiration supplied by the sufficient groundwater.More importantly, the transpiration rate of grass is higher than that of shrubs because shrubs are easier to survive in water-limited conditions (Yang et al., 2014;Wang et al., 2002;Wu, 2006).Therefore, in the vegetation rehabilitation process, the enhancement of the transpiration rate in Period IV was much higher than that in Periods I-III.Similar conclusions regarding increased E T due to the enhanced transpiration during the vegetation rehabilitation process were reported (Qiu et al., 2011;Yang et al., 2014;Sun et al., 2006;Li et al., 2009).Meanwhile, the regrowing grass could reduce the radiation absorbed by the soil and hence reduce soil evaporation.However, the interception of radiation by the grass canopy was expected to be smaller than that by the mixed shrub and grass canopy in Periods I-III because the leaf area index of grass was smaller than the sum of leaf area and stem area indexes of the mix of shrubs and grass.Therefore, the reduction in soil evaporation in Period IV might be small compared with the increase in soil evaporation in Periods I-III.
We noticed that the GWL decreased continuously from Period III to Period IV due to the enhanced E T by the regrowth of grass and relative low precipitation, and the regrowing grass has a higher transpiration rate than that of the native mixed shrub and grass.Therefore, we hypothesize that if the land use/cover condition of zone B continues to be grassland over the next several years, the groundwater level will decrease due to the larger consumption, making the soil water condition gradually become poorer for the growth of grass.Then, in this situation, the grassland is expected to degrade to shrubland in zone B because shrubs are easier to survive in water-limited ecosystems.Furthermore, in the next few years, potatoes will be planted in zone B. However, the water requirement of potato is more than 320 mm in the growing season (Qin et al., 2013;Liu et al., 2010) and the wa-ter consumption is more than that of natural grass (Qin et al., 2013(Qin et al., , 2014;;Hou et al., 2010).Thus, irrigation is necessary for planting potatoes during the growing season in waterlimited ecosystems (Liu et al., 2010;Fabeiro et al., 2001).Our results imply that the groundwater level might continue to decrease faster with the growth of potatoes in the future, which may lead to a more fragile ecosystem.

Conclusion
In this study, seasonal and inter-annual features of E T were analyzed.Daily E T was in a range from 0.0 to 6.8 mm day −1 during the four periods.NDVI was the predominant factor that influences the seasonal variations in E T , and vegetation greening had a positive effect on E T .During the land degradation process (Periods II-III), when natural vegetation (including leaves and branches), sand dunes, dry sand layers, and BSCs were all bulldozed, E T increased at a mild rate.During the vegetation rehabilitation process (Period IV) with less precipitation, E T increased at a faster rate than that in the degradation process.Our study demonstrated that when the land use/cover condition was changed by human activities, the underlying mechanisms that influence E T were complex, and vegetation type, topography, and soil surface characteristics may all contribute to the changes in E T .Furthermore, our results suggest that when we simulate the effects of land use/cover change on hydrological processes, the vegetation factor might not be the unique factor to parameterize; instead, the integrated effects of land surface and vegetation conditions should be considered.Our study also provides a scientific reference to the regional sustainable management of water resources in the context of intensive agricultural reclamation.

Data availability
Due to the shorter observational time, the measured data in our flux station are not directly open to the public.However, the readers can ask the authors for the data.
Competing interests.The authors declare that they have no conflict of interest.

Figure 2 .
Figure 2. Land use/cover conditions at the study site: (a) the natural land use/cover condition of shrubland (photo was taken on 6 August 2011); (b) the natural land use/cover condition of grassland (photo was taken on 7 September 2011); (c) the undisturbed zone (natural vegetation) and the disturbed zone (bare soil) in the land degradation process (photo was taken on 26 April 2013); (d) the undisturbed zone (natural vegetation) and the disturbed zone (grassland) during the vegetation rehabilitation process (photo was taken on 16 August 2014).

Figure 3 .
Figure 3. Diagrams of the wind rose and footprint: (a) wind rose of the study site by using half-hourly wind speed and wind direction data and (b) simulated footprint by ellipse (the long axis is 1682 m, and the short axis is 1263 m; zone A is the source area in which the land use/cover condition did not change, while zone B is the source area in which the land use/cover condition did change due to human activities; the white triangle is the flux tower).

Table 1 .
Daily air temperature (T a , • C), relatively humidity (R H , %), total sunshine duration (D S , h), soil water content of the root zone (θ r , m 3 m −3 ), the groundwater level (GWL, m), and total precipitation (P, mm) in 1954-2014 and in the growing season of each period (because there were some missing data in Period IV (from 12 September to 23 November 2014 and from 13 March to 22 April 2015), we excluded data in these two time ranges of Periods I-III and 1954-2014).

Figure 5 .
Figure 5. Seasonal and inter-annual characteristics of daily (a) evapotranspiration (E T , mm); (b) potential evapotranspiration (E TP , mm); (c) NDVI in zone A and zone B within the footprint; (d) the soil water stress of the root zone (f s ); and (e) the groundwater level (GWL, m) from 1 July 2011 to 30 June 2015.

Figure 6 .
Figure 6.The correlations between daily evapotranspiration (E T , mm) and its controlling factors: (a) daily potential evapotranspiration (E TP , mm); (b) daily weight-averaged NDVI (NDVI w ) within the footprint; and (c) daily soil water stress of the root zone (f s ) in Period I by excluding the data on rainy days (r: Pearson's correlation coefficient; T : t test significance).

Figure 7 .
Figure 7. Quantitative analysis of the correlations between (a) vegetation phenological change (NDVI w ) and daily normalized E T (f v = E T /(E TP × f s )) in Period I (excluding the data on rainy days and frozen days) and (b) the indicator of land use/cover change (D lu ) and total normalized E T (f v = E T /(E TP × f s )) in the growing season of each period.