A New Method of Estimating Groundwater Evapotranspiration at Sub-Daily Scale Using Water Table Fluctuations

: Riparian ecosystems fundamentally depend on groundwater, and accurate estimations of groundwater evapotranspiration (ET G ) are important for understanding ecosystem functionality and managing regional water resources. Over the past several decades, various methods have been proposed to estimate groundwater evapotranspiration based on water table ﬂuctuations. However, the majority of methods cannot resolve sub-daily variations in ET G . In this study, we proposed a new hydraulic theory-based ET G estimation method at a sub-daily time scale. To evaluate its performance, we employed a variety of measurements (i.e., water table levels, latent heat ﬂux and soil water contents) at a riparian forest ( T. ramosissima ) in Northwest China from 25 July to 10 October in 2017. The results indicated that the proposed method can successfully estimate ET G at both sub-daily ( R 2 = 0.75) and daily ( R 2 = 0.88) time scales, but the variations in the speciﬁc yield under different water table conditions should be carefully taken into account. In addition, we investigated the seasonal variations in water uptake source of the riparian plant, and found that it had strong plasticity in water usage during the study period. That is, it consumed approximately equal amounts of soil water and groundwater when soil moisture was available, and tended to consume more groundwater for survival as the soil moisture was depleted. To verify the seasonal patterns of the water uptake of the riparian forest, systematic isotope-based studies are needed in future study.


Introduction
Approximate 40% of the world's landmass is made up of drylands [1], where riparian ecosystems are very important in protecting biodiversity, maintaining ecological balance, and regulating ecosystem services [2]. Many riparian plants whose roots can tap perennial water tables depend on groundwater for their survival [3]. Therefore, accurate estimate of groundwater evapotranspiration (ET G ) is critical for understanding ecosystem functionality and managing regional water resources [4]. Up to now, various methods have been proposed to quantify ET G based on diurnal water table fluctuations [5][6][7][8][9][10]. Due to their cost-effectiveness and simplicity, these ET G estimation methods have been applied across diverse ecosystems [11][12][13][14][15].
Nevertheless, there are still some insufficiencies in the developments of the ET G estimation methods using water table fluctuations. First, the majority of methods assume that groundwater recovery rate is constant during a complete one-day period, and can only calculate ET G at a daily (or longer) time scale [5,6,9]. As first noted by Troxell [16] and evidenced by numerical simulations [7], the recovery rate is not constant over time, and must be estimated as a function of time. Until now, relatively little attention has been paid to estimate the time-dependent groundwater recovery rate with the exception of that of Loheide [7] and Gribovszki et al. [8]. Second, the key parameter in the water table fluctuation methods, known as specific yield, is difficult to estimate, because neither in situ nor laboratory measurements capture its variability [4]. Previous studies have clearly documented that it not only depends upon the soil hydraulic properties and the water table depth, but is different under different water table conditions (i.e., declining, stable or rising) [17][18][19][20][21]. Unfortunately, systematic investigations of the variations in specific yield under different water table conditions are critically limited in present studies [11]. Finally, identifying the water source of plants is important for understanding the ecohydrological processes of riparian ecosystems [22,23]. Isotopic techniques have been widely used in previous studies [24][25][26][27]. However, sampling the isotopes of various water sources (i.e., soil water, groundwater and plant tissues) is technically challenging, destructive, expensive and time-consuming, which may hamper the assessment of the temporal variations in plant water source during growing seasons [28]. Thus, it is attractive to identify plant water sources by comprehensively using water table fluctuations and flux measurements [12].
Interest in ecohydrology and the need for finer resolution estimates of groundwater consumption increases are increasing. Here, we tried to develop a hydraulic theory-based method to estimate ET G at both sub-daily and daily scales by taking the time dependency of recovery rate into account. Based on 78-day continuous field observations including water table fluctuations, soil moisture and latent heat fluxes over a riparian ecosystem in Northwest China in 2017, the specific objectives of the present study were to: (1) develop a new hydraulic theory-based ET G estimation method at both sub-daily and daily scales; (2) investigate the variations in specific yield under different water table conditions, and evaluate its impacts on estimations of ET G ; and (3) identify the plant water uptake source during the whole study period.

Site Descriptions and Data Collection
The study site is located in the lower reaches of the Heihe River Basin, Gansu Province, China (101.1374 • E, 42.0012 • N; 873 m a. s. l.) ( Figure 1). Situated in the hinterland of the Asian continent, the region has a continental climate with a mean annual precipitation of 42 mm and a mean annual potential evaporation of 2241 mm [29]. The soil texture profiles are silt loam with a clay interlayer. The riparian forest community is dominated by T. ramosissima with a mean height of 1.87 m and a density of 42 stem/hm 2 [30]. The field observation system was set up as part of the Heihe Watershed Allied Telemetry Experiment Research (HiWATER) project [31,32]. Latent heat fluxes were measured using the eddy covariance (EC) system at the height of 8.0 m, and post-processed by using EdiRe software (http://www.geos.ed.ac.uk/abs/research/micromet/EdiRe, accessed on 11 August 2021). Data gaps were filled using the look-up table method, and the imbalance of energy was corrected using the Bowen ratio method [33]. Continuous complementary measurements also included standard hydro-meteorological variables, such as rainfall, air temperature, relative humidity, wind speed/direction, net solar radiation, soil temperature and moisture at different depths (i.e., 0.02, 0.04, 0.1, 0.2, 0.4, 0.8, 1.2, 1.6 and 2.0 m). These data were logged every 30 min by a digital micro-logger [32][33][34] and are available from the National Tibetan Plateau/Third Pole Environment Data Center (http://data.tpdc.ac.cn, accessed on 11 August 2021). In addition, one groundwater level monitoring well was installed near the EC flux tower. The depth of the well is about 4 m and the screen is located from the bottom of the well up to 1.0 m below surface using an 80 mm diameter PVC tubing. The groundwater level was measured using an automated pressure transducer (HOBO Water Level Logger-U20, Onset Computer Corp, Bourne, Massachusetts, USA) at half-hour intervals. The leaf area index (LAI) was not measured during the experiment period, and was extracted from MOD15A (https://modis.ornl.gov/, accessed on 11 August 2021) at 500 m spatial and 4 day temporal resolutions. An average of eight surrounding pixels around the EC tower were used to represent the seasonal variations in LAI.

Water Balance in the Riparian Zone
Because of the scarce precipitation, surface evapotranspiration (ET; L T −1 ) mainly comes from two water sources at our site: groundwater and soil water in the root layer. Therefore, total ET over the riparian zone at a daily scale could be estimated by using the water balance method: where i represents the day of the year; and ET G and ET S (all in units of (L T −1 )) the daily ET from groundwater and soil water, respectively. ET S can be calculated from the changes in soil water content in the root zone: where n is the number of the total soil layer; d l lL] the thickness of the lth soil layer (l = 1, 2, . . . , n); θ l (i) and θ l (i−1) (all in L 3 L −3 ) are the soil water content in the lth soil layer for the current day i and previous day i−1, respectively. Inspired by previous studies i.e., [5][6][7][8]35], we here proposed a new hydraulic theory-based ET G estimation method using diurnal water table fluctuations.

Hydraulic Theory-Based ET G Estimation Method
At sub-daily scale, the change in water table in the riparian zone at time t (hour) is controlled by the flow rate of groundwater from the background to the riparian zone, q(t) [L T −1 ], and the groundwater evapotranspiration rate, ET g (t) [L T −1 ]: where S y (-) is the specific yield and WT [L] is the water table depth. Here, the lateral flow is defined as positive toward the riparian zone and the water table depth is defined as zero at the land surface and becomes increasingly negative as the water table becomes deeper. Using Darcy's law and the Dupuit approximation, the groundwater flow rate toward the riparian zone can be formulated as [36]: where K s [L T −1 ] is the saturated hydraulic conductivity of the aquifer system; L [L] is the distance from the background to the riparian zone; H 0 [L] and h [L] is the groundwater elevation in the background and the riparian zone, respectively ( Figure 2); and a 1 , a 2 and a 3 are hydrological coefficients (see details in Supporting Information), which are generally time-consuming and costly to measure directly [10]. We noticed that these coefficients can be derived from the water table records during times of zero ET g . Combining Equations (3) and (4) by setting ET g = 0, we arrive at: Thus, we established a relationship between the change in water table level over time dWT/dt and water table level at sub-daily scale, and the coefficients can be easily determined by using the quadratic curve fitting method (Figure 3a). In practice, the change in water table level over time (dWT/dt) was calculated every 30 min using the slope estimated from the water table data at the time of the estimate, as well as the data before and after that time [7]. At our site, we found that observed ETs were close to zero during the time from midnight to 8 a.m. in the morning and from 7 p.m. to midnight on the day of interest. Thus, we propose using observations during these time intervals to estimate the coefficients in Equation (5). Once these coefficients are determined, the groundwater flow rate q(t) during the daytime period can be determined from the corresponding water table records using Equation (4) (Figure 3b). By rearranging Equation (3), the value of ET g at time t can then be calculated as:  Integrating Equation (6) over one-day interval ∆t, we can obtain the accumulated daily ET G : where ET G [L T −1 ] is daily ET from groundwater ( Figure 3b). Thus, we can estimate groundwater evapotranspiration at both daily [ET G ; L T −1 ] and sub-daily [ET g ; L T −1 ] scales.

Determination of Specific Yield
The specific yield (S y ) is exceedingly difficult to estimate, and is highly variable in different water table conditions [4,5,11]. According to [12], we used the information available from the water balance method to estimate S y . In brief, the value of S y was determined by minimizing the residual sum of squares (RSS) of Equation (1): where T is the number of days during the study period; ET obs (i) is EC-measured evapotranspiration on the ith day (i = 1, 2, . . . , T); ET S (i) + ET G (i; β) represents estimated evapotranspiration from soil water and groundwater on the ith day using the water balance method (Equation (1)); and β is a coefficient that equals S y when RSS is minimized. In practice, we used a Monte Carlo (MC) algorithm to find the optimal value of S y . The pseudo-code of the algorithm is given below: Step 1: Calculating daily ET S from the soil moisture observations using Equation (2); Step 2: Selecting a random value of β from a prior interval of specific yield for silt loam soil (i.e., from 0 to 0.2), and calculating ET G using Equation (7) combined with selected β and diurnal water table records; Step 3: Calculating RSS using Equation (8); Step 4: Repeating the Steps 2 and 3 for 10,000 times, and selecting the value of β with minimum RSS as the optimal estimate of S y . The best 500 values of β with minimum RSSs were also used to calculate the standard deviation of S y .
To account for the influences of water table conditions on S y , we divided the entire observations into three different periods, and the optimal value of S y for each period was determined by using the MC algorithm. The three periods respectively represented a declining, constant and rising water table condition. Hereafter, we called this calibration procedure the 'dataset-by-dataset calibration'. In addition, we applied the MC algorithm on entire observations to obtain an optimal S y value for the whole study period, called the 'whole dataset calibration'. In this way, we can investigate the variations of S y under different water table conditions, and its impact on the performances of the proposed method throughout the study period.

Evaluating the Hydraulic Theory-Based ET G Estimation Method
After the values of S y were estimated, we used five statistical measures to evaluate the performance of the proposed method. These statistical measures were the coefficient of determination (R 2 ), slope, y-intercept, bias, root-mean square error (RMSE), and relative error (RE). The calculation formulas can be found in [37]. Among them, R 2 ranges between 0 and 1, with higher values indicating a good simulation result. Slope and y-intercept indicate how well the scatter plot between observed and estimated ET fits the 1:1 line. The values of bias, RMSE and RE illustrate the difference between observed and estimated ET, with lower values indicating a better model fit.

Environmental Variables
Detailed information on the seasonality of key environmental and biological variables is essential to understand seasonal variation in plant water uptake. The seasonal change in daily net solar radiation (R n ; W m −2 ), air temperature (T a ; • C), soil water content (θ; m 3 m −3 ), and leaf area index (LAI; m 2 m −2 ) are illustrated in Figure 4. During the study period (form 25 July to 10 October), the daily mean R n varied from 33.1 to 208.7 W m −2 with an average value of 124.8 W m −2 . The variation of mean daily T a has a similar trend to R n , varying from 5.2 to 29.1 • C with an average value of around 20.1 • C (Figure 4a). The leaf area index (LAI; m 2 m −2 ) showed a declining trend during the whole study period (Figure 4b). Interestingly, soil moisture (θ; m 3 m −3 ) only in the middle layers (80-160 cm) was observed to decrease gradually during the study period, and θ in other layers (0-80 cm and 160-200 cm) varied slightly (Figure 4c). The seasonal variation in water table can be divided into three periods (Figure 4d). From 25 July to 18 August, a generally declining trend of water level was observed. After that, the water table remained at a relatively constant level. In late September, the water table began to rise due to ceased plant water uptake.

Determination of Specific Yield
The values of specific yield (S y ) estimated using the minimum RSS method are shown in Figure 5. For the dataset-by-dataset calibration procedure, it was observed that S y exhibited significant variations under different water table conditions. That is, the maximum value (0.0359) was found during the water table declining period (25 July-18 August), followed by that (0.0206) during the water table stable period (19 August-22 September), and the minimum value (0.0137) occurred during the water table rising period (23 September-10 October). For the whole dataset calibration procedure, the value of S y was estimated to be 0.0251. This value was similar to that for the water table stable period (0.0206). This may be partly explained by the fact that the dataset during the water table stable period had a comparatively larger number of observations (i.e., 34 out of 77), and subsequently gained more weight in calculating the residual sum of squares.

Evaluating the Performances of the ET G Estimation Method
Having estimated S y as described above, we can obtain groundwater evapotranspiration at both a half-hourly (ET g ) and daily (ET G ) time step using Equations (6) and (7), respectively. To evaluate the performances of our proposed method, we compared the estimated daily ET using the water balance method with EC-measured values ( Figure 6). The results indicated that the estimated daily ET using S y calibrated by different datasets agreed well with EC-measured values. The points in the plots of measured-versus-estimated daily ET fell tightly along the 1:1 line (slope = 0.91, intercept of 0.31 mm day −1 and a correlation coefficient of 0.88), and the estimated daily ET fluctuated tightly with the measured values (bias= 0.06, RMSE = 0.85 and RE = 0.20) (Figure 6a,c). On the contrary, less satisfactory results were obtained by using S y calibrated by the whole dataset, with relatively low values of slope (0.66) and R 2 (0.76) and high RMSE (1.21) and RE (0.28) ( (Figure 6b). In general, the estimated daily ET G using S y calibrated by the whole dataset was slightly underestimated and overestimated during the water table declining period and rising period (Figure 6c), respectively. This seems to be due to the significant differences in estimates of S y between the dataset-by-dataset and multi-dataset procedures during these two periods.  Furthermore, we evaluated the performances of our proposed method at a half-hourly time step. Noticeably, the soil water evapotranspiration component in the water balance equation was not available at a half-hourly time step because the hour-to-hour variations in soil water content were generally not detectable. Thus, the slope of linear regression between EC-measured ET and estimated ET g at a half-hourly time step was expected to be less than 1 (Figure 7). Similar to the results observed at a daily time step, the estimated ET g using S y calibrated by the dataset-by-dataset procedure showed better correlations with measured ET than that using S y calibrated by the whole dataset (Figure 7a,b). In addition, it was observed that the diurnal pattern of estimated ET g and measured ET are very similar under different water conditions, despite that the magnitudes of ET g were different for the two S y calibration procedures (Figure 7c). Thus, it seemed that the method proposed here can successfully estimate groundwater evapotranspiration at different time steps, but the variations in S y should be carefully taken into account.

Water Use Pattern
To reveal the water uptake pattern of the riparian plants, we plotted the seasonal variations in ET and its different components (i.e., ET S and ET G estimated using S y calibrated by different datasets) in Figure 8. The results indicated that during the whole study period the cumulative ET S and ET G were 121 mm and 202 mm (Figure 8b), respectively. The ratio of ET that comes from groundwater (ET G /ET) and soil water (ET S /ET) was 0.62 and 0.38, respectively. Thus, groundwater was the main source for plant water consumption during the whole study period. Interestingly, we observed that the riparian plant had strong plasticity in water uptake during the whole study period. When soil moisture is available (i.e., from 25 July to 12 August), the plant consumed approximately equal amounts of soil water and groundwater resources (Figure 8a), and the ratios of ET G /ET and ET S /ET were 0.54 and 0.46 (Figure 8b), respectively. As the soil moisture was depleted (i.e., after 13 August), the plant tended to consume more groundwater for survival (Figure 8a), and the ratios of ET G /ET (>0.64) were larger than that of ET S /ET (<0.36) (Figure 8b).

Discussion
The specific yield (S y ) is of crucial importance for estimating ET G based on water table fluctuations because its error is translated directly to the final estimates [5,10,38]. In this study, the values of S y estimated using the minimum RSS method fell within the ranges of the literature reports for silt loam soil. For example, Johnson [39] reported that S y for silt loam varied from 0.01 to 0.2. Singhal and Gupta [40] documented that the values of S y for rocky silt loam ranged from 0.02 to 0.05. In addition, we also calculated the S y values using the method proposed by Loheide et al. [21], which ranged from 0.0297 to 0.0353 with a mean of 0.0339 during the whole study period (see Supporting Information Figure S1). Thus, we thought that the minimum RSS method was suitable to obtain proper S y values for ET G estimates in our site. Interestingly, we found that S y varied significantly under different water table conditions, and its value for a declining water table condition was about three times that of a rising water table condition. This can be explained by the fact that encapsulated air in the aquifer, which can be as high as 20% of soil porosity [41], reduces the value of S y during the water table rising periods. Therefore, it is important to take the variations in S y under different water table conditions into account for long-term groundwater evapotranspiration estimations.
We further compared the performance of our proposed method with the methods developed by Loheide [7] and Griovszki et al. [8]. The Gribovszki method uses the rate of the water table change dWT/dt [L T −1 ] to estimate the groundwater recharge rate q(t). For each day, and maximum and minimum q(t) were obtained by selecting the largest positive rate of dWT/dt and the mean of dWT/dt between midnight and 6 a.m., respectively; then it is assumed that q(t) behaves linearly between two consecutive estimations. Loheide [7] proposed to estimate q(t) as a function of the detrended water table. The function, set up using data between midnight to 6 a.m. of two subsequent days, is assumed to be approximately linear during each day. Overall, the performance of our method was very similar to the Gribovszki method ( Figure 9). However, we found that during the nighttime period (i.e., from 9 p.m. of the previous night to 6 a.m. in the morning) the estimated q(t) by the Gribovszki method was generally higher than the rate of dWT/dt (Figure 9a). Thus, unreasonable positive ET g was obtained by this method during the nighttime (Figure 9b). In addition, the maximum q(t) estimated by our method was generally higher and occurred earlier than that estimated by the Gribovszki method. Thus, the values of ET g estimated by our method were slightly higher than that of the Gribovszki method in the afternoon. On the contrary, the method proposed by Loheide [7] performed relatively unsatisfactorily. This may be mainly because the assumptions of the Loheide method (i.e., the head at the recovery source is constant or follows the general trend of the water table) may not be met at our site. Fahle and Dietrich [10] have systematically compared the performances of six widely used groundwater evapotranspiration estimation methods, and also found that the Loheide method performed considerably worse than other methods.
The results of our study revealed that the riparian plant exhibit substantial water uptake plasticity during different growing stages. When soil moisture was adequate, the riparian plant mainly used soil water. As soil water depleted, it then tended to depend on groundwater ( Figure 8). Similar water use strategies were also observed in other arid and semi-arid plant communities [24][25][26][27][42][43][44]. Furthermore, the estimated ratios of ET G /ET by our method agreed with the results using the stable isotopic methods. For example, Wu et al. [43] documented that the ratios of ET G /ET of T. ramosissima at the southern edge of the Gurbantonggut desert were more than 0.43 in summer and autumn. Zhao et al. [45] reported that the ratio of plant water uptake that from groundwater of T. ramosissima in the lower reaches of Heihe River Basin was 0.49 in wet seasons and increased to 0.90 in dry seasons. To assess the seasonal patterns of water uptake of the riparian plant, a systematic collection of stable isotopes of different water sources (i.e., soil water, groundwater and xylem water) is needed in the future studies.
Finally, it should be noticed that there are still some uncertainties in our method. First, we used information from the water balance method to help estimate the values of S y . In this study, the profile of soil moisture was measured only at one site, and it may be difficult to capture the spatial variability in soil moisture. Thus, there may be considerable uncertainty in the soil moisture portion of the water balance equation. To overcome this problem, a soil moisture monitoring network is planned to be set up at our site. Second, the periods from midnight to 8 a.m. in the morning and from 7 p.m. to midnight on the day of interest were selected for the recovery analysis at our site. However, this is a somewhat subjective choice of time period. To apply this method at other sites, a slightly different time period may be better for calibrating the parameters in Equation (5). Finally, if hydraulic redistribution is a significant process, plant water uptake from groundwater during the night is not zero. Thus, the fundamental assumption of diurnal water table fluctuation based methods will be violated [46].

Conclusions
In this study, we developed a new hydraulic theory-based method to estimate groundwater evapotranspiration at a sub-daily time scale. The diurnal pattern of estimated groundwater evapotranspiration exhibited similar trends to EC-measured evapotranspiration. Thus, it seemed that this method permits improvement in the estimation of groundwater evapotranspiration just by using diurnal water table fluctuations. Another important problem in applying water table fluctuation methods is to properly determine the values of specific yield (S y ). Our study revealed that the S y values determined by the minimum residual sum of squares method were consistent with those reported in the literature for similar soil, and varied considerably under different water table conditions. For long-term groundwater evapotranspiration estimations, the variations of S y under different water table conditions (i.e., falling, stable and rising) should be properly taken into account. In the future, stable isotope-based studies are needed to evaluate the accuracy of the method in estimating groundwater evapotranspiration rates.

Data Availability Statement:
The MOD15A products used in this study are available at https: //modis.ornl.gov/. The flux data are available from the National Tibetan Plateau/Third Pole Environment Data Center (http://data.tpdc.ac.cn).