Estimation of Groundwater Evapotranspiration of Different Dominant Phreatophytes in the Mu Us Sandy Region

Groundwater evapotranspiration (ETG) estimation is an important issue in semiarid areas for groundwater resources management and environmental protection. It is widely estimated by diurnal water table fluctuations. In this study, the ETG at four sites with different plants was estimated using both diurnal water table and soil moisture fluctuations in the northeastern Mu Us sandy region, in order to identify the groundwater utilization strategy by different dominant phreatophytes. Groundwater level was monitored by ventilatory pressure transducers (Solinst LevelVent, Solinst Canada Ltd.; accuracy ±3 mm), while soil moisture was monitored using EM50 loggers (Decagon Devices Inc., Pullman, USA) in K1 and K14 and simulated by Hydrus-1D in other observation wells. A significant spatial variation of ETG was found within a limited area, indicating a poor representativeness of site ETG for regional estimation. The mean values of ETG are 4.01 mm/d, 6.03 mm/d, 8.96 mm/d, and 12.26 mm/d at the Achnatherum splendens site, Carex stenophylla site, Salix psammophila site and Populus alba site, respectively, for the whole growing season. ETG is more sensitive to depth to water table (DWT) in the Carex stenophylla site than in the Achnatherum splendens site for grass-dominated areas and more sensitive to DWT in the Populus alba site than in Salix psammophila site for tree-dominated areas. Groundwater extinction depths are estimated at 4.1 m, 2.4 m, 7.1 m, and 2.9 m in the Achnatherum splendens site, Carex stenophylla site, Salix psammophila site and Populus alba site, respectively.


Introduction
In arid and semi-arid areas, restoring and protecting the environment is a huge project for the benefit of mankind [1]. While groundwater is a crucial source of ecological restoration because of the lower precipitation and higher evapotranspiration in arid regions. However, these areas also depend heavily on groundwater for industry, drinking water supply and agriculture [2,3]. Reasonable allocation of groundwater resources in all sectors is thus a major environmental issue.
Groundwater is a main water source for phreatophytes in arid and semiarid areas. Evapotranspiration partitioning shows that transpiration accounts for the majority of evapotranspiration [4,5]. Over the past few decades, vegetation coverage has displayed a large increase because of the multiple ecological restoration programs including the grain for green and afforestation in some arid areas of China [6]. However, ecological restoration by afforestation may increase groundwater depth and create potentially large ecological and water costs in arid and semiarid [7]. For example, Gribovszki et al. [8] found that the water table beneath forest was 0.4-0.5 m lower than that beneath grassland due to the afforestation in the Hungarian Great Plain. So, the groundwater consumption strategy of each phreatophyte should be considered in the process of afforestation.
In addition, groundwater evapotranspiration (ET G ) is a vital factor for groundwater balance analysis in arid and simi-arid regions. For example, almost 60% of natural groundwater discharge is through ET G in the Ordos Plateau in NW China [9]. Many studies usually estimate ET G indirectly, i.e., by the pan coefficient method [10], the remote sensing method [11] and the eddy covariance method [12]. The pan coefficient method is simple, but coefficients mainly are empirical with low precision. The remote sensing method is usually used to calculate evapotranspiration at a regional scale and the time step depends on remote sensing data. It cannot analyze the temporal and spatial variation of evapotranspiration at small scale. The eddy covariance method depends on expensive instruments, so it is not widely used. In addition, all of these methods estimate rather total evapotranspiration than ET G. Evapotranspiration partitioning is needed to exclude the contribution of soil water, which is a complicated process. Therefore, over the past few decades, the method of estimating evapotranspiration based on diurnal water table fluctuations, developed by White [13], has attracted more and more attention [14][15][16][17][18], because it can estimate ET G directly with fewer parameters and lower costs.
The Mu Us sandy region is a part of the Ordos Plateau in NW China and is used to be one of the most seriously desertified areas in China [19]. It has also undergone a significant vegetation restoration in recent decades [20]. While Cheng et al. [21] estimated the ET G at four sites in the center of the Mu Us sandy region, using only one observation at each site. However, the spatial variation of ET G at point scale under each vegetation is unclear and the representativeness of the ET G for one observation well is also uncertain. Multiple wells are expected to improve the accuracy of ET G estimation in each vegetation site [22]. Therefore, the objectives of this study are (1) to estimate the ET G from multiple observation wells based on the water table and soil moisture fluctuations at different sites; (2) to study the temporal and spatial variations of ET G ; and (3) to analyze the sensitivity of ET G of different vegetation to water table depth.

Site Description
The study site (39 • 16 54"~39 • 18 48" N, 108 • 47 40"~108 • 50 03" E) is located in the northeastern Mu Us sandy region, within a wetland of the Mukai Lake in Ordos City, Inner Mongolia, NW China. The geomorphology is mainly lake basin that contains some microrelief, such as wetlands, crescent sand dunes, and highlands ( Figure 1). Vegetation is distributed unevenly in the various microrelief. Vegetation coverage is higher in wetlands, up to about 70%. While, there is only sparse vegetation on sand dunes. The climate is semi-arid characterized by a low rainfall and high evapotranspiration. According to the long-term meteorological data (2000-2017) from the Wushenzhao meteorological station near the study area, the mean annual precipitation is~330 mm and over 70% of the annual precipitation falls from June to September. Average reference evapotranspiration (ET 0 ) calculated by the Penman-Monteith equation [23] is~1150mm/a. Average annual air temperature is about 7.0 • C, with a maximum of 22.4 • C in July and a minimum of −11.0 • C in January. Population density is low and arable lands near the study area are scarce, so the consumption of groundwater by human activities is negligible.
The main species of vegetation are Achnatherum splendens (Figure 2a

Monitoring Design
Fifteen shallow observation wells (K1 to K15) were drilled by a hand auger arou the lake (Figure 1). K1, K2, and K3 are in the Achnatherum splendens site; K4, K9, K10, K and K15 are in the Carex stenophylla site; K5 and K13 are in the Caragana korshinskii site; K7 are in the Salix psammophila site; K8 and K11 are located in sand dunes; K12 is in Populus alba site. The depth of the wells varies between 2 m and 4 m, depending on depth to water table (DWT) at each site. Each was constructed using a 50 mm PVC p screened over the entire subsurface length. Augured soil was backfilled around the P pipe and granular bentonite was then packed around land surface to avoid preferen flow.

Monitoring Design
Fifteen shallow observation wells (K1 to K15) were drilled by a hand auger around the lake (Figure 1). K1, K2, and K3 are in the Achnatherum splendens site; K4, K9, K10, K14, and K15 are in the Carex stenophylla site; K5 and K13 are in the Caragana korshinskii site; K6, K7 are in the Salix psammophila site; K8 and K11 are located in sand dunes; K12 is in the Populus alba site. The depth of the wells varies between 2 m and 4 m, depending on the depth to water table (DWT) at each site. Each was constructed using a 50 mm PVC pipe screened over the entire subsurface length. Augured soil was backfilled around the PVC pipe and granular bentonite was then packed around land surface to avoid preferential flow.
Groundwater level was continuously monitored from 1 April 2019 to 31 October 2019 at one-hour intervals by ventilatory pressure transducers with data logging capability (Solinst LevelVent, Solinst Canada Ltd.; accuracy ±3 mm) that automatically eliminated the effect of air pressure fluctuations by a vented cable to the wellhead. The well cap assembly of the transducers serves to prevent direct rainfall input and evaporation loss. In order to remove noise, a five-point moving average was used to smoothen groundwater

Monitoring Design
Fifteen shallow observation wells (K1 to K15) were drilled by a hand auger around the lake ( Figure 1). K1, K2, and K3 are in the Achnatherum splendens site; K4, K9, K10, K14, and K15 are in the Carex stenophylla site; K5 and K13 are in the Caragana korshinskii site; K6, K7 are in the Salix psammophila site; K8 and K11 are located in sand dunes; K12 is in the Populus alba site. The depth of the wells varies between 2 m and 4 m, depending on the depth to water table (DWT) at each site. Each was constructed using a 50 mm PVC pipe screened over the entire subsurface length. Augured soil was backfilled around the PVC pipe and granular bentonite was then packed around land surface to avoid preferential flow.
Groundwater level was continuously monitored from 1 April 2019 to 31 October 2019 at one-hour intervals by ventilatory pressure transducers with data logging capability (Solinst LevelVent, Solinst Canada Ltd.; accuracy ±3 mm) that automatically eliminated the effect of air pressure fluctuations by a vented cable to the wellhead. The well cap assembly of the transducers serves to prevent direct rainfall input and evaporation loss. In order to remove noise, a five-point moving average was used to smoothen groundwater level data [18]. The data were also calibrated frequently in the field with manual measurements of DWT.
Soil samples were collected from different depths and the particle sizes were determined by the sieving method for the sand fractions and by a soil hydrometer for the silt and clay fractions. Results showed that the soil texture of unsaturated zone at the study sites was mainly fine sand and homogeneous in the vertical direction at each well ( Table 1). The soil textures of all observation wells were divided artificially into two categories, with sand content higher than 90%, and with sand content less than 90%. Then, two observation sites, K1 and K14, were thought to be representative of the two soil textures. Soil moisture (θ s ) was monitored once an hour from 20 cm to 100 cm below ground surface with 20 cm intervals at K1 and from 20 cm to 60 cm with 10 cm intervals at K14 using EM50 loggers (Decagon Devices Inc., Pullman, USA). Meteorological data (i.e., solar radiation, air temperature, air humidity, and wind speed) were monitored at one-hour intervals by a meteorological station (HOBO-U30, Onset, USA) at the southeast of study site. Hourly reference evapotranspiration (ET 0 ) was calculated by the Penman-Monteith Equation (Equation (1)) using the meteorological data [23].
where ET 0 is reference evapotranspiration (mm·h −1 ); R n is net radiation (MJ·m −2 ·h −1 ); G is soil heat flux density (MJ·m −2 ·h −1 ); T h is mean hourly air temperature ( • C); ∆ is saturation slope vapor pressure curve at T h (kPa· • C −1 ); γ is psychrometric constant (kPa· • C −1 ); e • (T h ) is saturation vapor pressure at air temperature T h (kPa); e a is average hourly actual vapor pressure (kPa); µ is average hourly wind speed (m·s −1 ). It is worth noting that ET 0 is the evapotranspiration from the reference surface which is under a hypothetical grass reference crop with an assumed crop height of 0.12 m, a fixed surface resistance of 70 s·m −1 and an albedo of 0.23. The reference surface closely resembles an extensive surface of green, Soil hydraulic parameters were obtained through inverse modeling by Hydrus-1D, which was widely used by previous studies [24][25][26]. For K1 and K14, as the representative of the two soils, both groundwater level and soil moisture were monitored. Therefore, the meteorological data, water table and soil moisture data were used to inverse soil hydraulic parameters of the two soil textures in the study site.
For the inverse modeling, a 150 cm homogeneous soil profile was divided uniformly into 150 elements, with 151 nodes and two observations corresponding to the soil moisture monitoring depths (20 cm and 100 cm at K1; 20 cm and 40 cm at K14). The maximum number of iterations was set to 50. The van Genuchten-Mualem hydraulic model was selected and the soil hydraulic parameters of loamy sand in the model were used as the initial parameters because the soil texture is classified as loamy sand based on particle size analysis. The maximum and minimum values of each parameter, i.e., θ s , θ r , ϕ, n, and Ks, are set to limit the range of inversion. According to the particle analysis of soils in the study site (Table 1), the maximum values of ϕ, n and Ks were set to the default parameters of sand in Hydrus-1D code and the minimum values equaled to the default parameters of silt, as shown in Table 2. Data for inverse solution were soil moisture data measured in the field at the two observations depths. Time discretization was in hours and initial time step was set to 0.01. The upper boundary condition was defined as an atmospheric boundary condition. The potential evapotranspiration (ET p ) equals to ET 0 as recommended by Šimůnek et al. [27], and the potential evaporation (E p ) and potential transpiration (T p ) were estimated from ET p using Beer's law that partitions the solar radiation component of the energy budget via interception by the canopy [28,29] as follows in Equations (2) and (3): where k is a constant governing the radiation extinction by the canopy (-), the value is about 0.6 according to Belmans et al. [30]. LAI is the leaf area index and the values in grass site and tree site are 1.9 and 2.6, respectively. The lower boundary condition was variable pressure head that was monitored in the study site. The period of the inverse model was from 20 May 2019 to 5 June 2019 for K1 and from 24 May 2019 to 6 September 2019 for K14, when obvious diurnal fluctuations of water table and soil water moisture continuously occurred.
In the study site, the root depth was about 150 cm. An exponential root distribution model developed by Zeng [31] was used in our model, as follow in Equation (4): Water 2021, 13, 440 where Y is the cumulative root fraction from the surface to depth d; a and b are the empirical coefficients. For short grass, a and b are 10.74 and 2.61, respectively. The S-shaped water stress function of root water uptake was adopted [32].

Simulated Soil Moisture
Soil moisture was monitored at K1 and K14 and simulated by Hydrus-1D at other observation wells, because of the lack of monitoring equipment. The simulation domain is a 150 cm homogeneous soil profile and was divided uniformly into 150 elements with 151 nodes. The simulated period was from 1 April 2019 to 31 October 2019 and the soil hydraulic parameters were from the above-mentioned inverse modelling. The upper boundary condition was also atmospheric boundary condition. E p and T p were also estimated by Equations (2) and (3). The lower boundary condition was also defined as a variable pressure head condition. Root depth was 150 cm in the grass site and 200 cm in the areas with arboreous vegetation. The root water uptake model and root distribution model are same to the models used in the previous section. While the parameters a and b in the root distribution model are 7.07 and 1.95, respectively, for deciduous needleleaf tree, i.e., Salix psammophila, and are 5.99 and 1.96, respectively, for deciduous broadleaf tree, i.e., Populus alba, according to Zeng [31]. The observation points were set from the land surface to water table with 10 cm intervals to obtain soil moisture along the soil profile. To test the accuracy of the simulated data, the soil moisture of K1 (60 cm and 80 cm depth) and K14 (50 cm and 70 cm depth) were also simulated and then compared to the measured data.

Statistical Analysis
The root mean square error (RMSE), Nash-Sutcliffe efficiency (NSE) and percent bias (PBIAS) were used to compare the observed and simulated soil moisture. They are defined as Equations (5)-(7), respectively, as followed: where θ obs is the observed soil moisture in field, θ sim is the simulated soil moisture, n is the number of observations, aveθ obs is the average soil moisture of all the observed events. RMSE indicates the degree of dispersion of the simulated soil moisture relative to the observed moisture [33], and a higher value indicates a higher error in the simulated moisture. NSE represents the extent to which the simulated moisture and the observed moisture follow a 1:1 slope line [34]. It ranges from −∞ to 1 and higher values indicate a better agreement with the observed moisture [35]. PBIAS measures the average tendency of the simulated values to be larger or smaller than their observed ones. The optimal value of PBIAS is 0.0, with small values indicating accurate model simulation. Positive values indicate overestimation bias, and negative values indicate model underestimation bias.

ET G Estimation
White [13] initially developed a method to quantify daily ET G using diurnal water table fluctuations (hereafter referred to as the White method [36]), as follows in Equation (8): where ET G is the daily groundwater evapotranspiration (L·T −1 ), S y is the specific yield of sediments (−), r is the rate of water table recovery during nighttime (L·T −1 ) and s is the daily change in storage (L·T −1 ) In this study, a modified White method was adopted to estimate the daily ET G . For the method, daily hydrostatic equilibrium recovery (R e ) is proposed to correct the groundwater inflow estimated by the White method, taking the disequilibrium of soil water into consideration. The hydrostatic equilibrium at the capillary zone is disturbed at the moment when root water uptake commences in the morning, then the degree of deviation from hydrostatic equilibrium increases in daytime due to the higher ET rate than the hydrostatic equilibrium recovery rate. So, capillary water continues to be replenished by groundwater to lower the degree of deviation from hydrostatic equilibrium at night, resulting in groundwater inflow being underestimated by water table recovery in the White method. Thus, daily R e is estimated by extending the capillary water recovery rate at nighttime to the whole day, representing the error of daily groundwater inflow estimated by the White method. The procedure for calculating R e is quite similar to that of the recovery rate (r in Equation (8)) in the White method. The first step is the calculation of the deviation of total soil water (TSW) from hydrostatic equilibrium, (TSW eq − TSW s ), at two specified times (i.e., 22:00 and 6:00 in this study) during nighttime, respectively, then the difference between (TSW eq − TSW s ) of both points of time is determined, and finally the daily R e is calculated by multiplying by 24 h. Equations (9) and (10) are the mathematical expressions for R e : R e = 24 × [(TSW eq − TSW s ) t − (TSW eq − TSW s ) t+∆t ]/∆t (9) and where TSW eq is the TSW in the column for the corresponding DWT under hydrostatic equilibrium conditions (L); TSW s is the TSW computed from the field measured water content values (L); θ is water content (L 3 ·L −3 ); j is the soil moisture monitoring point of the soil column. Then the daily ET G was estimated as follows in Equation (11): The soil moisture under a hydrostatic equilibrium condition was calculated using the van Genuchten soil water retention properties (Equations (12) and (13)).
where h is the water pressure head (L), θ s and θ r are the saturated and residual soil moisture (L 3 ·L −3 ), respectively, ϕ is the air entry value also referred as bubble pressure (L −1 ), n is the pore size distribution parameter (−), m =1 − 1/n (−), and ϕ, n, and m are the empirical shape factors in the water retention function. It should be noted that specific yield (S y ) is the most crucial factor for ET G estimation by water table fluctuations [37][38][39]. S y is defined as the volume of water released under gravity from storage per unit cross sectional area per unit decline in water table [40]. It is affected by soil texture, depth to water table and drainage time [38]. So, the depthdependent and time-dependent S y are widely used to estimate ET G based on water table fluctuations [21,37].
In this study, the drainage time is not taken into account because the change of capillary water has been calculated. Therefore, the equation for depth-dependent S y developed by Duke [41] with van Genuchten parameters was used (Equation (14)) where d is the depth to water table (L).

Soil Hydraulic Parameters
The measured and simulated soil moisture from inverse modelling at two depths for each observation well (K1 and K14) are shown in Figure 3. The RMSE, NSE, PBIAS, and R-square (R 2 ) of measured and simulated data were 0.003 cm 3 /cm 3 , 0.99, −0.04%, and 0.99, respectively, for K1, while for K14, the RMSE, NSE, PBIAS, and R 2 of measured and simulated data were 0.015 cm 3 /cm 3 , 0.86, 1.22%, and 0.81, respectively. According to the statistical analysis, the simulated soil hydraulic parameters of K1 and K14, as showed in Table 3, were a good representation of the actual field parameters.
Water 2021, 13, x FOR PEER REVIEW 8 of

Soil Hydraulic Parameters
The measured and simulated soil moisture from inverse modelling at two depths f each observation well (K1 and K14) are shown in Figure 3. The RMSE, NSE, PBIAS, an R-square (R 2 ) of measured and simulated data were 0.003 cm 3 /cm 3 , 0.99, −0.04%, and 0.9 respectively, for K1, while for K14, the RMSE, NSE, PBIAS, and R 2 of measured and sim ulated data were 0.015 cm 3 /cm 3 , 0.86, 1.22%, and 0.81, respectively. According to the st tistical analysis, the simulated soil hydraulic parameters of K1 and K14, as showed in T ble 3, were a good representation of the actual field parameters.

Validation of the Simulated Soil Moisture
Soil moisture was also simulated at K1 and K14 to assess the accuracy of the modelle soil moisture, using the measured and simulated soil moisture at 60 cm and 80 cm for K and 50 cm and 70 cm for K14. As shown in Figure 4 and Table 4, the simulated soil mo ture has a good agreement with the measured data, with the mean RMSE of 0.012 cm 3 /cm at K1 and 0.019 cm 3 /cm 3 at K14. While the NSE and R 2 at K14 are 0.25 and 0.21, respe tively. The possible reason is that the sensitivity of the monitoring instrument decreas due to the observation depth being close to water table.

Validation of the Simulated Soil Moisture
Soil moisture was also simulated at K1 and K14 to assess the accuracy of the modelled soil moisture, using the measured and simulated soil moisture at 60 cm and 80 cm for K1 and 50 cm and 70 cm for K14. As shown in Figure 4 and Table 4, the simulated soil moisture has a good agreement with the measured data, with the mean RMSE of 0.012 cm 3 /cm 3 at K1 and 0.019 cm 3 /cm 3 at K14. While the NSE and R 2 at K14 are 0.25 and 0.21, respectively. The possible reason is that the sensitivity of the monitoring instrument decreases due to the observation depth being close to water table.
Water 2021, 13, x FOR PEER REVIEW

Temporal and Spatial Variations of Water Table
Obvious diurnal water table fluctuation occurred during the whole growin in the field. As shown in Figure 5a and b, groundwater table fell during daytim covered during nighttime. The variation pattern of water table at sub-daily scale caused by the evapotranspiration of vegetation. While solar radiation is a domin controlling plant water use in sub-daily scale [42,43], and further influences w variations (Figure 5c).
The changes in water table at seasonal scales are reflected in daily amplit the Carex stenophylla site, the average daily amplitudes of water table at K10 ar 2.3 cm, and 1.4 cm in the early, middle, and late of the growing season, respectiv ure 5a). While they are 2.3 cm, 4.1 cm, and 1.3 cm, respectively, at K12 in the Po site ( Figure 5b). Obviously, the daily amplitude of water table in the middle of t ing season is larger than that in the early and late of the growing season. It caused by the vegetation vitality that depends on air temperature (Figure 5d).

Temporal and Spatial Variations of Water Table
Obvious diurnal water table fluctuation occurred during the whole growing season in the field. As shown in Figure 5a and b, groundwater table fell during daytime and recovered during nighttime. The variation pattern of water table at sub-daily scale is mainly caused by the evapotranspiration of vegetation. While solar radiation is a dominant factor controlling plant water use in sub-daily scale [42,43], and further influences water table variations (Figure 5c).
The changes in water table at seasonal scales are reflected in daily amplitudes. For the Carex stenophylla site, the average daily amplitudes of water table at K10 are 1.7 cm, 2.3 cm, and 1.4 cm in the early, middle, and late of the growing season, respectively (Figure 5a). While they are 2.3 cm, 4.1 cm, and 1.3 cm, respectively, at K12 in the Populus alba site (Figure 5b). Obviously, the daily amplitude of water table in the middle of the growing season is larger than that in the early and late of the growing season. It is mainly caused by the vegetation vitality that depends on air temperature (Figure 5d).
The patterns of water table fluctuation in the field are also different for various DWT. In the Achnatherum splendens site, the mean daily amplitude of water table is 0.2 cm with an average DWT of 1.85 m at K2 from 21 May to 25 May (Figure 6a). At K3, the mean daily amplitude of water table is 0.6 cm with an average DWT of 0.71 m in the same period. In addition, the mean daily amplitude of water table at K6 and K7 are 5.8 cm and 1.8 cm with average DWT of 0.99 m and 1.67 m, respectively, from 1 July to 5 July in the Salix psammophila site (Figure 6b). Obviously, the daily amplitude of water table with shallow water table is larger than that with deep water table in the same vegetation condition, mainly resulting from the decrease of root distribution and increase of specific yield with DWT. The patterns of water table fluctuation in the field are also different for various DWT. In the Achnatherum splendens site, the mean daily amplitude of water table is 0.2 cm with an average DWT of 1.85 m at K2 from May 21 to May 25 (Figure 6a). At K3, the mean daily amplitude of water table is 0.6 cm with an average DWT of 0.71 m in the same period. In addition, the mean daily amplitude of water table at K6 and K7 are 5.8 cm and 1.8 cm with average DWT of 0.99 m and 1.67 m, respectively, from July 1 to July 5 in the Salix psammophila site (Figure 6b). Obviously, the daily amplitude of water table with shallow water table is larger than that with deep water table in the same vegetation condition, mainly resulting from the decrease of root distribution and increase of specific yield with DWT.
Vegetation condition is also a vital influencing factor of water level fluctuation pattern. For example, the daily amplitude of water table (3.7 cm) at K12 in the Populus alba site is larger than that (2.0 cm) at K14 in the Carex stenophylla site with the same average DWT of 0.70 m (Figure 6c). Furthermore, K10 is in the Carex stenophylla site and a sand dune at a distance of 20 m from K10 is free of vegetation. There is obviously water table fluctuation at K10 but no water table fluctuation at K11 (Figure 6d). The water table shows no fluctuation at K5 and K13 either, indicating that water table fluctuation is not influenced by other factors, such as tidal [44,45] and barometric pressure [46].

Temporal Variations of ETG
ETG was estimated except for the wells with no diurnal water table fluctuations, i.e K5, K11, and K13, indicating that Caragana korshinskii is a plant independent of groun water. There were no ETG data on rainy days due to the inapplicability of the estimate method. Average daily ETG and the standard deviation of ETG at all the observation wel are showed in Figure 7. Temporal variations of ETG are obvious in the growing seaso Vegetation condition is also a vital influencing factor of water level fluctuation pattern. For example, the daily amplitude of water table (3.7 cm) at K12 in the Populus alba site is larger than that (2.0 cm) at K14 in the Carex stenophylla site with the same average DWT of 0.70 m (Figure 6c). Furthermore, K10 is in the Carex stenophylla site and a sand dune at a distance of 20 m from K10 is free of vegetation. There is obviously water table fluctuation at K10 but no water table fluctuation at K11 (Figure 6d). The water table shows no fluctuation at K5 and K13 either, indicating that water table fluctuation is not influenced by other factors, such as tidal [44,45] and barometric pressure [46].

Temporal Variations of ET G
ET G was estimated except for the wells with no diurnal water table fluctuations, i.e., K5, K11, and K13, indicating that Caragana korshinskii is a plant independent of groundwater. There were no ET G data on rainy days due to the inapplicability of the estimated method. Average daily ET G and the standard deviation of ET G at all the observation wells are showed in Figure 7. Temporal variations of ET G are obvious in the growing season. The average ET G is small (1.5 mm/d) in the early April and gradually increases to peak values (10.6 mm/d) at the end of July, then gradually decreases to 3.5 mm/d at the end of October. The variation trend is consistent with daily average air temperature (Figure 7), indicating that seasonal ET G variation mainly depends on air temperature. However, the higher air temperature corresponds to the lower ET G in April. It may result from the less vitality of vegetation due to the less leaf area in the early growing season [42].

Temporal Variations of ETG
ETG was estimated except for the wells with no diurnal water table fluctuations, i.e K5, K11, and K13, indicating that Caragana korshinskii is a plant independent of ground water. There were no ETG data on rainy days due to the inapplicability of the estimate method. Average daily ETG and the standard deviation of ETG at all the observation wel are showed in Figure 7. Temporal variations of ETG are obvious in the growing seaso The average ETG is small (1.5 mm/d) in the early April and gradually increases to pea values (10.6 mm/d) at the end of July, then gradually decreases to 3.5 mm/d at the end October. The variation trend is consistent with daily average air temperature (Figure 7 indicating that seasonal ETG variation mainly depends on air temperature. However, th higher air temperature corresponds to the lower ETG in April. It may result from the le vitality of vegetation due to the less leaf area in the early growing season [42].

Spatial Variations of ETG
Spatial variation is significant indicated by the high standard deviation of the ET The standard deviation shows a similar variation pattern to the average ETG, which small in the early and late growing season and large in the middle growing season. Th large deviation (4.2 mm/d on average) in the middle of the growing season mainly d pends on vegetations types and DWT of the observation wells. For example, the mean ET in June was 10.9 mm/d at K12 site with Populus alba, which was 1.7 times that of K14 (6

Spatial Variations of ET G
Spatial variation is significant indicated by the high standard deviation of the ET G . The standard deviation shows a similar variation pattern to the average ET G , which is small in the early and late growing season and large in the middle growing season. The large deviation (4.2 mm/d on average) in the middle of the growing season mainly depends on vegetations types and DWT of the observation wells. For example, the mean ET G in June was 10.9 mm/d at K12 site with Populus alba, which was 1.7 times that of K14 (6.4 mm/d) in the Carex stenophylla site with the same average DWT of 0.65 m. On the other hand, the mean ET G was 2.9 mm/d at K2 site with a mean DWT of 2.03 m and 4.1 mm/d at K3 site with a mean DWT of 0.87 m in June for the Achnatherum splendens site. The spatial variations of ET G were also significant although the distance of two observations was small, for instance, the mean ET G of K10 was 5.8 mm/d over the growing season in the Carex stenophylla site while there was no ET G at K11 in a sand dune even though it is only 20 m away from K10.
The higher standard deviation of ET G is also observed from the wells at the same vegetation site (  The large spatial difference of ET G at different wells reflects that ET G estimated by one observation well can only represent point-scale ET G . However, Butler et al. [42] used a point-scale ET G to quantify the ET G of an area of 150 m 2 . According to our analysis, it may cause a large error to study regional ET G by a point value and the mean ET G of multiple wells in an area will be more representative, even though in the same vegetation site. Therefore, more observation wells are needed in an area to obtain a region-scale ET G . It is obvious that the ET G in the tree sites is about twice as much as that of the grass sites. The same relationship has been found in previous studies [21,22]. For example, Satchithanantham et al. [22] found that the average ET G at trees sites and grass sites were 4.1 mm/d and 2.1 mm/d, respectively, in the Assiniboine River drainage basin in southwestern Manitoba, Canada. Although the average ET G of trees and grass also shows two to one ratio, both of them are much less than that of our study. There are several possible reasons. Firstly, the average DWT in the Assiniboine River drainage basin (1.4 m for trees sites and 1.0 m for grass sites) was larger than that in our study site (1.1 m and 0.5 m for trees and grass sites, respectively). While groundwater absorbed by vegetation decreases with DWT [43]. Secondly, the annual precipitation in the Assiniboine River drainage basin was 474 mm, being 44% larger than that of our study site (330 mm) due to the climate difference. Vegetation will use more soil water from rainfall infiltration than groundwater because the soil water is rich in nutrients and oxygen [47,48]. Finally, the ET G estimated method in our study considers the daily hydrostatic equilibrium recovery in capillary zones that is considered to be a part of ET G .

Sensitivity of ET G to Influencing Factors
ET 0 reflects the combined effect of individual meteorological variables because it is determined by solar radiation, air temperature, relative humidity, and wind speed (Equation (1)). As shown in Figure 8, ET G has a linear positive correlation with ET 0 in each observation well. The slope of the linear expression of them are considered to be the parameter that reflects the relationship between ET G and ET 0 , and they are different among the wells. The slope in K6, K7, and K12 are bigger than one (ranging from 1.03 to 1.78), indicating that the daily ET G in these wells are generally larger than daily ET 0 . The possible explanation is that the vegetation in K6, K7, and K12 are trees that have larger ET G than grass. While K4 is close to the lake with dense grass, resulting in the larger ET G . The slope in the other wells ranges from 0.46 to 0.91. It reflects that the proportion of ET G and ET 0 shows obvious spatial difference. However, the coefficient of determination (R 2 ) of the linear expressions is not suitable for all wells, ranging from 0.06 to 0.67, resulting from other factors influencing ET G . For example, R 2 are 0.06 and 0.08 in K2 and K7, respectively, where DWT is larger.
wells. The slope in K6, K7, and K12 are bigger than one (ranging from 1.03 to 1.78), cating that the daily ETG in these wells are generally larger than daily ET0. The pos explanation is that the vegetation in K6, K7, and K12 are trees that have larger ETG grass. While K4 is close to the lake with dense grass, resulting in the larger ETG. The s in the other wells ranges from 0.46 to 0.91. It reflects that the proportion of ETG and shows obvious spatial difference. However, the coefficient of determination (R 2 ) o linear expressions is not suitable for all wells, ranging from 0.06 to 0.67, resulting other factors influencing ETG. For example, R 2 are 0.06 and 0.08 in K2 and K7, respecti where DWT is larger. In order to obtain the relationship between ETG and DWT, a ratio of ETG to (ETG/ET0) are used to exclude the influence of climatic conditions. For various vegeta types, the exponential decay function between ETG/ET0 and DWT are obtained as sh in Figure 9, but the parameters of these functions are different. The maximum ETG/E 1.1 for the grass sites (the Carex stenophylla site and the Achnatherum splendens site), w the maximum ETG/ET0 is 3.6 for the tree sites (the Populus alba site and the Salix psamm ila site), resulted from the larger ETG in the tree sites. The coefficient of the indepen variable (CIV) in the exponential function can reflect the decay rate of the ETG/ET0 In order to obtain the relationship between ET G and DWT, a ratio of ET G to ET 0 (ET G /ET 0 ) are used to exclude the influence of climatic conditions. For various vegetation types, the exponential decay function between ET G /ET 0 and DWT are obtained as shown in Figure 9, but the parameters of these functions are different. The maximum ET G /ET 0 is 1.1 for the grass sites (the Carex stenophylla site and the Achnatherum splendens site), while the maximum ET G /ET 0 is 3.6 for the tree sites (the Populus alba site and the Salix psammophila site), resulted from the larger ET G in the tree sites. The coefficient of the independent variable (CIV) in the exponential function can reflect the decay rate of the ET G /ET 0 with DWT. For grass, ET G /ET 0 decreases more rapidly with DWT in the Carex stenophylla site (CIV = −2.35) than in the Achnatherum splendens site (CIV = −1.55), indicating that ET G is more sensitive to DWT in the Carex stenophylla site. While ET G /ET 0 decreases more rapidly with DWT in the Populus alba site (CIV = −2.66) than in the Salix psammophila site (CIV = −0.98) for the tree sites, suggesting that ET G in the Populus alba site is more sensitive to DWT.
The groundwater extinction depth is assumed to be reached when ET G /ET 0 decreases to 0.5% [49]. According to the relationship of ET G /ET 0 and DWT in our study site, the groundwater extinction depth are 4.1 m, 2.4 m, 7.1 m, and 2.9 m in the Achnatherum splendens site, Carex stenophylla site, Salix psammophila site and Populus alba site, respectively. It reflects that roots of Achnatherum splendens and Salix psammophila are deeper and better able to adapt to lower water levels.
to 0.5% [49]. According to the relationship of ETG/ET0 and DWT in o groundwater extinction depth are 4.1 m, 2.4 m, 7.1 m, and 2.9 m in the A dens site, Carex stenophylla site, Salix psammophila site and Populus alba sit reflects that roots of Achnatherum splendens and Salix psammophila are d able to adapt to lower water levels.

Conclusions
In this study, ETG of four vegetation types was estimated based on table and soil moisture fluctuations in the northeastern Mu Us sandy r to the variation of water table fluctuations at multiple wells, the influ barometric pressure could be excluded. ETG in all the observation well cant temporal and spatial variation in each vegetation site, indicating th tiveness of ETG from one observation well is poor and the mean ETG of an area will be more representative. The average ETG of multiple obse 4.01 mm/d, 6.03 mm/d, 8.96 mm/d and 12.26 mm/d at the Achnatherum sp stenophylla site, Salix psammophila site, and Populus alba site with the ave m, 0.49 m, 1.30 m, and 0.74m, respectively, in the whole growing season The main influencing factors of ETG are ET0 and DWT. The ETG is DWT in the Carex stenophylla site than that in the Achnatherum splendens s and for tree sites, ETG in the Populus alba site is more sensitive to DWT Salix psammophila site, according to the relationship between ETG/ET0 an tion, the groundwater extinction depths are also predicted to be 4.1 m, 2.9 m in the Achnatherum splendens site, Carex stenophylla site, Salix psam Populus alba site, respectively. This typical geomorphic unit in the study be able to represent the evapotranspiration characteristics of the ground vegetation covered area in the Mu Us Sandy region, northern China.

Conclusions
In this study, ET G of four vegetation types was estimated based on the diurnal water table and soil moisture fluctuations in the northeastern Mu Us sandy region. According to the variation of water table fluctuations at multiple wells, the influence of tidal and barometric pressure could be excluded. ET G in all the observation wells shows a significant temporal and spatial variation in each vegetation site, indicating that the representativeness of ET G from one observation well is poor and the mean ET G of multiple wells in an area will be more representative. The average ET G of multiple observation wells are 4.01 mm/d, 6.03 mm/d, 8.96 mm/d and 12.26 mm/d at the Achnatherum splendens site, Carex stenophylla site, Salix psammophila site, and Populus alba site with the average DWT of 1.19 m, 0.49 m, 1.30 m, and 0.74 m, respectively, in the whole growing season.
The main influencing factors of ET G are ET 0 and DWT. The ET G is more sensitive to DWT in the Carex stenophylla site than that in the Achnatherum splendens site for grass sites, and for tree sites, ET G in the Populus alba site is more sensitive to DWT than that in the Salix psammophila site, according to the relationship between ET G /ET 0 and DWT. In addition, the groundwater extinction depths are also predicted to be 4.1 m, 2.4 m, 7.1 m, and 2.9 m in the Achnatherum splendens site, Carex stenophylla site, Salix psammophila site, and Populus alba site, respectively. This typical geomorphic unit in the study area is thought to be able to represent the evapotranspiration characteristics of the groundwater-dependent vegetation covered area in the Mu Us Sandy region, northern China.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.