One of the overarching objectives of the CHEESEHEAD19 study was to determine how spatial heterogeneity of the land surface impacts the Earth’s surface energy balance. The deployment of the EC tower network over a heterogeneous forested ecosystem provides a unique opportunity to test several hypotheses regarding the causes of energy balance non-closure.

First, we investigate the spatial and temporal variations in Imb across the CHEESHEAD19 domain over the course of a seasonal transition. We ask how similar is energy balance across the different sites? Are there areas which have consistently high/low Imb compared to other sites? And what factors are most strongly correlated with increased Imb? It is expected that the different surface characteristics across the domain will influence Imb. We test the hypothesis that flux towers with greater heterogeneity within their flux footprint will measure greater Imb than more homogeneous sites.

Additionally, recent CHEESEHEAD19 airborne flux results showed that mesoscale eddies were persistently contributing to H _S and H _L throughout the study domain, with diurnal and seasonal variations warranting further investigation (Paleri et al., 2022). Here we build on those results with the more continuous tower fluxes measurements to evaluate the hypothesis that mesoscale atmospheric structures arising from surface heterogeneity are an important cause of energy balance non-closure.

Lastly, the tower network also provided an opportunity to test the effectiveness of using spatial forms of EC to measure energy balance. These have been proposed as a potential way to reduce energy balance residuals by incorporating contributions from mesoscale motions better than conventional, single-site EC. We ask whether spatial forms of EC can reduce energy balance non-closure; evaluating the hypothesis that spatial EC will reduce non-closure by better incorporating spatial contributions to H _S and H _L not captured by traditional, single-site EC.

Energy balance components were continuously measured during the CHEESEHEAD19 field campaign at a network of field sites across the domain. This included 17 EC sites from the NSF Lower Atmosphere Observing Facility (LAOF) Integrated Surface Flux System (ISFS). Each site had an EC tower ranging from 32 m above ground level (AGL; 12 towers) to 25, 12, and 3 m AGL (1, 2, and 2 towers respectively). These were sited in a stratified random grid pattern within the domain, over a variety of surface types, including broadleaf deciduous forests, needleleaf evergreen forests, mixed forests, and wetlands (Figure 1). Each tower sampled three-dimensional wind velocity, temperature, and moisture at 20 Hz to determine surface-atmosphere fluxes (momentum flux, H _S , H _L ). Each tower was also instrumented with a net radiometer (for R _N ) and multilevel temperature and humidity measurements (for S _atmo ). Additionally, all sites had four levels of ground heat flux plates and thermometers to measure G and soil storage (S _soil ). An additional three sites with EC towers were excluded from this study because they did not have multi-level ground flux plates and thermometers to measure G and S _soil . Detailed information on instrument models/manufacturers can be found on the National Center for Atmospheric Research website at https://www.eol.ucar.edu/content/isfs-operations-cheesehead.

Temporal sensible and latent heat fluxes (Eqs 2, 3) from EC towers were calculated as H s = ρ a ¯ c p w ′ T ′ ¯ (2) H L = ρ a ¯ L v w ′ q ′ ¯ (3) where ρ a ¯ (mol m⁻³) is the mean dry air density, c _p (J kg⁻¹ K⁻¹) is specific heat capacity of dry air, w (m s⁻¹) is vertical wind speed, T (K) is dry air temperature, L _v (J kg⁻¹) is latent heat of vaporization, and q (kg kg⁻¹) is specific humidity. The dry air temperature was calculated from the sonic temperature after correction for the effect of water vapor on air density and speed of sound (Schotanus et al., 1983). Overbars correspond to time averages and primes indicate fluctuations about the mean. For temporal fluxes the fluctuation components for w, T, and q were calculated by subtracting the 30-min mean from the instantaneous 20 Hz values (e.g., w ′ = w − w ¯ ) on a site-by-site basis. Because some sites experienced directionally-dependent angles of attack, a single (not sectorial) planar fit coordinate rotation was applied to the wind vector for all analyses (except where otherwise stated).

While EC is generally considered a reliable method for measuring fluxes, there are certain limitations and uncertainties. To reduce measurement uncertainty, we applied several standard quality control techniques to the flux measurements. Specifically, we rejected flux intervals based on sensor diagnostic flags, integral turbulence characteristics tests (Foken and Wichura, 1996; Mauder and Foken, 2004), and nonstationarity (Blomquist et al., 2014). Uncertainties related directly to instrument biases (e.g., Frank et al., 2013) were not examined because individual sites lacked the instrument redundancies required to quantify such issues. While such uncertainties are likely to be partially responsible for observed energy balance non-closure, they are not expected to qualitatively change the energy balance forcings investigated herein.

In addition to traditional, temporal EC, we calculated fluxes using both spatial EC and spatio-temporal EC. These were used to test whether spatial forms of EC could better represent the energy balance over heterogeneous ecosystems. They were calculated as H S   s p a t i a l / s p a t i o − t e m p o r a l = ρ a c p w ″ T ″ (4) H L   s p a t i a l / s p a t i o − t e m p o r a l = ρ a L v w ″ q ″ (5) where angle brackets represent spatial averaging and double prime refers to fluctuations around the spatial mean. The key difference between spatial and spatio-temporal EC is in the way the velocity and scalar fluctuation components are calculated. For spatial EC the fluctuation terms are calculated by subtracting 20 Hz values by the 20 Hz spatial mean across all sites (e.g., w ″ = w − w ), resulting in 20 Hz fluxes. For spatio-temporal EC the fluctuation terms are calculated by subtracting 20 Hz values by the 30-min spatial mean across all sites (e.g., w ″ = w − w ¯ ). For this study the spatial EC fluxes were time-averaged from 20 Hz frequency to 30 min to achieve temporal equivalence in comparisons with temporal and spatio-temporal EC fluxes. Note that timestamps at all sites were continuously monitored and adjusted by a network time protocol software, using a GPS reference clock with a precise pulse-per-second signal. Times are expected to be correct across sites to within 50 microseconds.

One additional term calculated to identify spatial patterns in the surface energy fluxes was dispersive flux. Like spatial EC, it represents a covariance arising from spatial correlations of measurements over space, but is calculated with time-averaged quantities. The term can be added to single-site, temporal EC fluxes to produce a total, spatially averaged covariance (Raupach and Shaw, 1982). It is calculated using the equations: H S   d i s p = ρ a ¯ c p w ¯ ″ T ¯ ″ (6) H L   d i s p = ρ a ¯ L v w ¯ ″ q ¯ ″ (7)

In the simplest terms, a dispersive flux is the mean spatial flux over a period of time, resulting from spatially localized deviations from the mean flow. Dispersive fluxes represent a portion of the influence of larger-scale atmospheric motions on a flux area; scales that are difficult to measure by a single-site eddy covariance tower. However, this does not account for larger-scale divergence of the mean flows (Metzger, 2018; Morrison et al., 2021; Morrison et al., 2023).

In addition to flux measurements, the ISFS towers of the NSF LAOF were designed to measure the G, S _soil , and S _atmo terms. Soil heat fluxes (G) were measured using heat flux plates (HFT; REBS) installed at a depth (z _p ) of 5 cm. The details of this process are described in Oncley et al. (2007). Soil storage (Eq. 8) values in the layer above the heat flux plates were calculated using the equation: S s o i l = − c s o i l z p d T s o i l d t (8) where soil temperature (T _soil ) was calculated as the linear average of 4 temperature sensors (TP01; Hukseflux) evenly spaced through the first 5 cm of soil. The soil heat capacity (c _soil ) was calculated as a combination of water and dry soil heat capacities ( c s o i l = ρ w a t e r   q s o i l   c w a t e r + ρ d r y _ s o i l   c d r y _ s o i l ) weighted using soil moisture data (EC-5; Decagon) and bulk soil density measurements.

Atmospheric storage of heat from T and q (S _atmo = S _air-T + S _air-q ) below the EC flux sensor are represented by the Eqs 9, 10: S a i r − T = ∫ 0 z m ρ a c p d T d t d z (9) S a i r − q = ∫ 0 z m ρ a L v d q d t d z (10) where z _m is flux measurement height. For the calculation these integrals were discretized, using T and q measurements at z _m and from additional TRH sensors installed at 2 and 10 m AGL. Unlike soil measurements, T and q properties in the air are variable over short timescales. Therefore, to avoid introducing an unreasonable amount of noise, the T and q time derivatives were calculated as the difference between the means of the first and last 5-min periods of each 30-min flux interval divided by 1,500 s (i.e., 25 min).

At densely forested sites it was expected that some amount of energy would also be stored in plant biomass (S _bio ). Therefore, instruments aimed at estimating this additional storage term were deployed. S _bio was estimated from tree temperature measurements obtained at five ISFS sites from September 2019 to April 2020. At each site (NE2, NE3, NE4, SW2, SW4) two Onset Hobo MX-2304 were installed into two different trees at 2 m AGL. S _bio was calculated as S b i o = b i o m a s s ∙ c p ∙ d T d t (11) where c _p is the specific heat of wood and dT/dt is the change in temperature over time measured by the temperature sensors and biomass representing the mass of the forest. Biomass was estimated by using Wisconsin DNR forest lidar data to estimate mean tree height for each site (mean tree height = digital surface model–digital elevation model), then estimating diameter at breast height (DBH) from mean tree height using allometric equations and constants specific to the tree species at each site (Desai et al., 2007). Then mean DBH and mean tree height were used to calculate biomass per m² using allometric equations and constants.

The evaluation of spatial and spatio-temporal EC was incorporated into a larger sensitivity analysis, which evaluated the influence on energy balance closure of 1) EC type, 2) tower selection, and 3) coordinate rotation. Figure 2 shows all permutations of processing techniques for which H _S and H _L were calculated. The flux calculations influenced no other part of the energy balance equation (i.e., R _N , G, or S). Therefore, we interpreted the improvement to energy balance by observing how Imb varied for each permutation of techniques.

Preliminary analyses showed that the selection of suitable towers within the tower network can influence the calculated fluxes. Therefore, in the sensitivity analysis we calculated fluxes using 1) all towers and 2) towers ≥25 m AGL.

Coordinate rotation refers to the rotation of the wind vector used to calculate vertical velocity fluctuation used in the flux calculation. Two variations were applied to all EC types–single rotation and planar fit. A third coordinate rotation, a double rotation, was applied to temporal EC. Single rotation refers to an azimuthal rotation into the mean horizontal wind. Double rotation refers to an azimuthal rotation into the mean horizontal wind and a subsequent elevational rotation so that the mean vertical wind speed over 30 min equals zero. This rotation was excluded from both forms of spatial EC because setting w ¯ = 0 for all sites removes the spatial deviations in w expected to contribute to the spatially-calculated flux. Planar fit is also a double rotation, but one in which the wind vector is rotated onto the plane of the long-term (i.e., over weeks or months) mean wind field (mean crosswind component v ¯ = 0; but w ¯ is not necessarily 0). It is generally appropriate in temporal EC when a site has a consistent directionally-dependent angle of attack, because it quantifies vertical wind deviations from the long-term mean (Wilczak et al., 2001).

During CHEESEHEAD19 the magnitude of the energy imbalance (Imb) varied diurnally (Figure 3F) and seasonally (Figure 4C). Most sites exhibited a distinct diurnal trend of negative Imb at night and larger, positive Imb during the daytime, leading to a net daily positive Imb. Imb also varied across space, with some sites exhibiting characteristically high Imb throughout the experiment while other sites exhibited low Imb (Figure 3F).

In addition to changing diurnally, the energy balance components changed seasonally (Figure 4C). Positive Imb values decreased throughout the study period from their highest values in June (∼100 W m⁻²) to their lowest values in October (∼0 W m⁻²). These magnitudes correlate with the reduction in R _N throughout the season (Figure 4B). This results in the energy balance ratio staying constant near 0.8 through the seasonal transition (i.e., Imb changes in direct proportion to R _N ).

An investigation across the flux sites revealed relationships between canopy characteristics (mean canopy height and standard deviation) and energy balance components. Energy balance residual had positive relationships with both canopy height and canopy standard deviation (Figure 5). This appears at least partially because sites with shorter, scrubbier vegetation measured greater ground flux and soil storage, which accounted for a sizable portion of what otherwise would have been residual. This larger G and S _soil appeared to be the result of direct incoming solar radiation reaching the surface, causing greater heating than occurs at the surface under more dense canopies. These sites experienced G + S _soil sums that peaked at ∼150 W m⁻² during the day. This was an order of magnitude higher than peak G + S _soil sums at more densely forested sites.

The investigation into biomass storage revealed that S _bio also varied across the sites, with peak daytime storage ranging from 1 to 11 W m⁻² and a mean of 5 W m⁻² (Figure 6). This represents roughly 6.5% of the mean daytime peak in Imb of 76 W m⁻². This result supports previous findings which show that accounting for biomass heat storage improves energy balance closure, but that it cannot explain the entirety of the missing energy (e.g., Lindroth et al., 2010).

A major challenge facing energy balance studies is the different scales at which the measurements are taken. In an idealized scenario the footprints of the H _S , H _L , R _N , and G measurements would precisely match. In practice, there are orders of magnitude difference between them, with H _S and H _L footprints on the order of 100,000 m², R _N roughly 1,000 m², and G roughly 1 m². If the environmental conditions are homogeneous then the different size footprints should not influence the resulting energy balance characterizations. However, when a region is heterogeneous the different measurement footprints can cause mischaracterization of the energy balance for that region. This is particularly difficult when measuring over a forested site, as the installation of a tower generally requires some degree of open space (i.e., a gap in the trees). This means that the downward-facing components of the net radiometer are (at least partially) measuring unforested surfaces (e.g., grass, dirt, road) and will subsequently have a R _N not precisely representative of the flux footprint.

An additional issue related to these open areas is that, at some sites, the soil atop the heat flux plates that measured G received direct sunlight for part of the day. This is expected for grassland sites with limited canopy cover. Such sites have larger downward daytime G than forested sites (e.g., Figure 4B). However, for several forested sites the open patches near the tower resulted in large daytime G measurements uncharacteristic of those measured below an unbroken forest canopy. This is expected to artificially reduced Imb and is likely responsible for the large range across sites of midday Imb seen in Figure 3F.

Our investigation into spatial patterns in energy balance across the domain tested the hypothesis that flux towers with greater footprint heterogeneity would measure greater Imb than more homogeneous sites. The data supported this hypothesis. It was found that sites with higher canopies and more variable canopies (i.e., greater standard deviation) had larger energy balance residuals, both in magnitude and ratio. This corroborated several previous findings, including those of an EC tower network experiment in the Heihe River basin in China which found that in an oasis-desert environment surface heterogeneity within the flux footprints correlated with reduced energy balance closure (Xu et al., 2017; Zhou and Li, 2019). Additionally, across the continental range of sites included in FLUXNET, sites with more heterogeneous land surface (e.g., >σ ² _EVI) were found to have greater energy balance residuals (Stoy et al., 2013).

Part of the explanation for the positive relationship between vegetation height and variability appears to be due to measurement footprint discrepancies described above. Sites with shorter vegetation had less canopy cover. This led to larger G measurements, from a greater amount of incoming solar radiation reaching the surface. This G (interpreted as energy departing the exchange surface into the Earth) is subtracted from R _N to get the available energy side of Eq. 1, resulting in a reduction of the magnitude and ratio of the energy balance residual. However, for sites with larger canopies, the small footprint of the G measurement in the open canopy area near the towers was not representative of the entire flux footprints of the H _S and H _L measurements. Therefore, the combined atmospheric heat flux was likely lower for these sites than it would have been under homogeneous conditions matching the G footprint. This also helps explain why the two wetland sites did not match the trends of the forested sites (Figure 5). For the wetland sites the G measurement footprint was generally a better match to the R _N footprint and flux footprints.

The investigation of biomass storage showed that the omission of biomass storage plays a role in the existence of energy balance non-closure. Across five sites and 1 month of measurements, it was found that biomass storage accounted for 6.5% of the magnitude of daytime Imb. On days with large air temperature swings, daytime biomass storage values were frequently between 20 and 40 W m⁻². The relatively large biomass storage magnitudes on such days were due to large biomass dT/dt derivatives (Eq. 11) resulting from the positive relationship between biomass temperature and air temperature. While biomass storage is often omitted from energy balance studies for practical reasons, a thorough accounting of energy balance requires its inclusion. Here we show that while it is not the singular cause of non-closure, it appears to be an important component, especially under certain conditions.

The strength of mesoscale eddies was positively correlated with Imb (Figure 9). However, the factors that contribute to Imb appear to act consistently even under conditions with less contribution from mesoscale eddies. This was demonstrated by the consistent energy balance ratio across the seasonal transition (Figure 4C). One example of a proportional cause of non-closure in this dataset was time averaging in the flux calculation.

In traditional EC, time averaging is necessary to provide enough statistical information to calculate a representative flux. Ideally, the window is long enough to include contributions from all scales of turbulent eddies, but short enough to exclude larger scale, non-turbulent contributions (e.g., advection, diurnal cycles). This is the stationarity criterion. While no window length works perfectly for all situations, EC networks (e.g., Ameriflux, NEON, FLUXNET, etc.) generally set 30-min averaging windows as standard.

The failure of cospectral curves in Figure 10 to go to zero at lower frequencies suggests that the 30-min averaging window was responsible for eliminating some actual low frequency contribution to the flux. The magnitude of the missing low frequency contribution was greater on days with more developed mesoscale eddies (i.e., days with high w * ). However, the percentage of the missing contributions with respect to Imb was roughly equivalent to 10% for both high and low w * conditions. Because 30-min averaging windows are a common standard, this suggests that the high-pass filter caused by the averaging window may be a persistent cause of a small, but non-negligible, fraction of Imb across many field measurement campaigns.

Our primary objective of this study was to quantify the energy balance influence and improvements of specific data processing techniques. The unique experimental design of the CHEESEHEAD19 project enabled testing both spatial and spatio-temporal EC, as well as dispersive flux contributions to temporal EC fluxes. It should be noted that the experimental setup notably did not lend itself to addressing all possible causes of non-closure. For example, the lack of flux measurements at multiple heights and the absence of side-by-side towers precluded the calculation of vertical flux divergence and advective fluxes. The omission of these fluxes is expected to contribute to energy balance non-closure (Finnigan et al., 2003; Oncley et al., 2007; Morrison et al., 2021; 2023). Nevertheless, the testing of spatial forms of EC over a domain with the spatial scale of the CHEESEHEAD19 tower network represents a novel scientific endeavor.

Previous studies have suggested that accounting for dispersive flux may reduce energy balance non-closure (Kanda et al., 2004; Mauder et al., 2020; Mauder et al., 2021; Morrison et al., 2023). There are theoretical reasons for including dispersive fluxes. Namely, that they represent one contribution of mesoscale motions to the energy balance that cannot be measured by single EC towers (Mauder et al., 2020; Mauder et al., 2021). Fortunately, this term can be incorporated with straightforward calculations (i.e., Eqs. 6, 7) using data from a flux tower network. Because it relies solely on mean meteorological terms it can also be measured without flux towers, by having a network of basic meteorological stations (Mauder et al., 2008).

The dispersive fluxes calculated using the CHEESEHEAD flux network were relatively small compared to other energy balance components included in the analysis. The hourly bin medians of the diurnal cycle for both H _{S disp} and H _{L disp} (from ≥25 m towers) peaked at roughly 10 W m⁻² and −10 W m⁻² (Figure 11A), respectively, but averaged to 7 and −3.5 W m⁻² over the daily cycle. The existence of relatively small dispersive fluxes agrees with the findings of Morrison et al. (2023), who found that dispersive fluxes accounted for a negligible portion of the surface energy budget. The considerable degree of variability in the dispersive fluxes is potentially due to the variable nature of secondary circulations under different meteorological forcings (Paleri et al., 2022).

The temporal flux of the spatial means (Figure 11B) showed larger magnitude and more consistent diurnal cycle compared to dispersive fluxes. The flux direction in the diurnal cycle matched that of the temporal flux measurements from individual sites, albeit with values an order of magnitude lower. Engelmann and Bernhofer (2016) state that this term may collect information about larger and/or slower structures. However, the theoretical case for this term has not been made. Here we found that the term was equivalent to the difference between spatial and spatio-temporal fluxes. Therefore, it may play a useful role in the decomposition of components contributing to spatio-temporal flux measurements (i.e., those contributions that are purely spatial vs. those that arise from the time-averaging procedure in the spatio-temporal flux calculation).

In order to test the ability of spatial forms of EC to improve energy balance closure we conducted a sensitivity analysis in which energy balance residuals were calculated for different permutations of EC type (temporal, spatial, and spatio-temporal), tower selection (all towers, towers ≥25 m), and coordinate rotation (single rotation, double rotation, and planar fit). The inclusion of the latter two categories was to ensure that technical considerations were not responsible for Imb differences observed between EC types. An important consideration in the interpretation of the CHEESEHEAD19 sensitivity analysis is that reduction of energy balance closure does not necessarily correspond to improved (i.e., “truer”) flux measurements. This is because any change that increases H _S or H _L will have the effect of reducing Imb. And there are processing techniques that increase H _S and H _L that may not be theoretically defensible. But when the processing technique is both supported by theory and has the effect of reducing energy balance non-closure it provides a better justification for choosing that specific technique.

The coordinate rotation of the wind vector was found to be the category with the least impact on Imb (Figure 12A). It was observed that averaged over the course of the study period (and across all sites) the difference in Imb between single rotation and planar fit EC fluxes was negligible (<1 W m⁻²). However, a closer inspection of H _S and H _L revealed a range of influences on actual fluxes. Some sites, the ones with relatively minor angle of attack directional dependence, had nearly no difference between single rotation and planar fit fluxes, because the planar fit was effectively orthogonal to gravity just like single rotation. The sites with stronger directional dependence did witness up to 10% differences between the flux measurements calculated using different coordinate rotations. However, except for one site (NW4), the differences were equally distributed around a common mean (i.e., no coordinate rotation showed a consistent bias high or low). This resulted in the long-term analysis of Imb to reveal no impact from coordinate rotation. While further analyses on the specific deviations on Imb under different conditions will be performed, this analysis revealed that the type of coordinate rotation used is not a primary contributor to Imb.

With respect to EC type, we hypothesized that spatial forms of EC would reduce Imb by better incorporating contributions from mesoscale eddies. This hypothesis was informed by an LES analysis which demonstrated this result (Xu et al., 2020). We found that the spatial forms of EC tested in this field experiment did not support the stated hypothesis. In fact, we observed a substantial increase in Imb for spatial EC, averaging 50 W m⁻² greater during peak daytime conditions compared to temporal EC (Figure 12C). This corresponds to the results of Engelmann and Bernhofer (2016), who also observed increased Imb for spatial EC. In addition to spatial EC, we tested spatio-temporal EC fluxes for their ability to improve energy balance closure. We found that this method did improve closure (i.e., reduced Imb) compared to spatial EC, but still produced larger Imb values than temporal EC. The difference between spatio-temporal and temporal EC-derived Imb was reduced when flux calculations were performed using towers of similar height (i.e., the midday difference of 20 W m⁻² with all towers was reduced to 11 W m⁻² when selecting for tall towers; Figure 13). However, the persistently larger Imb for both spatial and spatio-temporal EC under all permutations shows that, under the test conditions, spatial forms of EC did not perform better at balancing the energy budget.

In their experiment, Engelmann and Bernhofer (2016) performed an additional EC technique which they call the combined EC method. For that calculation they added the spatial EC fluxes to the additional flux derived from the temporal covariance of the spatial mean values. This resulted in a large Imb reduction compared to spatial EC and a modest Imb reduction compared to their reference temporal EC fluxes. We also tested this approach and found that while it reduced Imb relative to spatial EC, it still produced Imb values slightly larger than those calculated using temporal EC. Interestingly, we found that the combined EC method of Engelmann and Bernhofer (2016) produced the same values as directly-calculated spatio-temporal fluxes (Eqs. 4, 5). This suggests that the combined EC method is a spatio-temporal flux calculation simply decomposed into its constituent parts (i.e., the spatial and temporal EC terms).

With this in mind, we can see that the spatio-temporal fluxes of the 9-tower Eddy Matrix of Engelmann and Bernhofer (2016) did slightly improve energy balance closure, while the 17-tower spatio-temporal fluxes from the CHEESEHEAD19 showed slightly poorer energy balance closure compared to temporal EC. One possible cause of this is the practical limitations of deploying a tower network over a large, heterogeneous domain, including greater distances between the sites, the variety of land surface types, and the varying flux measurement heights relative to canopy. Such physical characteristics may complicate the ways in which turbulence is measured using spatial and spatio-temporal EC. Indeed, the theoretical requirements for spatial EC suggest that in order for a spatial array of instruments to provide meaningful flux data, the measurements should be collected within a homogeneous turbulence field (Stull, 1988). This requirement would appear more difficult to fulfill in the broad domain of CHEESEHEAD19 compared to the small, idealized setup of Engelmann and Bernhofer (2016). However, such shortcomings (if relevant) would likely be resolved by the inclusion of all terms in the continuity equation, including horizontal advection and vertical flux divergence (Metzger, 2018; Morrison et al., 2021; 2023). Unfortunately, this could not be tested here due to the limitations of the experimental design. In spite of this, the spatial and spatio-temporal EC results from the CHEESEHEAD19 domain often show energy fluxes of similar magnitude compared to temporal EC fluxes. This shows that, at least to some degree, spatial forms of EC are capable of measuring over larger domains.