Evaporation and Transpiration From Multiple Proximal Forests and Wetlands

Climate change is intensifying the hydrologic cycle and altering ecosystem function, including water flux to the atmosphere through evapotranspiration (ET). ET is made up of evaporation (E) via non‐stomatal surfaces, and transpiration (T) through plant stomata which are impacted by global changes in different ways. E and T are difficult to measure independently at the ecosystem scale, especially across multiple sites that represent different land use and land management strategies. To address this gap in understanding, we applied flux variance similarity (FVS) to quantify how E and T differ across 13 different ecosystems measured using eddy covariance in a 10 × 10 km area from the CHEESEHEAD19 experiment in northern Wisconsin, USA. The study sites included eight forests with a large deciduous broadleaf component, three evergreen needleleaf forests, and two wetlands. Average T/ET for the study period averaged nearly 52% in forested sites and 45% in wetlands, with larger values after excluding periods following rain events when evaporation from canopy interception may be expected. A dominance analysis revealed that environmental variables explained on average 69% of the variance of half‐hourly T, which decreased from summer to autumn. Deciduous and evergreen forests showed similar E trajectories over time despite differences in vegetation phenology, and vapor pressure deficit explained some 13% of the variance E in wetlands but only 5% or less in forests. Retrieval of E and T within a dense network of flux towers lends confidence that FVS is a promising approach for comparing ecosystem hydrology across multiple sites to improve our process‐based understanding of ecosystem water fluxes.

and water transport to the atmospheric boundary layer to which the land surface is coupled (Oren et al., 1999).Understanding how different ecosystems regulate water supplies reveals insight into how much water is available for groundwater recharge and runoff (Anderson et al., 2017), and how ecosystem management decisions impact water cycling.
E and T can be measured using lysimeters, leaf-level gas exchange measurements, sap flow, and other techniques (Kool et al., 2014), but these measurements are difficult to scale up to the ecosystem level (Berkelhammer et al., 2016).Ecosystem ET is readily measured using eddy covariance (EC), but most studies do not seek to partition EC measured-ET into its components, which results in a paucity of ecosystem-scale E and T observations to understand the impacts of land use and land management changes on water cycling (Stoy et al., 2019).Water balance (Liu et al., 2016), machine learning (Granata et al., 2020), remote sensing (Martens et al., 2017), land surface models, and hydrologic models (Sun et al., 2017; X.-J.Zhang et al., 2014) can all be used to help estimate the contributions of E and T to ET at different spatial and temporal resolutions, but often rely on assumptions that require further validation.A large contributor to our uncertainty in understanding ET across multiple scales is the lack of ground-based E and T observations at the ecosystem scale (Rigden et al., 2018).
Most EC studies focus on a single ecosystem, paired ecosystems in close proximity, or multiple ecosystems across large spatial extents.Few studies have measured multiple ecosystems in close proximity, which is critical for understanding how different ecosystems contribute to the water balance of heterogeneous landscapes (Chu et al., 2021;Sun et al., 2017).We seek to address this by estimating E and T directly at the ecosystem scale at multiple sites within a diverse landscape using a dense array of EC towers within a 10 × 10 km area (Butterworth et al., 2021).We do so using an approach developed to partition carbon and water fluxes from high frequency EC measurements called flux variance similarity (FVS) (Scanlon & Kustas, 2012;Scanlon & Sahu, 2008).FVS assumes that stomatal and non-stomatal fluxes independently conform to Monin-Obukhov similarity theory and has been successfully applied to study E and T across multiple global ecosystems (Kustas et al., 2018;Perez-Priego et al., 2018;Rana et al., 2018;Scanlon & Kustas, 2012;Sulman et al., 2018;Wagle et al., 2020) but less frequently in forests and wetlands (Klosterhalfen et al., 2019;Sulman et al., 2018), leaving opportunities to better understand the pathways by which different ecosystems use water.
Here, we use FVS to directly partition E and T from ET measurements across multiple forest and wetland ecosystems in northern Wisconsin, USA.We briefly describe the study ecosystems and the conditions that result in successful partitioning of the FVS algorithm.We focus our analysis on quantifying the response of E and T from the study ecosystems to the diverse climatic conditions and phenological changes experienced during the measurement period across the transition from summer to autumn.We are interested to see how the algorithm detects changes in ecosystem water flux and expect to see patterns, namely that: (a) the wetlands and forested sites will partition E and T differently across the study period and that wetlands will support more E.More specifically, (b) T will dominate ET during summer due to high energy inputs and leaf area and E will dominate ET during autumn in the forested ecosystems due to decreased T, and that (c) T will differ little amongst forest ecosystems following the notion that it is a "conservative" quantity that is relatively insensitive to forest type (Oishi et al., 2010;Roberts, 1983).We focus our analysis on E and T, and their dynamics across ecosystems and time including large precipitation events and note that the extensive publicly available observations from the study ecosystems can be used for many lines of enquiry as outlined in Butterworth et al. (2021).

Study Sites
The Chequamegon Heterogeneous Ecosystem Energy-balance Study Enabled by a High-density Extensive Array of Detectors 2019 (CHEESEHEAD19) was designed to understand how the atmospheric boundary layer responds to heterogeneity in land surface structure and function (Butterworth et al., 2021).It deployed one of the highest density networks of EC measurements of surface-atmosphere energy fluxes to date: 19 flux towers in a 10 × 10 km area among an extensive array of land surface and atmospheric measurements in northern Wisconsin, USA, centered around the existing 447-m WLEF tower AmeriFlux site (US-PFa) (45.9459, −90.2723) (Figure 1, Table 1) near Park Falls, Wisconsin, USA in the Dfb (warm-summer humid continental) Köppen climate zone (Beck et al., 2018).The 3 July-13 October 2019 measurement period allowed us to capture how seasonal changes in vegetation structure and function impact E and T. Mean annual temperature is 4.33°C and mean annual precipitation is 823 mm with significant precipitation across all seasons.
Tower placement within the study domain followed a stratified random grid pattern with towers placed on average 1.4 km away from their nearest neighboring tower and 3.5 km from the WLEF tall tower (Davis et al., 2003;Xu et al., 2017).This partial randomization-taking distance to road, United States Forest Service (USFS)-owned land, and appropriate tree gap for tower into consideration-resulted in spatial variability in vegetation height and ecosystem type including observations in challenging flux measurement conditions.The study ecosystems included four evergreen needleleaf forests (ENF) sites dominated by Pinus spp., Picea mariana and/or Larix spp., two wetland sites (WET), six aspen (Populus tremuloides)-dominated sites, one maple (Acer saccharum)-dominated site, and four mixed forest sites, one of which was near a lake, and most of which are a mix of aspens and pines (Tables S1 and S9).All mixed and deciduous forests are denoted as deciduous broadleaf forests (DBF) for     the purposes of this study given that their change in leaf area across the seasonal transition will differ from ENF sites.Towers were mounted 32 m above ground level at the forest sites and between 1 and 3 m above ground level at the wetland sites (Table 1).Of the total 19 EC towers, 17 used an open path infrared gas analyzer rather than closed path systems to which an additional density term needs to be applied to the Fluxpart algorithm (Stoy et al., 2019), described below.We evaluated all 17 sites and excluded three sites that did not meet the set criteria of at least 20% physically realistic solutions from the FVS algorithm (Table 2).We further excluded one site whose measurement record began in mid-July rather than late June to not have to either gap fill data beyond its measured range to make the data record lengths compatible, or to shorten the time series of the remaining sites for compatibility and limit our data record as a consequence.The remaining 13 sites included 3 ENF, 2 WET, and 8 DBF ecosystems.We provide an extended description of all sites in Table S1 as well as vegetation classification fraction in the 90% flux footprint in Table S9.(Campbell Scientific, Inc.) to measure turbulent fluxes, a NR01 four-component radiometer (Hukseflux, Delft, The Netherlands), and SHT85 aspirated air temperature (Ta) and relative humidity (RH) sensors (NCAR) above the plant canopies from which VPD was calculated.Ta and RH were also measured at 2 m and at mid-canopy in the forests.Soil measurements included four-level soil temperature measurements (NCAR) at 0.6, 1.9, 3.1, and 4.4 cm depths, EC-5 SM probes (Decagon, Pullman, WA) at 5 cm, and HFT heat flux plates at 5 cm in forested sites and buried in mats in the wetlands.Precipitation was measured at a SURFRAD station (Augustine et al., 2000), a ground-based measurement system for continuous long-term measurements of climatic data, located in a grass field within the CHEESEHEAD19 study domain at 45.9458, −90.2944.From these observations we also calculated a one-dimensional water balance as the cumulative sum of precipitation minus ET for the study period.

Eddy Covariance and Flux Partitioning
The EC method is widely used to measure carbon dioxide and water vapor fluxes worldwide (Baldocchi, 2014).EC takes high-frequency (10-20 Hz) measurements of three-dimensional wind velocity, CO 2 , and H 2 O concentrations in the roughness sublayer over the plant canopy.Assuming that turbulent mixing is sufficient, half-hourly (or hourly) net ecosystem-scale flux measurements can be calculated.However, to improve our process-level understanding of ecosystems and to improve models, EC fluxes need to be partitioned into their separate parts: photosynthesis and respiration for carbon dioxide fluxes and E and T for water.This underlies the need for effective and accurate flux partitioning methods (Kool et al., 2014;Stoy et al., 2019).
A promising approach for partitioning ET into T and E is flux-variance similarity (FVS) (Scanlon & Kustas, 2010).FVS assumes that stomatal and non-stomatal turbulent fluxes conform independently to Monin-Obukhov Similarity Theory.For a brief conceptual description, assume that there are two end-member scenarios for a parcel of air transported from the surface: one that interacts only with stomata and one that does not.An eddy transported away from a surface that is respiring CO 2 and evaporating water through pathways other than stomata will have deviations from mean CO 2 concentration (c') and water vapor concentration (q') that are positively correlated.An eddy of air that interacts with a surface with open stomata will have a negative relationship between c' and q' that is determined by plant water use efficiency (WUE) with more water vapor and less carbon dioxide (due to photosynthesis) than surrounding air.WUE can be used to establish a relationship between the variance of CO 2 associated with stomatal intake and the correlation between stomatal and non-stomatal CO 2 exchange processes (Figure 2).
In this way, evapotranspiration (ET) can be partitioned into E and T by matching observed correlations of q' and c' that represent a combination of stomatal and non-stomatal processes to those that correspond to purely stomatal or non-stomatal processes (Scanlon & Sahu, 2008).Subsequent work by Skaggs et al. (2018) noted an algebraic solution to terms that had previously been solved numerically (Palatella et al., 2014).The analytical solution was incorporated into an open-source Python 3 module, Fluxpart (Skaggs et al., 2018), that implements the FVS method to partition E and T. The original FVS approach of Scanlon and Kustas (2010) used a simple WUE formulation, which assumes a constant ratio between leaf-internal CO 2 concentrations and atmospheric CO 2 concentration of 0.7 for C3 plants following Campbell and Norman (1998).Fluxpart also implements other methods for WUE, including models in which intercellular CO 2 varies as a function of VPD (Skaggs et al., 2018), and a model in which optimal stomatal behavior is assumed in response to VPD (Scanlon et al., 2019).We compared the constant ratio and optimality models and chose the constant ratio approach because it had a greater amount of successfully partitioned measurements, which we chose to analyze based on data acquisition (see Table S8 in Supporting Information S1 for all results).While Wagle et al. (2023) outline the benefits of using the optimality-based model for partitioning, they also found a substantially lower number of successful partitioned outputs.
Fluxpart delivers predictions of latent heat flux due to evaporation (LEe) and transpiration (LEt) from the EC latent heat flux measurements (LE) with units of W m −2 .E, T, and ET in mm half hr −1 were calculated by dividing their corresponding LE by the latent heat of vapourization using the equation of Henderson-Sellers (1984) and the density of water.It is important to note that the Fluxpart algorithm is not always able to make a calculation for different reasons.FVS and the EC methodology assume the surface-atmosphere exchange is well-represented by a turbulent flux, which is not always the case, especially at night (Zhu et al., 2006).Additionally, the FVS method is applicable only when a negative carbon flux exists due to photosynthesis, and positive carbon and water fluxes exist due to respiration, transpiration and evaporation, and when EC measured fluxes, variances, and correlations for CO 2 and water satisfy several constraints (see Scanlon et al. (2019), Equation 13).FVS results tend to match known fluxes generated using large eddy simulations when there is a clear separation of CO 2 and H 2 O sinks and sources (Klosterhalfen et al., 2019), and our results are subject to these limitations.As a consequence of these limitations, and due to the fact that EC time series contain multiple gaps (e.g., due to weather or instrumentation, etc.), we gap filled missing E, T, and ET data to create continuous time series over the study period.
Eddy covariance data and FVS outputs were gap filled using REddyProc (Wutzler et al., 2018) for R (R Core Team, 2023).REddyProc inputs half-hourly EC data and performs quality checks and data filtering based on measured flux and friction velocity to discard data collected under insufficient turbulence.LEe and LEt were gap filled using marginal distribution sampling, which replaces missing values based on observations measured within 1 hour of a total of seven adjacent days if climate conditions are similar under the assumption that these fluxes, like ET, vary little from day-to-day with similar climate forcing.These gap filled LEe and LEt time series were used to gap fill LE.
We analyzed the success rate of FVS partitioning to site and measurement characteristics using linear regression to explore if there were certain circumstances that led to greater partitioning success.Partitioning success was analyzed against leaf area index (LAI), canopy height, (z), and instrument height (h), as well as the distance between them (z − h) and their ratio (z/h).

Uncertainty Analysis
Eddy covariance measurements, like all measurements, contain uncertainty, which must be quantified for a robust estimate of net fluxes (Goulden et al., 1996).We used the approach of Richardson et al. (2006) which assumes that surface-atmosphere fluxes taken under similar climatic conditions at the same time on consecutive days should be approximately equal and uses differences between these measurements to estimate a random error, which trends to follow a Laplacian (double exponential) distribution (Hollinger & Richardson, 2005).We performed a sensitivity analysis on the conditions considered to be sufficiently similar by Richardson et al. (2006)-photosynthetic photon flux density differences of less than 37.5 μmol m −2 s −1 , air temperature differences less than 3°C, and wind speed (WS) differences less than 1 m s −1 -and found little justification to use other values.We added a wind direction threshold of a difference of no greater than 30° to account for heterogeneous vegetation surrounding some towers (Butterworth et al., 2021).We then used this "daily differencing" approach to calculate the standard deviation of measured or FVS-partitioned ET, T, and E by calculating the linear relationship between the magnitude of the flux and its standard deviation, and applying this linear model to all observations, which was taken to be the random error of the measurements.
Gap filling uncertainty was calculated using the standard deviation of the marginal distribution sampling calculated by REddyProc output (Wutzler et al., 2018) for ET, T, and E. We combined gap filling uncertainty for the gap filled quantities with random uncertainties of the measurements to create continuous vectors of the standard deviation of ET, T, and E at half-hourly time intervals.The resulting standard deviation estimates were autocorrelated, and uncertainty was propagated by estimating the effective number of observations after accounting for autocorrelation (nEff) using "lognorm::computeEffectiveAutoCorr" (Wutzler, 2021) in R. We used the observed and gap filled flux measurements to compute nEff because the magnitude of measurement uncertainty is a function of the magnitude of the flux (Richardson et al., 2006), and because this approach provided a more conservative estimate of the standard error of the mean flux value.We computed the standard error of the mean for each latent heat flux (LEe, LEt, and LE) for each ecosystem for which sufficient partitioned values were obtained (described further in Section 3), by first calculating the variance for random and gap filling uncertainty (x) separately using Where   2  is the mean value of the variance that describes the random or gap filling uncertainty of each flux (LEe, LEt, and LE).The total standard deviation was then calculated by summing variances.Significant differences in the mean value for each site were calculated using one-way ANOVA with Tukey's HSD post-hoc test.The Bonferroni filter was applied to adjust the 95% significance level for the 66 comparisons that resulted from analyzing differences amongst 13 sites: (1−0.95)/66= 0.00076.We also averaged observations by vegetation type to simplify the visual display.All site-level statistical analyses are applied to the raw and not the filtered data.

Dominance Analysis
To find a parsimonious explanation for the amount of variability in ET, T, and E determined by different environmental drivers, we chose a dominance analysis to account for the strong correlation between driving variables (Budescu, 1993).We used the dominance_analysis library in Python to calculate the relative importance of the drivers of the Penman-Monteith equation: available energy (Q, taken here to be Rn minus soil heat flux), VPD, WS (which impacts atmospheric and boundary layer conductance), SM, and temperature (which impacts the slope of the saturation vapor pressure curve) to ET, E, and T. Air temperature was chosen for ET and T calculations and 0.6 cm soil temperature was used in the dominance analysis for E.

Ecosystem Structure
We collected site-scale observations including LAI and z as noted (and the standard deviation of the latter), and soil data from the NRCS Soil Survey to categorize the dominant soil type for most of the flux footprint at each site in order to capture soil texture on an ecosystem scale.Of the 17 sites with open path gas analyzers, 5 were sandy soils, one was a loam, and the rest were sandy loams (Table 1).

Results
Rn (Figure 3a) and Ta (Figure 3b) decreased from summer to autumn with oscillations at approximately weekly time scales as a consequence of synoptic-scale meteorological drivers observed during the measurement period (Butterworth et al., 2021;Desai et al., 2021).A maximum temperature of 27.7°C was observed on 5 July with a minimum of 2.3°C on 5 October.There was an extended period of above-average air temperatures for over a week beginning 16 September; such positive temperature anomalies during autumn are a common feature of the climate of the upper midwestern U.S. The temperature difference between the top of the canopy and the bottom of the canopy decreased throughout the season (not shown) suggesting more energy reaching the ground as leaves begin to fall despite decreasing Rn throughout the seasonal shift.
The study year, 2019, had frequent rain events during the measurement period (Figure 4a): Precipitation across the state of Wisconsin deviated from average by a positive 39.1 mm in July, negative 6.1 mm in August, positive 83.3 mm in September, and positive 63.5 mm in October (NOAA 2022).The average daily precipitation event over the measurement period was 8.34 mm day −1 , excluding trace rainfall, with a maximum event of 52 mm on 3 September.The longest period without a rain event was 8 days in late-August.Soil moisture at most sites dropped to its lowest point in late-August at the end of this rain-free period (Figure 4b and Figure S1 in Supporting Information S1).The highest SM content followed the storm events in September, peaking on 13 September (Figure 4b).The one-dimensional water balance reached positive values during the measurement period in early September (Figure 4c) after the large precipitation events at most sites.Atmospheric VPD was similar among the forested ecosystems where it averaged (3.47 kPa) but was lower in the wetland ecosystems (1.32 kPa) (Figure 4d).

FVS Partitioning
Results relied on physically realistic solutions based on the constraints of the algorithm.The FVS algorithm satisfied the constraints of physically realistic solutions 40% of the time at SE6, an ENF site (Table 2), but never did at NE1.It satisfied these constraints less than 20% of the time at NW2 and SE5.These sites with <20% partitioning success were removed from subsequent analyses in the interest of studying sites with more available data as noted.NE3 was also left out of analysis, despite higher partitioning success, due to late start of data in mid-July that complicated comparisons.Subsequent analyses were performed for the remaining 13 sites as noted (Table 1).
The lack of partitioning success for most sites is due mainly to missing data, which accounted for between 15% and 42% of invalid partitioning data, depending on the site.The rest of the errors (Table 2) were largely due to an inability to satisfy realistic solutions and error messages indicating that the observation did not align with Monin-Obukhov similarity theory.Partitioning success increased with a greater difference between instrument height and canopy height (R 2 = 0.30), but no relationship between partitioning and position of measurements  1) during late-June to mid-October 2019.
relative to the roughness sublayer (using instrument height to canopy height ratio as a proxy) was found, nor between partitioning success and LAI, instrument height, or canopy height.

Evapotranspiration
ET generally decreased from summer to autumn (Figure 5a), in response to fluctuations in Rn at approximately weekly time scales that correspond to frontal weather systems (Figure 3a).The forested ecosystems maintained higher average ET than the wetland sites in mid-September (Figure 5a) after the rainy period in early September (Figure 4a).The ET (as well as E and T) uncertainty analysis indicated that cumulative fluxes from most sites were significantly different from most other sites (Table S2), but there was no significant difference amongst forests dominated by DBF or ENF trees (Figure 5) when comparing cumulative sums.Wetland sites had significantly lower average ET over the measurement period (244 mm), than forested sites (312 mm) (P < 0.0001, Figure 5b).

Evaporation
FVS-partitioned E across all sites was more aseasonal than perhaps expected given leaf fall in the DBF ecosystems during autumn (Figure 6a).The sum of E for the wetland sites during the measurement period (129 mm) was not significantly different from the mean of the forested sites (134 mm) (Figure 6b).Mean LEe was significantly different among most sites, using pairwise comparison, but LEe at SW1 and NE2 was not significantly different to most other sites due to their larger standard errors (Table S3).The site with highest E (178 mm), SW2, was dominated by aspen, while the site with the lowest E (100 mm) was SW4, a mixed hardwood forest.As with ET, forested ecosystems maintained higher average E than the wetland sites in mid-September (Figure 6a) after the rainy period in early September (Figure 4a).

Transpiration
T decreased steadily, on average, throughout the study period along with the decline in Rn (Figures 3a and 7a).Mean LEt was significantly different across most ecosystems (Table S4).SW2, an aspen stand DBF, had the highest cumulative T (209 mm), while NW1, a pine and spruce/fir forest, had the lowest amongst forest sites (163 mm).Wetland sites (121 mm) transpired significantly less than the forested sites (182 mm) across the measurement period (P < 0.0001).There was a period of increased transpiration in mid-September following the large precipitation events (Figure 4a) that was more prominent in the forested sites.

Uncertainty Analysis
The range of uncertainties in the FVS-partitioned water fluxes, calculated as the combination of random measurement uncertainty and gap filling uncertainty, ranged from 3.0% to 6.3% for LE (Table S5 in Supporting Information S1), 4.1% to 7.8% for LEt (Table S6 in Supporting Information S1), and 5.3% to 10.5% for LEe (Table S7 in Supporting Information S1).

The T/ET Ratio
T/ET (Figure 8) decreased, on average, over the course of the study with a notable decrease in early September during the wet period and subsequent recovery thereafter.Average T/ET for all sites over the course of the entire study was 52% with an average of 53% at both DBF and ENF sites and 45% at the wetland sites (Figure 9).

Dominance Analysis
The environmental variables Q, VPD, WS, SM, and temperature combined explained 73% of the variability in ET (Figure 10a) and 67% of the variability in T (Figure 10c), but only 49% of the variability in E (Figure 10b), on average, at the site scale.Q explained more of the variability in all water fluxes than other variables and explained 46% of the variability in ET on average, and 39% of the variability in E, but only 32% of the variability in T. Instead, VPD explained 25% of the variability of T at the site scale, but only 18% of ET and 6% of the variability in E, on average.Interestingly, VPD explained 13% of the variability in E in the wetland sites but on average only 5% in the forested ecosystems.Temperature explained on average 6% and 7% of the variability in ET and T, but only 2% of the variability in E. Soil moisture and WS explained 1% or less of the variability in all water fluxes during the measurement period.

Relationships Between Water Fluxes and Ecosystem Structure
Both canopy height and its standard deviation were significantly related to the cumulative sum of T but not E (p < 0.04), but not after excluding wetland sites with shorter canopy heights (Table 1).Sand and silt were unrelated to the cumulative sum of E and T, and loam was significantly and negatively related to T (p = 0.01), but again not after excluding wetland sites that had loamier soils.Rather, edaphic factors were significantly related to water fluxes after large rain events.For example, on the day of the large rain event (3 September 2019), E was significantly and negatively related to silt and loam content (p < 0.025) and positively related to sand content (p = 0.02).

Relationships Between Water Fluxes and Precipitation Events
To examine the impacts of canopy wetness on E/T partitioning, we recalculated T/ET ratios after excluding all days with more than 0.5 inches of precipitation and the day after as periods indicative of having a wet canopy following findings that temperate forests canopies largely dry after ∼24 hr (Klemm et al., 2002).The T/ET ratio increased by 10% on average across all sites during dry canopy conditions (Figure 11).

Overview
Partitioning ET to increase our understanding of its components is becoming increasingly important for water resource management in a changing climate.While T and E can be readily measured at point or individual tree scales, it remains difficult to measure at ecosystem scale across multiple sites.This study attempts to address this using a dense array of 17 EC towers in a 10 × 10 km area of which 13 ecosystems had a higher proportion of Fluxpart algorithm convergence and data availability.We expected that wetlands and the different types of forests will partition E and T differently during the seasonal transition from early July to early October in northern Wisconsin and expected T to dominate ET in the summer and decrease into autumn with little differentiation between ecosystem types.We also expected that E should dominate in DBF forested ecosystems in autumn due to loss of transpirable surfaces and T will differ little amongst forest ecosystems following relative insensitivity to forest type (Oishi et al., 2010;Roberts, 1983).
We found that ET decreased (Figure 5a) when energy inputs into our energy-limited ecosystems decreased (Figure 3a), as expected.This seasonal decline in ET was dominated by a corresponding decline in T (Figure 7a), whereas E was largely aseasonal across the measurement period (Figure 6a) with little difference among ecosystem types in its cumulative sum (Figure 6b) despite different responses to wet conditions that were related to soil texture.The cumulative sum of T across the measurement period differed little among forests as anticipated (Figure 7b), despite differences in forest composition (Butterworth et al., 2021;Murphy et al., 2022), and wetlands supported less cumulative T (Figure 7b).We describe the FVS partitioning outcomes that resulted in these findings before describing seasonal patterns of fluxes and their responses to micrometeorological variability across different ecosystems.Figure 11.The median of T/ET for all sites from 3 July to 13 October 2019 with ("ALL," darker shade) and without days that include at least 0.5 cm of rain and the subsequent day ("NO RAIN," lighter shade).Evergreen needleleaf forests sites are denoted in red, deciduous broadleaf forests sites in black, and wetlands in blue.

FVS Partitioning
Towers were placed in a quasi-random nature to address the CHEESEHEAD19 study objectives to understand the role of landscape heterogeneity in mesoscale atmospheric dynamics (Butterworth et al., 2021).As a result, some towers were located in less-than-ideal flux measurement terrain including heterogeneous ecosystems and at the aquatic/terrestrial interface, which contributed to the lack of Fluxpart algorithm success in some cases.For example, partitioning results were not obtained at NE1 because the Fluxpart models for estimating WUE require the instrument height to be above the canopy height (Table 1).The Fluxpart algorithm often failed at the wetland sites due to "negative VPD," which is consistent with the challenges of precisely calibrating humidity measurements near saturation (Meyer et al., 2008) noting that the near-surface air at the wetland ecosystems had lower VPD on average (Figure 4d).
Partitioning success is influenced by scalar-scalar correlations such as those between q and c, which in turn are influenced by sink-source distribution, height (atmospheric surface layer, roughness sublayer), surface heterogeneity, and canopy density (Klosterhalfen et al., 2019;Skaggs et al., 2018;Zahn et al., 2022).An addition to the inherent heterogeneity of the landscape, the senescence of (primarily) deciduous foliage produced additional patchiness to the landscape which can be an important factor in the validity of scalar similarity and may have interfered with scalar-scalar relations (Williams et al., 2007).We examined the relationship between partitioning success and site characteristics and found no significant correlation between LAI, the canopy height or instrument height.There was a weak positive correlation (R 2 = 0.30) between the distance between the canopy and instrument height and partitioning success.Klosterhalfen et al. (2019) found that partitioning success was correlated positively to the instrument to canopy height ratio, overall canopy height, instrument height, and LAI using large eddy simulation outputs, which differ from field observations of heterogeneous ecosystems.We found no correlation between partitioning success and the canopy to instrument height ratio, and no correlation to canopy or instrument height individually.If sampling heights are too far from the canopy, turbulence would likely fully mix the air parcel and distinguishing the scalars via partitioning may therefore no longer be accurate (Klosterhalfen et al., 2019;Zahn et al., 2022).Zahn et al. (2022) recommends an instrument to canopy height ratio (z/h) of less than or equal to 2, to fall outside of the roughness sublayer but not be too far above the canopy for flux partitioning algorithms like FVS.Of our 17 total sites, 7 fell outside of this range, one of which was below canopy height.FVS also had the lowest rate of convergence in a methodological comparison because of frequent failure to satisfy its assumptions, compared to conditional sampling methods for partitioning (Zahn et al., 2022); the main challenge arose from approximations used to represent the correlation coefficient between carbon and water vapor components.Despite the limitations of any flux partitioning approach, FVS provided novel insights into E and T in this study that help us understand the ecohydrology of multiple different ecosystems within a landscape.
We also note that results presented here used the assumption that Ci/Ca = 0.7 for C3 vegetation following the original FVS derivation (Scanlon & Kustas, 2012), because the approach that uses optimality theory to derive WUE developed by Scanlon et al. (2019) yielded lower partitioning success as described in Table S8 in Supporting Information S1.The constant Ci/Ca assumption often yields similar T:ET values to the optimality theory-based approach (Wagle et al., 2023) and additional comparisons and tests against independent data in well-instrumented ecosystems are likely to yield valuable insights into the best approach for water flux partitioning (Stoy et al., 2019).Ancillary measurements like sapflow and lysimeters were not available in the current project to compare against EC (Perez-Priego et al., 2017;Poyatos et al., 2016), due in no small part to its large scope and the complexity of operating 19 EC towers and other instrumentation necessary to address the overarching hypotheses of the CHEESE-HEAD19 project that the lack of EC balance closure is due to mesoscale atmospheric features that arise due to surface heterogeneity (Butterworth et al., 2021).However, many studies have been conducted in the vicinity of the study area using sap flux in the past (Ewers et al., 2002;Mackay et al., 2002Mackay et al., , 2007;;Tang et al., 2006).Mackay et al. (2007), for example, found that T in wetland sites was more sensitive to VPD and radiation in summer months, and we found that E in wetland sites was more related to VPD than forested ecosystems, which helps us understand controls over its variability, which are often rather muted and less related to microclimatic variability (Figure 10).

Controls Over Evaporation
The daily sum of E varied little throughout the shift from late June until mid-October (Figure 6b).This finding is supported by Paul-Limoges et al. (2020) who also found quasi-constant daily E below the canopy throughout seasonal changes in a deciduous forest in Switzerland.E is strongly associated with changes in SM under water limited conditions (Or & Lehmann, 2019;Perez-Priego et al., 2018).When the moisture content of the soil surface is close to saturation, atmospheric conditions control E. Our observations typically fell between saturation and water limitation during the measurement period (e.g., Figure 4b) such that we found weak relationships between both edaphic variables and atmospheric drivers when considering seasonal sums.We did find relationships immediately after large rain events-namely that E was positively related to sand content but negatively related to silt and loam content, and that the T/ET ratio increased by some 10% when excluding days with rain and the subsequent day (Figure 11).
An explanation follows from recent advances in soil E modeling.Based on the widely used equation proposed by Gardner (1959) which suggests that evaporation decreases as a function of the square root of time following precipitation events, Brutsaert (2014) proposed an exponential decay function of SM drawdown.Or and Lehmann (2019) introduced a model based on the notion that soil is an "evaporative capacitor" with recharge and discharge based on precipitation statistics and soil physical characteristics.In this model, stage-I of drying soil is governed by atmospheric conditions, as capillary flow through the soil is sufficient to satisfy evaporative demand and E is only limited by available energy in the gradient between upper layers of the soil and the atmosphere.During stage-II, E becomes primarily a function of soil water content and soil hydraulic properties: evaporation rates drop significantly as the soil continues to dry and the system becomes water-limited.The transition between the first and second stage is dictated by storage, the product of average water content and the characteristic length (L c ) of the soil-the limiting depth for capillary extraction of water from deeper soils-which varies by soil texture, with shorter lengths for both coarse and fine-textured soils (Or & Lehmann, 2019).The characteristic length is considered a good indicator to encompass soil hydraulic properties' effect on evaporation (Schneider et al., 2021).
We can examine our observations in the context of these models.Frequent precipitation inputs in our study domain are consistent with a situation where the time-varying component of E during dry-downs played a minor role in the observed time-series, on average.The range in soil textures across the study domain indicated that we would expect to see evaporative differences at some sites during different climatic conditions.We found that while mean E at most sites was statistically different from each other, there was very little overall difference in cumulative LEe between them despite the variability in soil type.This is consistent with the notion that the ecosystems were rarely if ever water-limited in a way that would influence evaporation rates due to L c during the study period.Soil moisture was variable throughout the summer but tended to increase through autumn (Table S1), keeping characteristic lengths short and E largely under atmospheric control.This is supported by Schneider et al. (2021) who found no difference in E between different soil types, while other studies have found distinct differences between evaporation and soil texture at water-limited sites (Merlin et al., 2016).We find that differences in E are situational following large precipitation events in the study ecosystems that did not sum to significant seasonal differences in our case (Figure 6).
During autumn, as total Rn decreases (Figure 3a), more sunlight can penetrate the canopy as leaves begin to fall during senescence, especially in the DBF forests.We posit that declining leaf area compensated for declining Rn through the seasonal transition such that changes in subcanopy Rn were muted, and as a result the subcanopy Rn declined at a lower rate than above-canopy Rn.Coupled with sporadic precipitation events which maintained SM levels in the later months, E/ET increased as T declined (Figure 7a).This resulted in a situation where E is quasi-constant throughout the season while T declined steadily due to its close coupling with Rn.We also found that VPD explained a larger proportion of the variance in E in wetlands than forests (Figure 10) which is interesting because VPD is characteristically lower in wetlands near at the study sites (Figure 4d).These findings are consistent with the notion that evaporation from standing water in the wetlands is more related to atmospheric demand than in forests where standing water is less frequently observed.

Controls Over Transpiration
T is driven by differences in atmospheric water and soil water potentials but, unlike E, is regulated through plant stomata.Q was the strongest determinant of T, which explained 34% of the average variance across all sites in our energy-limited ecosystems (Figure 10).High VPD typically causes plants to close their stomata to minimize water loss (Monteith, 1995), but high VPD conditions tend to co-occur with high radiation therefore making it difficult to separate their effects (Grossiord et al., 2020), hence our use of dominance analysis to separate the effects of covarying drivers.T tends increase as VPD increases to a certain threshold depending on the ecosystem (Ficklin & Novick, 2017;Franks et al., 1997;Marchin et al., 2016;Sulman et al., 2016;Will et al., 2013).T positively covaried with VPD in the study ecosystems, but it infrequently reached the limiting threshold commonly taken to be ∼1 kPa (Körner, 1995;Oren et al., 1999).VPD remained below ∼1 kPa 90% of the time or more across the study ecosystems during the measurement period.

T/ET
Average T/ET from the July to October study period was 52% across all sites, with an average of 53% among the forested sites and 45% among the wetland sites.These numbers fall within the wide range established by many previous global studies of various ecosystems of 40%-90% (Good et al., 2015;Schlesinger & Jasechko, 2014;L. L. Wang et al., 2014;Wei et al., 2017;Zhou et al., 2016).This study is on the lower end of the previously reported range of 40%-86% for temperate forests (Schlesinger & Jasechko, 2014) and is lower than estimates generated using isotopic approaches on the order of 64% (Good et al., 2015).However, Paschalis et al. (2018) found a significant decline in T/ET when moving from dry to wetter areas-about a 10% decrease-which may explain why we see lower T/ET than some other studies in our energy-limited ecosystems, as precipitation during 2019 exceeded annual averages compared to previous years.L. Wang et al. (2014) showed an exponential relationship between LAI and T/ET indicating that vegetative control over T/ET occurs over the lower LAI range.However, when only natural vegetation sites (with LAI > 1 m 2 m −2 ) were considered, there was a negligible dependence of T/ET on LAI (Paschalis et al., 2018).LAI data collected for five mixed forest sites on 25 June when LAI was measured (NE2, NE4, SE3, SE6, SW4) revealed no relationship between T/ET and LAI; however, all sites had LAI > 1 m 2 m −2 .This finding is supported by Fatichi and Pappas (2017) and Berkelhammer et al. (2016) who found that LAI matched seasonal dynamics of T/ET but not diel, daily or annual timescales.In contrast, Nelson et al. (2020) found that T/ET varied much more between sites than different years at the same site, indicating more reliance on site characteristics than climatic variables.This finding is supported by other recent studies, who have also found a strong relationship between LAI and T (Paul-Limoges et al., 2022).
Wetland T/ET is significantly altered by structural factors such as open water, plant species and diversity, as well as environmental factors like diurnal fluctuations in air and water temperature and water table depths (Drexler et al., 2004;Eichelmann et al., 2018).Between wetland sites, the most important factor affecting ET levels is the proportion of open water versus vegetation cover (Eichelmann et al., 2018).Shorter vegetation, such as tussock grasses which can often be found in wetlands, optimize leaf structure to minimize water loss in light rich environments (Givnish, 1988) resulting in lower transpiration but greater WUE compared to forested sites.As the area of open water increases, E increases more than T, as we found evidence of decreased T/ET in wetlands compared to the forested sites.

Uncertainty Analysis
Uncertainties in the FVS-partitioned water fluxes ranged between 3.0% and 10.5%, which was similar to or slightly larger than random measurement uncertainties (3%-6%) from other forested ecosystems (Goulden et al., 1996;Oren et al., 2006;Stoy et al., 2006).Earlier approaches assumed that LE can be gap filled with a high degree of success due to their relatively predictable response to environmental conditions (Falge et al., 2001).For example, REddyProc tends to estimate gaps in latent heat flux measurements with a low degree of uncertainty (Foltýnová et al., 2020).Nevertheless, future efforts should examine the ability of gap filling methods to accurately simulate LEt and LEe.More LE gap filling methods are being developed (e.g., Khan et al., 2021) but it needs to be determined if they will offer improved predictive skill.We elected not to estimate the "spatial" uncertainty of estimating heterogeneous ecosystems using a single point measurement (Oren et al., 2006) as we did not have multiple towers in individual ecosystems to do so and sought to avoid unconstrained estimation.We also had little basis for estimating bias errors including potential underestimation of LE due to the known challenges of lack of energy balance closure of EC measurements (Foken, 2008;Leuning et al., 2012;Stoy et al., 2013) that the CHEESEHEAD19 experiment as a whole is designed in part to address (Butterworth et al., 2021).However, independent approaches are converging on the notion that underestimated turbulent flux terms may be due more to sensible than latent heat fluxes (Charuchittipan et al., 2014;Gerken et al., 2018;Mauder et al., 2020;Paleri et al., 2022).Ongoing efforts to parameterize underestimated turbulent fluxes require estimates of atmospheric boundary layer heights (Mauder et al., 2021), which are difficult to quantify using individual flux towers, but the extensive land surface and atmospheric observations available from CHEESEHEAD19 (Butterworth et al., 2021) provide opportunities to better-understand potential bias uncertainties in EC measurements due to mesoscale atmospheric motions.
It is important to note that the unknown uncertainties discussed above and different methods to estimate WUE in the FVS algorithm also impact results.FVS may overestimate soil E (Klosterhalfen et al., 2019) while conditional sampling methods may underestimate it in comparison to chamber measurements (Thomas et al., 2008;Zahn et al., 2022), with implications for our understanding of E and T. Ongoing efforts to improve water flux partitioning including optimality solutions for WUE in flux variance similarity (Scanlon et al., 2019) and new conditional sampling methodologies (Zahn et al., 2022) will ideally improve the accuracy of these methods to create defensible estimates of water flux terms (Wagle et al., 2023).
Interception remains a key uncertainty in flux partitioning as it can account for up to 15%-30% of incident precipitation to the atmosphere (Crockford & Richardson, 2000), but attempts to capture its magnitude through modeling remains a challenge (De Kauwe et al., 2013).This remains a critical area of research that needs further exploration (Stoy et al., 2019).The result that T/ET increases after excluding periods when canopies are likely wet after rain events suggest that interception is partitioned by FVS largely into E despite similar source areas of interception and T (Figure 11).If the increase in E/ET after rain events is largely attributable to interception, the magnitude of intercepted evaporation at our sites during the study period is of similar magnitude to other studies (e.g., Crockford & Richardson, 2000).

Implications for Forests, Wetlands, and Water Management
Our results have implications for forest hydrology, but few studies have looked at the influence of forest management practices on evapotranspiration (Komatsu & Kume, 2020).Practices like forest thinning are commonly used as a water saving strategy but are understudied, particularly in wet environments (Sun et al., 2017).Forest thinning, which influences ET partitioning with increased light penetration due to a more open canopy, alters the microclimate of the ground layer.Sun et al. (2017) found that ET decreased by 20% after thinning, yielding an increase in water yield to the watershed, as supported by other studies (Dung et al., 2012;Hawthorne et al., 2013).Most forest hydrology research involves changes in runoff rather than change in ET; this omission-due in part to the relative difficulties in measuring ET versus runoff-can hinder accurate modeling and therefore optimal management (Komatsu & Kume, 2020).
Wetlands are essential climate regulators, but are globally declining three times faster than forests, with most losses due to land use change, agriculture and climate change (Finlayson & Davidson, 2018;Granata et al., 2020).Considering both wetland losses and wetland restoration efforts underway-such as converting crop lands back to wetlands-understanding their hydrologic impacts is imperative.Differences in evapotranspiration from wetland sites (especially with open water) are significant compared to drained agricultural sites, with much higher ET in the wetlands (Eichelmann et al., 2018).Not only is ET higher, but the use of water is much different; T/ET has found to be between 31% and 37% in wetlands (slightly smaller than our observations here) but on the order of 70% in the cropping systems that often replace them (Wang-Erlandsson et al., 2014;Wei et al., 2017).Additionally, increasing global temperatures could have significant implications for evaporative loss from wetlands as it has been shown that air and water temperatures are strong drivers of nighttime ET, which is dominated by E in these systems (Eichelmann et al., 2018(Eichelmann et al., , 2022) ) as we also find here.When wetlands are drained there are detrimental effects for water storage and groundwater inflation (van der Kamp & Hayashi, 2009).Improving our understanding of the contributions of E and T to ET is essential to understanding land use and climate change impacts on water cycling in these systems.

Conclusion
We applied FVS to partition ET using high frequency data from EC towers to investigate the role of vegetation and seasonal dynamics in E and T. On average T accounted for some 53% of ET at forested sites and 45% for wetlands, emphasizing a lower contribution of E in forests.E was relatively aseasonal and independent of ecosystem type throughout the study period due to the frequent precipitation but differed after large precipitation events as a function of soil type.T is highly correlated with climatic variables, especially compared to E, and varied significantly between wetlands and forested ecosystems.Wisconsin has seen a 15% increase in annual precipitation since 1950, with most extreme increases dominated by seasonal transitions (i.e., spring and fall) and this trend is expected to continue.Partitioning ET into its components during a seasonal transition across multiple ecosystems provides new insights for understanding how different ecosystems use water, with implications for hydrologic modeling in an era of rapid land use and climate change.CHEESEHEAD19 was primarily supported from the National Science Foundation Award 1822420 in addition to support from the DOE Ameriflux Network Management Project and NOAA.Brian Butterworth was additionally supported by the NOAA Physical Sciences Laboratory.NOAASURFRAD measurements were provided by NOAA's Climate Program Office, Climate Observations and Monitoring Program Project NA20OAR4310338.We thank NCAR's Earth Observing Laboratory, sponsored by the National Science Foundation for their operational, technical, and scientific support.We thank Anita Thompson for insightful comments on the manuscript.This project would not be possible without the support provided by the U.S. Forest Service and Wisconsin Educational Communications Board.We wish to acknowledge the Indigenous Nations on whose ancestral lands and ceded territories the study took place and recognize that the infrastructure used for this project is built on Indigenous land.We recognize the Ojibwe people as past, present, and future caretakers of this land whose stewardship of the region was interrupted through their physical removal by the 1830 Indian Removal Act and through US Assimilation policies explicitly designed to eradicate Indigenous language and ways of being until the 1970s.

Figure 1 .
Figure 1.A map of 17 of the 19 CHEESEHEAD19 eddy covariance study sites with open path infrared gas analyzers in a 10 × 10 km area in northern Wisconsin, USA.Sites are labeled with abbreviations used in the CHEESEHEAD19 project (Butterworth et al., 2021) (Table 1).Black markers indicate deciduous broadleaf forests, red indicates evergreen needleleaf forests, blue indicates wetland sites, and the red star indicates the WLEF tower (GIS & Wisconsin State Cartographer's Office, 2022, Esri, 2023).

Figure 2 .
Figure 2. A schematic representation of the basics of flux variance and correlation between stomatal and non-stomatal CO 2 and water vapor exchange processes.

Figure 3 .
Figure 3. Thermodynamic variables including the (a) half hourly net radiation (Rn) and (b) air temperature (Ta) in black and soil temperature (Ts) in red for a representative tower, NE2 (Figure 1, Table1) during late-June to mid-October 2019.

Figure 4 .
Figure 4. Hydrological variables during the 3 July to 13 October 2019, measurement period.(a) Precipitation measured at the SURFRAD station, (b) soil water content for a representative site, NE2 (Table1) (See FigureS1in Supporting Information S1 for more sites), (c) the one-dimensional water balance (the cumulative sum of evapotranspiration subtracted from cumulative sum of precipitation) at all sites averaged by ecosystem type, and (d) daily averaged vapor pressure deficit by ecosystem type.Evergreen needleleaf forests sites are denoted in red, deciduous broadleaf forests sites in black, and wetlands in blue.

Figure 5 .
Figure 5. (a) The daily median evapotranspiration (ET) from 3 July to 13 October 2019, by ecosystem type and (b) the cumulative sum of evapotranspiration by ecosystem type.Evergreen needleleaf forests sites are denoted in red, deciduous broadleaf forests sites in black, and wetlands in blue.

Figure 6 .
Figure 6.(a) The daily median evaporation (E) from 3 July to 11 October 2019 and (b) cumulative sum across all sites.Evergreen needleleaf forests sites are denoted in red, deciduous sites in black, and wetlands in blue.

Figure 7 .
Figure 7. (a) The daily median transpiration and (b) cumulative sum of transpiration from 3 July to 13 October 2019.Evergreen needleleaf forests sites are denoted in red, deciduous broadleaf forests sites in black, and wetlands in blue.

Figure 8 .
Figure 8.The daily mean of T/ET for all sites from 3 July to 13 October 2019.Evergreen needleleaf forests sites are denoted in red, deciduous broadleaf forests sites in black, and wetlands in blue.

Figure 9 .
Figure9.The seasonal median of T/ET (color block) and E/ET (empty block) and for all sites from 3 July to 11 October 2019.Evergreen needleleaf forests sites are denoted in red, deciduous broadleaf forests sites in black, and wetlands in blue.

Figure 10 .
Figure 10.Kernel density estimates of the amount of variability in (a) ET, (b) E, and (c) T explained by available energy (Q), vapor pressure deficit, air (for ET and T) or soil (for E) temperature (Temp.),wind speed, soil moisture, and the total variance explained by all selected variables.

Table
). Black markers indicate deciduous broadleaf forests, red indicates evergreen needleleaf forests, blue indicates wetland sites, and the red star indicates the WLEF tower (GIS & Wisconsin State Cartographer's Office, 2022, Esri, 2023).
Note.Italicized sites did not pass data availability thresholds as described in the text.

Table 1
All CHEESEHEAD19 Eddy Covariance Sites With Open-Path Infrared Gas Analyzers With Site Descriptions and Characteristics 10.1029/2022WR033757 4 of 23

Table 2
The Likelihood of Fluxpart Success by Site and Error Codes for the Cause of Failure