Using High-Frequency Water Vapor Isotopic Measurements as a Novel Method to Partition Daily Evapotranspiration in an Oak Woodland

: Partitioning evapotranspiration ( ET ) into its constituent ﬂuxes (transpiration ( T ) and evaporation ( E )) is important for understanding water use e ﬃ ciency in forests and other ecosystems. Recent advancements in cavity ringdown spectrometers (CRDS) have made collecting high-resolution water isotope data possible in remote locations, but this technology has rarely been utilized for partitioning ET in forests and other natural systems. To understand how the CRDS can be integrated with more traditional techniques, we combined stable isotope, eddy covariance, and sap ﬂux techniques to partition ET in an oak woodland using continuous water vapor CRDS measurements and monthly soil and twig samples processed using isotope ratio mass spectrometry (IRMS). Furthermore, we wanted to compare the e ﬃ cacy of δ 2 H versus δ 18 O within the stable isotope method for partitioning ET. We determined that average daytime vapor pressure deﬁcit and soil moisture could successfully predict the relative isotopic compositions of soil ( δ e ) and xylem ( δ t ) water, respectively. Contrary to past studies, δ 2 H and δ 18 O performed similarly, indicating CRDS can increase the utility of δ 18 O in stable isotope studies. However, we found a 41–49% overestimation of the contribution of T to ET ( f T ) when utilizing the stable isotope technique compared to traditional techniques (reduced to 4–12% when corrected for bias), suggesting there may be a systematic bias to the Craig-Gordon Model in natural systems.


Introduction
Evapotranspiration (ET) is an important process in the terrestrial hydrological cycle that accounts for evaporation (E) from the soil, open water, and canopy-intercepted water and transpired (T) water from vegetation. ET constitutes a large percentage of the water cycle in most environments, and up to 95% in arid environments, deeming it an important water flux at a variety of spatial scales [1,2]. An improved understanding of ET partitioning and the contribution of T to ET (f T ) is necessary for refined water resource management and quantification of vegetative responses to climate change and alterations to the carbon cycle. For instance, complex dynamics exist between carbon fluxes and precipitation at multiple scales, and the water-use efficiency of plants varies with ecosystem aridity and vegetation cover [3,4]. Partitioning of ET is necessary because T is seen as a desirable aspect of the water cycle that allows vegetation to grow, while water that is evaporated is generally seen as being "lost" from the system. days where soil and twig sampling provided proxy values for the isotopic compositions of soil (δ e ) and xylem (δ t ) water; (2) to create δ e and δ t predictive models for partitioning daily ET on non-sampling days using a combination of monthly field samples, environmental, and micrometeorological data; (3) to compare the efficacy of δ 2 H and δ 18 O in partitioning ET.

Site Description
The study site is located within management unit 75 of Lyndon B. Johnson (LBJ) National Grasslands in Wise County, TX, just north of the city of Decatur, TX. It is a core terrestrial site in the National Ecological Observatory Network's (NEON) Southern Plains domain (D11). LBJ is considered part of the Cross Timbers ecoregion of Texas, and wooded areas consist mostly of mixed oaks (Quercus spp.), elms (Ulmus crassifolia and Ulmus alata), eastern red cedar (Juniperus virginiana), and sugarberry (Celtis laevigata) with a mean canopy height of 3.96 m. The LBJ study site is situated at approximately 250 m above sea level, and the EC tower is located precisely at 33.40123 • N, −97.57 • W. Dominant tree species in the immediate vicinity of the eddy covariance tower are blackjack oak (Quercus marilandica) and post oak (Quercus stellata).
Soils consist of Keeter very fine sandy loams (1-6 percent slopes) and Duffau-Weatherford Complex (3-8 percent slopes) [22]. The mean annual temperature is 18 • C and peak temperatures usually occur in August, while mean annual precipitation (MAP) is 840 mm·year −1 with a bimodal distribution (rainy seasons in late spring and fall) [23]. Mean annual potential evapotranspiration (PET) and runoff are 1785 mm·year −1 and 8 mm·year −1 , respectively [24]. The site has undergone prescribed burning for two consecutive years (2017 and 2018) during spring as a strategy to control the thick understory of greenbrier (Smilax rotundifolia), though the site was not burned at any time during the study period from May 2019 to January 2020. The total period can be divided into the growing/vegetation season (May to the end of October) and the dormancy period (November to January).

Micrometeorology
At LBJ, NEON operates an EC tower that is fitted with instruments at four different heights in the woodland canopy, with one set of sensors extending above the canopy as well. These instruments include a sonic anemometer, net radiometer, open-path CO 2 /H 2 O analyzer, as well as tubing that collects vapor at each level and directs it inside a sheltered unit for isotopic analysis on the Picarro™ CRDS.
Actual evapotranspiration (ET EC ) values were converted from the given latent heat flux (λE, W·m −2 ) values to mass units (mm·day −1 ) using the latent heat of vaporization of water [25] for every 30 min and then summed over 24 h. Days where there were not at least 40/48 data points available were excluded from the analysis.
The flux footprint area was modeled for May of 2019 using the Flux Footprint Prediction Model (FFP) (found at https://geography.swansea.ac.uk/nkljun/ffp/www/) [26]. Input data for the model include wind speed, wind direction, friction velocity, Obukhov length [27], displacement height, measurement height, and standard deviation of lateral velocity fluctuations after rotation.
Due to EC sensor issues at the beginning of the study, the Penman-Monteith reference evapotranspiration (ET 0 ; mm·day −1 ) was calculated to provide ET estimates. The FAO-56 method [28] was utilized at a daily scale. If there were incomplete data (<40 data points per day) for any of the input parameters, then that day was excluded from the analysis. A summary of NEON data products used in this study, along with their associated instrumentation, can be found in Table S1 [23].

Sap Flow Sensors
To directly measure transpiration in the EC tower footprint, eleven oak trees were fitted with thermal dissipation sap flow sensors [14,15] on 11 May 2019 and ran continually until 31 January 2020. All instrumented trees were either Quercus marilandica or Quercus stellata with a minimum, Water 2020, 12, 2967 4 of 19 maximum, and average diameter at breast height (DBH) of 21.5, 54.5, and 32.6 cm, respectively (Table S2). All trees had branches above sensor installation height (~1.3 m) and were not shaded by any nearby tree. Sap flow sensors were custom made in the Moore Ecohydrology Lab at Texas A&M University [14,29]. Data was collected every 30 s and later averaged over 30 min intervals on a CR1000 datalogger (CR1000, Campbell Scientific Inc., Logan, Utah).
The sapwood depth was determined by taking tree cores using an increment borer and partially immersing the fresh cores in a safranin-fucsin dye [30][31][32]. Where the dye clearly moved from the bottom to the top of the core was considered to be actively conducting sapwood. All trees with sensors had a sapwood radius greater than the sensor depth of 2 cm (Table S2) [33]. Sapwood area ranged from 0.022 to 0.149 m 2 , with an average of 0.056 m 2 per tree.

Baseline Corrections
It was determined after initial calculations of VPD and sap flow processing that nighttime flows were occurring at the site. VPD did not reach zero on a nightly basis, and there was not a consistent baseline for ∆T m which introduces error into the original thermal dissipation equation. To account for this error and adjust ∆T m for nighttime flows, a "baseline smoothing" process was utilized where daily minimum VPD was plotted with daily precipitation ( Figure 1A) and days that had <0.1 kPa VPD and zero precipitation were considered days where calculated ∆T m was accurate. The ∆T m for remaining days were linearly gap-filled to create a smooth baseline. These new ∆T m were then used to re-process the sap flow data.
Water 2020, 12, x FOR PEER REVIEW 4 of 20 S2). All trees had branches above sensor installation height (~1.3 m) and were not shaded by any nearby tree. Sap flow sensors were custom made in the Moore Ecohydrology Lab at Texas A&M University [14,29]. Data was collected every 30 s and later averaged over 30 min intervals on a CR1000 datalogger (CR1000, Campbell Scientific Inc., Logan, Utah). The sapwood depth was determined by taking tree cores using an increment borer and partially immersing the fresh cores in a safranin-fucsin dye [30][31][32]. Where the dye clearly moved from the bottom to the top of the core was considered to be actively conducting sapwood. All trees with sensors had a sapwood radius greater than the sensor depth of 2 cm (Table S2) [33]. Sapwood area ranged from 0.022 to 0.149 m 2 , with an average of 0.056 m 2 per tree.

Baseline Corrections
It was determined after initial calculations of VPD and sap flow processing that nighttime flows were occurring at the site. VPD did not reach zero on a nightly basis, and there was not a consistent baseline for ∆ which introduces error into the original thermal dissipation equation. To account for this error and adjust ∆ for nighttime flows, a "baseline smoothing" process was utilized where daily minimum VPD was plotted with daily precipitation ( Figure 1A) and days that had <0.1 kPa VPD and zero precipitation were considered days where calculated ∆ was accurate. The ∆ for remaining days were linearly gap-filled to create a smooth baseline. These new ∆ were then used to re-process the sap flow data.

Transpiration Scaling
Scaling up sap flow measurements to the footprint of the EC tower was done using surveys from trees in 29 individual 400 m 2 vegetation plots [23]. Within the 29 plots, there were 915 live individuals total for which the stem diameter and tree height were recorded. Using DBH, the basal area was then calculated for each individual tree, and an allometric equation was applied to calculate the sapwood area for each individual based on the basal area (y = 0.6075x + 0.0004, R 2 = 0.92). The total sapwood area per plot was divided by the total plot area to get sapwood area per unit ground area (m 2 m −2 ) for each plot, and a threshold value of 0.2 m 2 of sapwood area was used to differentiate between woodland and grassland plots ( Figure S1). The average sapwood area per unit ground area for all of the woodland plots (x = 0.0012 m 2 m −2 ± 0.0006, n = 19) was used to calculate the total sapwood area in the EC tower airshed. Transpiration (mm·day −1 ) was then calculated by multiplying the average sap-flux density of all sap flow trees (kg m −2 day −1 ) by the average sapwood area per unit ground area. This bottom-up approach has been utilized in other studies [5,8] and validated by a review on methods for scaling up T measurements from tree to stand level [34]. It should be noted that this T estimate was based on oak allometry and sap flux, but there were other species present in the plots that may have different sapwood area to basal area ratios. Note that no understory species were included in this T estimate, including any woody vegetation with DBH < 0.5 cm.

Method Framework
Given ET = E + T, the respective contributions of E and T to ET can be calculated using the isotopic mass balance approach, assuming a two-source model where δ ET , δ E , and δ T are equal to the isotopic composition of evapotranspiration, soil evaporation, and plant transpiration, respectively. The ratio of T to ET can then be calculated as To use this method, the isotopic composition of water in three "end members" (i.e., soil, vegetation, and atmosphere) are needed, and the resulting ratios of each can be applied to total ET to get values of the individual components [35,36]. Differences in the isotopic composition of water in vegetation and in soil occur from fractionation during the evaporation process. The soil evaporating front changes with soil texture, soil moisture, and other environmental conditions, but is generally at the surface when the soil is saturated, and between 0.2-0.3 m in depth when the soil is not saturated [37][38][39]. Water is not fractionated during the transpiration process under steady-state conditions, so the isotopic composition of transpired water can be assumed equal to that of the xylem water [18,40,41]. Steady-state conditions are usually met around midday in field conditions [42]. Figure 2 shows the stepwise process for partitioning ET using the stable isotope method.

δET, δE, and δT
The Keeling plot approach was used for obtaining δET [35]. Measurements of the concentration and isotopic composition of water vapor at varying heights in the woodland canopy are sufficient to fit the linear model used in the Keeling plot approach [43][44][45]. Data for this approach were available from NEON [23]. Water vapor was collected by tubing at five different heights on a 30 min time step in the woodland canopy and circulated into the instrument hut where it was then analyzed for δ 18 [35]. Measurements of the concentration and isotopic composition of water vapor at varying heights in the woodland canopy are sufficient to fit the linear model used in the Keeling plot approach [43][44][45]. Data for this approach were available from NEON [23]. Water vapor was collected by tubing at five different heights on a 30 min time step in the woodland canopy and circulated into the instrument hut where it was then analyzed for δ 18 O and δ 2 H via cavity ringdown spectroscopy (Picarro L2130-i, Picarro Inc., Santa Clara, CA, USA) in near-real-time. Water vapor concentrations at the same heights were included in the bundled EC data products. Only values from 10:30-14:00 h were used in the analysis.
Values of δ E were attained from a combination of field samples and the Craig-Gordon Model (CGM) [46] where δ e is the isotopic composition of liquid water at the soil evaporating front, α * is the equilibrium fractionation factor (1.0098 and 1.084 at 20 • C for 18 O and 2 H, respectively [47]), h is relative humidity, δ v is the isotopic composition of the background atmospheric water vapor, α k is the isotopic fractionation factor (0.9755 and 0.9723 for 2 H and 18 O, respectively [48,49]), ε eq is equal to 1000(1 − 1/α*), and ε k is equal to 1000 (αk − 1). Values for δ T were attained directly through a sampling of twig water under steady-state assumptions, as has been done in other studies [35]. Under this scenario, it is assumed that the isotopic composition of the source water and transpired water are the same (δ X = δ T ).

Soil and Twig Sampling
To obtain isotopic values for soil water at the evaporating front (δ e ) and xylem water (δ t ), monthly soil and twig samples were taken. There was a total of nine sampling events during the study period that were carried out approximately every four weeks (Day of Year [DOY] 143, 170, 204, 234, 266, 294, 329, 351 in 2019 and DOY 24 in 2020). Soil samples (n = 4 per sampling day) were collected at midday (between 10:30-14:00 h) every month during the period of the study (May 2019-January 2020) to capture seasonal variability in δ E . Samples were collected at 0.2 to 0.3 m depth when the soil was unsaturated, and taken from the surface when the soil was saturated. Soil saturation was determined at the site on sampling days before cores were taken. Samples were taken using a hand auger and immediately transferred into 12 mL glass vials and sealed with a cap and parafilm. To prevent evaporation, vials were transported upside down at~34 • C to a freezer within 8 h. Samples were kept frozen until being processed at the Stable Isotopes for Biosphere Science (SIBS) Laboratory at Texas A&M University.
Twig samples (n = 4 per sampling day) were collected on the same days as soil samples during the same time window to satisfy steady-state conditions for xylem water sampling and to capture seasonal variability in δ T . Samples were taken using a pole pruner from mature branches on trees containing sap flow sensors and were approximately three to five centimeters in length and less than 19 mm in diameter. The bark was then stripped from the twigs, and the twigs were immediately transferred into 12 mL glass vials and stored in the same manner as the soil samples.

Cryogenic Water Extractions and Analysis
All samples underwent water extractions using a cryogenic (liquid nitrogen) vacuum distillation system. Soil samples were extracted for approximately 75 min and twigs for approximately 55 min. The extracted water was then transferred to scintillation vials via pipet and stored at 34 • C until analysis. The hydrogen and oxygen isotope analyses were performed using a Thermo Scientific High Temperature Conversion/Elemental Analyzer (TC/EA) coupled to a Conflo IV and a Thermo Scientific Delta V Advantage IRMS at the Stable Isotopes for Biosphere Science Laboratory, Texas A&M University. One micro liter (1 µL) of the sample was injected into TC/EA using a PAL autosampler. The injected sample was converted to H 2 and CO gas by pyrolysis reaction through a glassy carbon tube filled with glassy carbon chips and heated at 1370 • C. The H 2 and CO gas were separated by a 2 m packed gas chromatograph and were analyzed for the hydrogen and oxygen isotope ratios, respectively, in the Delta V Advantage IRMS. One sample was injected three times, and reported values are an average of these triplicate injections ± 0.01% .

f T Analysis for Sampling Days
To determine how f T compares between the stable isotope method and the residual of ET (ET EC or ET 0 ) and T, only data from days when the soil and twig samples were collected (to provide surrogate values for the Keeling plots and CGM) were analyzed initially for the growing season (May-October 2019). Four methods were used to calculate f T : (1) T/ET EC , (2) T/ET 0 , (3) δ 2 H-based stable isotope method, and (4) δ 18 O-based stable isotope method. This approach was used to provide insight as to how various partitioning methods compare to each other, how δ 2 H and δ 18 O compare within the stable isotope method, as well as how the various methods might change over time and with changes in seasonality. Samples were taken for November 2019 to January 2020 as well; however, it was assumed that f T would be zero during this time, due to the trees having entered dormancy and the leaves abscising.

Predictive Model Development
An advantage of our study relative to others like it is the high data resolution and data availability across spatial and temporal scales that are collected by NEON. Due to the high-resolution data that is available throughout the growing season, two models were developed for predicting δ e and δ t for all days in the study period based on data from the sampling days. As obtaining soil and twig samples on a daily basis is cumbersome, expensive, and logistically infeasible, this inherently limits the utility of the stable isotope method by limiting the temporal resolution of surrogate values used in the CGM and Keeling plots. However, this barrier may be removed by predicting the isotopic composition of δ e and δ t (both δ 2 H and δ 18 O) based on other environmental variables that are easier to measure and do not require on-site presence.
Single and multiple linear regression models were tested using environmental data provided from NEON (air temperature, VPD, soil moisture, precipitation, soil heat flux, net radiation, water vapor concentration, water vapor isotopic composition, wind speed, and relative humidity) to predict the measured isotopic composition of δ e and δ t for sampling days during the growing season. Values from the δ t predictive models were used directly for δ T in the isotope method for partitioning ET, and values from the δ e predictive models were used in the CGM to get δ E. f T was then determined every day from 11 May to 31 October for both δ 2 H and δ 18 O.

Environmental Data
The majority of the growing season was characterized by high air temperature (>20 • C), stable wind speed, high daytime VPD (>0.3 kPa), and low soil moisture (<0.20 cm 3 cm −3 ) following storms in May and June (Figure 3). Soil moisture fell consistently between July and September, and this same period experienced the highest daily air temperatures and the least amount of precipitation during the study period. Daytime VPD also reached its maxima in August, indicating peak summer conditions.

Environmental Data
The majority of the growing season was characterized by high air temperature (>20°C), stable wind speed, high daytime VPD (>0.3 kPa), and low soil moisture (<0.20 cm 3 cm −3 ) following storms in May and June ( Figure 3). Soil moisture fell consistently between July and September, and this same period experienced the highest daily air temperatures and the least amount of precipitation during the study period. Daytime VPD also reached its maxima in August, indicating peak summer conditions.  However, with the onset of the dormant season at the end of October and the beginning of November, the air temperature dropped markedly, and daily average wind speed showed higher fluctuations without the tree canopy acting as a buffer around the EC tower ( Figure 3A). Due to changes with wind speed and air temperature dynamics, daytime VPD also began to drop, relative to growing season conditions ( Figure 3B). The beginning of the dormant season was also characterized by a few large rain events that elevated soil moisture levels relative to the summer, and smaller, more frequent rain events maintained these higher levels throughout the end of the study ( Figure 3C).
Water vapor isotopes also had distinct characteristics during the growing and dormant seasons. While trees were still transpiring, temperatures were still warm, and wind speeds were stable, average water vapor δ 2 H and δ 18 O was also relatively stable ( Figure 4). As trees began to senesce in October and November, resulting in less stomatal conductance, the boundary layer conditions changed, and this is reflected in the higher degree of fluctuations of water vapor isotopes in the dormant season relative to the growing season.
Water vapor isotopes also had distinct characteristics during the growing and dormant seasons. While trees were still transpiring, temperatures were still warm, and wind speeds were stable, average water vapor δ 2 H and δ 18 O was also relatively stable ( Figure 4). As trees began to senesce in October and November, resulting in less stomatal conductance, the boundary layer conditions changed, and this is reflected in the higher degree of fluctuations of water vapor isotopes in the dormant season relative to the growing season.

Water Fluxes
Due to technical issues, actual evapotranspiration (ETEC) was unreliable for the beginning of the study period ( Figure 5). However, there was a moderate agreement for daily evapotranspiration estimated by Penman-Monteith (ET0) and the eddy covariance method for the remainder of the growing season (2.25 ± 0.89 mm day −1 and 2.60 ± 0.92 mm day −1 , respectively). The eddy covariance method provided sufficient data 40% of the time from 11 May 2019 to 31 January 2020 (107/265 days in the study period), and ET0 provided sufficient data 88% of the time (233/265 days).

Water Fluxes
Due to technical issues, actual evapotranspiration (ET EC ) was unreliable for the beginning of the study period ( Figure 5). However, there was a moderate agreement for daily evapotranspiration estimated by Penman-Monteith (ET 0 ) and the eddy covariance method for the remainder of the growing season (2.25 ± 0.89 mm·day −1 and 2.60 ± 0.92 mm·day −1 , respectively). The eddy covariance method provided sufficient data 40% of the time from 11 May 2019 to 31 January 2020 (107/265 days in the study period), and ET 0 provided sufficient data 88% of the time (233/265 days).  Total precipitation during the measurement period prior to dormancy (11 May to 31 October 2019) was 315.77 mm, or about 38% of MAP, which was 36% below normal for that time of year based on data from the last 50 years [50]. Total T, ET 0 , and ET EC for this same time period were 102.2 mm (12% of MAP), 391.5 mm (47% of MAP), and 452.4 mm (54% of MAP), respectively, for that partial growing season. Average transpiration for the study period was 0.59 ± 0.18 mm·day −1 , and trees were officially deemed dormant on DOY 325 (21 November 2019) based on sap flow data, just before the November sampling event. The average ET 0 and ET EC for the dormant season (November-January) was 0.94 ± 0.55 and 0.44 ± 0.28 mm·day −1 .

Global Meteoric Water Line and Keeling Plot
The isotopic composition (δ 18 O and δ 2 H) of atmospheric and precipitation samples fell closely along the global meteoric water line (GMWL), as was expected ( Figure 6). The soil samples fell to the right of the GMWL, indicating preferential evaporation of lighter isotopologues of water and soil enrichment with the heavier isotopologues. These results indicate the evaporative fractionation of soil water and an appropriate degree of separation between the signals of twig and soil water isotope ratios. The average of the inverse water concentration between 10:30-14:00 h was used against the isotope ratio (both δ 18 O and δ 2 H) during the same time period to construct turbulent mixing relationships, or Keeling plots. In examining Keeling Plot relationships, we found that the inverse water vapor concentration was closely related to isotopic ratios during all sampling events (p < 0.05) except July (δ 2 H and δ 18 O) and October (δ 18 O) with R 2 ranging up to 0.86 (Table 1). As 1/[H2O] increased, heavier isotopes tended to decrease. The monthly predictions of δET from Keeling plots varied widely throughout the study period for both H and O.  The average of the inverse water concentration between 10:30-14:00 h was used against the isotope ratio (both δ 18 O and δ 2 H) during the same time period to construct turbulent mixing relationships, or Keeling plots. In examining Keeling Plot relationships, we found that the inverse water vapor concentration was closely related to isotopic ratios during all sampling events (p < 0.05) except July (δ 2 H and δ 18 O) and October (δ 18 O) with R 2 ranging up to 0.86 (

f T for Soil and Twig Sampling Days
The average f T calculated from the stable isotope method with data from only sampling days was 0.85 ± 0.27 for δ 18 O and 0.7 ± 0.55 for δ 2 H. On two of the sampling days where the Keeling Plots had significant slopes (May and September) f T exceeded 100% of total ET for both δ 2 H and δ 18 O. Additionally, the δ 2 H July sampling event had a f T < 0, though the Keeling Plot for this event was not significant. The average f T calculated from ET 0 was substantially lower (0.25 ± 0.1) and never went below 0 or above 1 (Figure 7). Additionally, the δ 2 H July sampling event had a fT < 0, though the Keeling Plot for this event was not significant. The average fT calculated from ET0 was substantially lower (0.25 ± 0.1) and never went below 0 or above 1 (Figure 7).

Predictive Models for Growing Season and Dormant Season Sampling Days
It was determined that average daytime VPD (between 6:00-18:00 h) and average soil moisture at 0.16 m depth were the single best predictors of δe and δt, respectively (Figure 8). Soil moisture at other depths was tested; however, 0.16 m performed the best and was empirically logical as it is the depth in-between saturated (surface) and unsaturated (0.2 to 0.3 m) sample depths. Equations from Figure 8 were then used to estimate daily values of δe and δt for use in the stable isotope method to partition ET. All growing season models were significant at the 95% confidence level, with the exception of the δe δ 2 H model. A total of 6 data points were used in the δt models, but the July sampling event was excluded from the δe models as its Keeling Plot had an insignificant slope (Table  1), and its isotope ratio values were outliers compared to the other five months.

Predictive Models for Growing Season and Dormant Season Sampling Days
It was determined that average daytime VPD (between 6:00-18:00 h) and average soil moisture at 0.16 m depth were the single best predictors of δ e and δ t , respectively (Figure 8). Soil moisture at other depths was tested; however, 0.16 m performed the best and was empirically logical as it is the depth in-between saturated (surface) and unsaturated (0.2 to 0.3 m) sample depths. Equations from Figure 8 were then used to estimate daily values of δ e and δ t for use in the stable isotope method to partition ET. All growing season models were significant at the 95% confidence level, with the exception of the δ e δ 2 H model. A total of 6 data points were used in the δ t models, but the July sampling event was excluded from the δ e models as its Keeling Plot had an insignificant slope (Table 1), and its isotope ratio values were outliers compared to the other five months. For the dormant season, the same approach was utilized for δe and δt modeling during the growing season. The additional data from November 2019-January 2020 was added to the four models, and new regressions were conducted. The δt models were mostly unaffected and continued to be significant for both δ 2 H and δ 18 O (Figure 9). However, for the δe models, the slope of the regression line switched directions, the R 2 values dropped markedly, and the p-value for the previously significant δe δ 18 O model jumped to 0.75. Due to these reasons, δe and δt were not modeled for the dormant season.  For the dormant season, the same approach was utilized for δ e and δ t modeling during the growing season. The additional data from November 2019-January 2020 was added to the four models, and new regressions were conducted. The δ t models were mostly unaffected and continued to be significant for both δ 2 H and δ 18 O (Figure 9). However, for the δ e models, the slope of the regression line switched directions, the R 2 values dropped markedly, and the p-value for the previously significant δ e δ 18 O model jumped to 0.75. Due to these reasons, δ e and δ t were not modeled for the dormant season. For the dormant season, the same approach was utilized for δe and δt modeling during the growing season. The additional data from November 2019-January 2020 was added to the four models, and new regressions were conducted. The δt models were mostly unaffected and continued to be significant for both δ 2 H and δ 18 O (Figure 9). However, for the δe models, the slope of the regression line switched directions, the R 2 values dropped markedly, and the p-value for the previously significant δe δ 18 O model jumped to 0.75. Due to these reasons, δe and δt were not modeled for the dormant season.

Keeling Plot
To predict δ ET for the entire growing season, atmospheric isotope data and water vapor concentration data for every day from 11 May 2019 to 31 October 2019 were collected and filtered for data only during 10:30-14:00 for each of the five heights in EC tower. This data was then input into Rstudio and underwent linear regression (Keeling plots/turbulent mixing relationships) for every day there were sufficient data, for both δ 18 O and δ 2 H. Relevant parameters were extracted for each day, including sample size (N), slope, intercept (δ ET ), R 2 , and a p-value, which are reported in Table S3. Of the 164 days in the growing season, there were sufficient EC tower data to construct Keeling plots for 131 days (80%). Of these 131 days with complete data, there were significant slopes for 112 days (86%) using δ 2 H, and 89 days (68%) using δ 18 O. Keeling plots were not constructed for the dormant season, due to issues with δ e modeling and unstable boundary layer conditions.

Craig-Gordon Model
Of the 164 days in the growing season, the CGM successfully ran for 148 (90%) of days. In this case, δ 18 O and δ 2 H performed equivalently in terms of producing values for δ E , and the limiting factor for running the CGM was obtaining complete water vapor isotope data from the EC tower and CRDS. Sixteen days (10%) were missing or had incomplete water vapor isotope data. Relative humidity data were complete for the entire season, as well as soil moisture used to predict δ e .

Daily f T for Entire Study Period
Days with insignificant Keeling plots (p-value > 0.05) were excluded for f T using the stable isotope method, along with days having an f T < 0 or >1.0 for all methods. Day 131 was excluded, due to site setup and disturbance. For the first approach (T/ET EC ), the average f T was 0.36 ± 0.31, and there was a total of 70 days (43%) with valid and complete data. For the second approach (T/ET 0 ), the average f T was 0.28 ± 0.11, and there was a total of 140 days (86%) with valid and complete data. Both δ 2 H and δ 18 O had an average f T of 0.77 ± 0.20, though there were only valid data for 79 and 56 of the 164 days (48% and 34%), respectively. The results from raw calculated f T are reported in Figure 10. To predict δET for the entire growing season, atmospheric isotope data and water vapor concentration data for every day from 11 May 2019 to 31 October 2019 were collected and filtered for data only during 10:30-14:00 for each of the five heights in EC tower. This data was then input into Rstudio and underwent linear regression (Keeling plots/turbulent mixing relationships) for every day there were sufficient data, for both δ 18 O and δ 2 H. Relevant parameters were extracted for each day, including sample size (N), slope, intercept (δET), R 2 , and a p-value, which are reported in Table  S3. Of the 164 days in the growing season, there were sufficient EC tower data to construct Keeling plots for 131 days (80%). Of these 131 days with complete data, there were significant slopes for 112 days (86%) using δ 2 H, and 89 days (68%) using δ 18 O. Keeling plots were not constructed for the dormant season, due to issues with δe modeling and unstable boundary layer conditions.

Craig-Gordon Model
Of the 164 days in the growing season, the CGM successfully ran for 148 (90%) of days. In this case, δ 18 O and δ 2 H performed equivalently in terms of producing values for δE, and the limiting factor for running the CGM was obtaining complete water vapor isotope data from the EC tower and CRDS. Sixteen days (10%) were missing or had incomplete water vapor isotope data. Relative humidity data were complete for the entire season, as well as soil moisture used to predict δe.

Daily fT for Entire Study Period
Days with insignificant Keeling plots (p-value > 0.05) were excluded for fT using the stable isotope method, along with days having an fT <0 or >1.0 for all methods. Day 131 was excluded, due to site setup and disturbance. For the first approach (T/ETEC), the average fT was 0.36 ± 0.31, and there was a total of 70 days (43%) with valid and complete data. For the second approach (T/ET0), the average fT was 0.28 ± 0.11, and there was a total of 140 days (86%) with valid and complete data. Both δ 2 H and δ 18 O had an average fT of 0.77 ± 0.20, though there were only valid data for 79 and 56 of the 164 days (48% and 34%), respectively. The results from raw calculated fT are reported in Figure 10.  On average, f T calculated from T/ET EC or T/ET 0 were 41% and 49% lower than those calculated from the stable isotope method. Minimum f T calculated from the stable isotope method was nearly always greater than the maximum value calculated by traditional methods (T/ET EC , T/ET 0 ). To remove this potential systematic bias from the estimates, f T from the stable isotope method was bias-corrected to ET 0 (Figure 11 left) and ET EC (Figure 11 right) using the average difference as the trends (direction and slope) were similar, despite the absolute magnitude being different. The average ET 0 bias-corrected stable isotope f T was 0.40 ± 0.15, and the ET EC bias-corrected f T was 0.32 ± 0.15, which are only 4-12% greater than traditional f T , as opposed to 41-49% before normalization.
Water 2020, 12, x FOR PEER REVIEW  15 of 20 On average, fT calculated from T/ETEC or T/ET0 were 41% and 49% lower than those calculated from the stable isotope method. Minimum fT calculated from the stable isotope method was nearly always greater than the maximum value calculated by traditional methods (T/ETEC, T/ET0). To remove this potential systematic bias from the estimates, fT from the stable isotope method was biascorrected to ET0 (Figure 11 left) and ETEC (Figure 11 right) using the average difference as the trends (direction and slope) were similar, despite the absolute magnitude being different. The average ET0 bias-corrected stable isotope fT was 0.40 ± 0.15, and the ETEC bias-corrected fT was 0.32 ± 0.15, which are only 4-12% greater than traditional fT, as opposed to 41-49% before normalization.

Discussion
This study demonstrated the utility of using a combination of stable isotopes, sap flux, and eddy covariance techniques to partition ET in an oak woodland and introduces a novel methodology using CRDS water vapor isotopes. The stable isotope technique has benefited from improved technological advancements with cavity ringdown spectroscopy and high-temporal resolution vapor collection systems that work in unison with eddy covariance and sap flux systems. However, the stable isotope method for partitioning ET did not compare favorably with more traditional techniques, such as using the difference of total ET from eddy covariance or ET0 along with T estimated from sap flow measurements, on any sampling day. There was a 41-49% overestimation of fT in this system when utilizing the stable isotope technique compared to ETEC or ET0 relative to sap-flux-based T. When using the average difference to normalize fT for either δ 18 O or δ 2 H the overestimation was reduced to 4-12%, which is within the range found from other studies [10]. This suggests that there may be a systematic bias to the CGM, which leads to the overestimation of fT in natural systems, and this bias should be investigated thoroughly in the future.
When comparing δ 18 O and δ 2 H within the stable isotope method, there was much agreement between the two, which contrasts with results and conclusions from other studies [10,21,32]. Typically, δ 2 H is utilized over δ 18 O, due to its higher degree of sensitivity. However, the high temporal-resolution data provided by the CRDS was able to provide the necessary data resolution for δ 18 O to have a similar performance to δ 2 H. Both δ 18 O and δ 2 H were able to satisfy requirements for the CGM 90% of the time, and δ 18 O only slightly underperformed in the Keeling Plot approach Figure 11. Stable isotope f T bias-corrected to (left) ET 0 ; (right) ET EC . Data for ET EC before DOY 241 were invalid and excluded from the analysis.

Discussion
This study demonstrated the utility of using a combination of stable isotopes, sap flux, and eddy covariance techniques to partition ET in an oak woodland and introduces a novel methodology using CRDS water vapor isotopes. The stable isotope technique has benefited from improved technological advancements with cavity ringdown spectroscopy and high-temporal resolution vapor collection systems that work in unison with eddy covariance and sap flux systems. However, the stable isotope method for partitioning ET did not compare favorably with more traditional techniques, such as using the difference of total ET from eddy covariance or ET 0 along with T estimated from sap flow measurements, on any sampling day. There was a 41-49% overestimation of f T in this system when utilizing the stable isotope technique compared to ET EC or ET 0 relative to sap-flux-based T. When using the average difference to normalize f T for either δ 18 O or δ 2 H the overestimation was reduced to 4-12%, which is within the range found from other studies [10]. This suggests that there may be a systematic bias to the CGM, which leads to the overestimation of f T in natural systems, and this bias should be investigated thoroughly in the future.
When comparing δ 18 O and δ 2 H within the stable isotope method, there was much agreement between the two, which contrasts with results and conclusions from other studies [10,21,32]. Typically, δ 2 H is utilized over δ 18 O, due to its higher degree of sensitivity. However, the high temporal-resolution data provided by the CRDS was able to provide the necessary data resolution for δ 18 O to have a similar performance to δ 2 H. Both δ 18 O and δ 2 H were able to satisfy requirements for the CGM 90% of the time, and δ 18 O only slightly underperformed in the Keeling Plot approach when compared to δ 2 H (68% and 86% of the study period). This may be in part to the more sensitive fractionation coefficient for δ 2 H relative to δ 18 O, which increases the potential for significant results, or there may have been evaporation in soil and twig sample vials during transit from the field to the lab (for IRMS analysis), which could have led to the error. However, these results complement those from the GMWL, indicating that the stable isotope approach can successfully be utilized within the framework provided in this study. The limiting factors in constructing Keeling plots at the daily scale using the CRDS and the EC tower are having a robust data collection system and a field crew that can maintain and repair the instruments in the case of malfunction or failure. While δ 2 H slightly outperformed δ 18 O, which was expected, δ 18 O performed more strongly than previous studies have suggested. This indicates that the high temporal-resolution data collected by the CRDS is truly necessary to make this approach viable at the daily scale for both stable isotopes of water.
A severe limiting factor to the utility of the stable isotope approach for partitioning ET, however, is the ability to obtain soil and xylem water samples from the field at fine enough resolutions to examine multi-day or seasonal trends in δ e and δ t . For this study, field samples were only feasibly obtained monthly, which limited the data available for partitioning initially to one day per month. However, due to partnership with NEON, a diverse suite of data was available at a much higher resolution, and soil and xylem water were successfully modeled using 30 min VPD and soil moisture data, respectively. This enabled us to partition ET from what was initially only six sampling events, to 79 and 56 days for δ 2 H and δ 18 O, respectively. While this was still only 56% and 48% of all days in the growing season, the stable isotope method was markedly improved by utilizing a combination of modeling and field sampling events. Our model also coincides with results from recent work that established tree water status, driven by soil water potential and atmospheric conditions, was the main reason for source water partitioning and isotopic fractionation during a 7-week lysimeter experiment in Switzerland [51].
While 0.16 m soil depth worked best during the growing season for this study, other depths may work better under different precipitation regimes. For instance, if the soil constantly remained unsaturated during the growing season, then 0.2 to 0.3 m may be more appropriate to model δ e dynamics. Similarly, under flooded conditions, it would be more appropriate to model δ e from soil moisture at the surface. While source water for trees was not independently identified during this study, that may be advantageous in the future to verify the assumption that δ X = δ T . Furthermore, modeled δ t in this study was based solely on xylem water from Quercus spp. twig samples. However, these values may change if other species in the tower footprint were sampled instead. A more accurate estimation of both T and δ t may be achieved by measuring sap flow in additional species and sampling those species for xylem water simultaneously. Steady-state assumptions were justified in this study, due to sampling in the mid-day and afternoon conditions, but there is an increasing body of evidence that in the early morning and evening, steady-state conditions are no longer present. This should be taken into consideration for future work; however, some studies have also shown that f T calculated both ways is similar at the end of the day [36].
When dormant season data was added to the δ e and δ t models, δ t remained significant as the water inside of the trees were not being utilized during dormancy; however, the δ e models became unreliable. This suggests that evaporation dynamics during the dormant season are different than during the growing season, and that these dynamics were not successfully captured in the CGM. Furthermore, atmospheric water vapor isotopes at the onset of and during the dormant season are considerably more variable than during the growing season. This is likely due to the physical barrier that leaves give at the top of the canopy, as well as the regulating properties of active transpiration on the boundary layer micrometeorology [52].
It should be noted that forests and woodlands inherently have more variability and heterogeneity when compared to engineered or agricultural systems like croplands, and this may result, in part, account for variance observed in the boundary layer. However, the combination of insignificant δ e models and sporadic water vapor concentrations just before and after trees lose their leaves suggests that the CGM cannot be utilized during the dormant season in such systems.
While the application of δ e and δ t modeling may be used to interpolate between sampling events during the growing season, the utility of the approach is severely limited during the dormant season, due to shifts in boundary layer conditions and evaporation dynamics at the soil surface. Additional experiments with the CGM using high temporal-resolution data collection systems similar to this study are warranted to reconcile discrepancies for CGM performance over changes in seasonality. While in this study, the assumption that T was negligible during the dormant season may be appropriate, this may not be the case in other systems.

Conclusions
This study demonstrated a novel methodology for partitioning ET using a combination of the high-resolution stable isotope, sap flux, and eddy covariance techniques over multiple seasons in a natural oak woodland site. To the best of our knowledge, this is the first report of utilizing a CRDS for ET partitioning at this timescale and under these conditions, and only the second report of direct partitioning comparisons between stable isotope derived f T and f T from sap flow and eddy covariance. Furthermore, we have developed predictive models for estimating the isotopic composition of the soil (δ e ) and xylem (δ t ) water (both components necessary for the stable isotope-based ET partitioning method) on days when field samples cannot be taken, by utilizing common environmental and micrometeorological data. With the high-frequency data made available through NEON, we found that δ 2 H and δ 18 O produced similar results, where, in the literature, δ 2 H is typically only utilized. The utility of the stable isotope method for partitioning ET has been markedly improved with the with the integration of high-temporal resolution data collection systems, such as an EC tower coupled with a CRDS, along with field sampling of soil and twigs. However, further refinement of the methodology, particularly in natural systems, is still needed to reconcile potential biases inherent in this approach, and to test the utility of the proposed δ e and δ t models under wide-ranging conditions. We recommend that future investigations continue to refine the use of high frequency stable isotope data for novel application in partitioning ET, given the advantages of a single integrated whole-canopy estimation.