Mechanistic links between underestimated CO2 fluxes and non-closure of the surface energy balance in a semi-arid sagebrush ecosystem

The surface energy balance non-closure problem in eddy covariance (EC) studies has been largely attributed to the influence of large-scale turbulent eddies (hereafter large eddies) on latent and sensible heat fluxes. However, how such large eddies concurrently affect CO2 fluxes remains less studied and mechanistic links between the energy balance non-closure and CO2 fluxes are not well understood. Here, using turbulence data collected from an EC tower over a sagebrush ecosystem during two growing seasons, we decomposed the turbulence data into small and large eddies at a cutoff frequency and analyzed their contributions to the fluxes. We found that the magnitude of CO2 fluxes decreased concurrently with decreased sensible and latent heat fluxes (and thus increased energy balance non-closure), primarily caused by large turbulent eddies. The contributions of such large eddies to fluxes are dependent not only upon their magnitudes of vertical velocity (w) and scalars (i.e. temperature, water vapor density, and CO2 concentration), but also upon the phase differences between the large eddies of w and scalars via their covariances. Enlarged phase differences between large eddies of w and these scalars simultaneously led to reductions in the magnitudes of both CO2 and heat fluxes, linking the lower CO2 fluxes to energy balance non-closure. Such increased phase differences of large eddies were caused by changes in the structures of large eddies from unstable to near neutral conditions. Given widespread observations in non-closure in the flux community, the processes identified here may bias CO2 fluxes at many sites and cause upscaled regional and global budgets to be underestimated. More studies are needed to investigate how landscape heterogeneity influences CO2 fluxes through the influence of associated large eddies on flux exchange.


Introduction
The eddy covariance (EC) method has been widely used for assessing the turbulent exchanges of heat, water vapor, and CO 2 between land surfaces and the atmosphere for decades (Baldocchi et al 2001, Wilson et al 2002, Liu et al 2005, Wang et al 2010, Novick et al 2018. Most studies using the EC method indicate that the sum of the measured sensible (H) and latent (LE) heat fluxes fails to balance the available energy (R n − G 0 , where R n is the net radiation and G 0 is the ground heat flux) (Wilson et al 2002, Foken 2008, Leuning et al 2012, Stoy et al 2013, Gao et al 2017a, Reed et al 2018. The average ratio of H+LE to Rn -G is about 85%, known as the nonclosure problem of the surface energy balance. The lack of surface energy balance closure raises the concern that the measured CO 2 fluxes (Fc) may also be in error (Wilson et al 2002).
In case of similar sources and sinks for heat, water vapor, and CO 2 , atmospheric turbulent exchange processes that lead to underestimated H and LE may also cause a reduction in CO 2 fluxes. Indeed, previous studies have indicated a possible link between energy balance non-closure and CO 2 fluxes (e.g. Twine et al 2000, Wilson et al 2002. Using data measured at multiple FLUXNET sites in diverse ecosystems and distinct climate conditions during summer seasons, Wilson et al (2002) revealed that the magnitude of Fc increased as the energy balance closure improved, after taking into account the environmental variables. However, mechanistic links between underestimated H and LE (i.e. non-closure of surface energy balance) and Fc have not been well explored.
Numerous mechanisms have been proposed to explain the non-closure of the energy balance, including violations in the assumptions of steady state and horizontally homogeneous landscapes for the EC measurements, errors in measurements of each component of the energy balance, mismatch in instrument footprints, advection, and inadequate sampling of large-scale turbulent eddies (hereafter large eddies) (Wilson et al 2002, Mauder et al 2007, Foken 2008, Leuning et al 2012, Wohlfahrt and Widmoser 2013, Eder et al 2015, Gao et al 2016, 2017a, Zhou et al 2018. Of these factors, large eddies are one of the primary contributors to underestimated H and LE (Kanda et al 2004, Foken et al 2011, Gao et al 2016, 2017a, Zhou et al 2018. These large eddies can be induced by a variety of sources, such as horizontal advection, detached and attached eddies, convective motions, and secondary circulations, casting low-frequency variations in turbulence time-series data (e.g. vertical velocity w, air temperature T, and water vapor density q) (Zhang et al 2010). Low-frequency motions (i.e. large eddies) contribute to variances of the corresponding quantities and covariances between w and scalars (i.e. fluxes). However, the covariances between large eddies of w and those of scalars are dependent not only upon the magnitudes of large eddies of w and scalars but also upon the phase differences between them. The phase difference can be understood as a measure of the synchronization of two signals with smaller phase differences suggesting better synchronization (Li and Bou-Zeid 2011). Using a dataset measured over a flood-irrigated cotton field where LE is roughly commensurate with R , n Gao et al (2017a) found that enlarged phase differences between large eddies of w and q lead to decreased LE and thus decreased surface energy balance closure. However, it has not been clear whether large eddies, which lead to the underestimated H and LE and thus the non-closure, also cause a reduction in CO 2 fluxes, and whether these previously identified mechanisms regarding the influence of phase differences on heat fluxes are applicable to CO 2 fluxes.
In this study, we aim to provide a mechanistic explanation of the links between underestimated H and LE (and thus surface energy balance non-closure) and underestimated Fc. Here the underestimated CO 2 fluxes refer to the underestimated magnitude of net ecosystem exchange of CO 2 (NEE), since negative values of Fc indicate carbon uptake. We examine how large eddies affect both heat fluxes and Fc simultaneously, and thus how surface energy balance nonclosure is linked to underestimated CO 2 exchanges. We show that larger phase differences between large eddies of w and CO 2 concentration (c) are well associated with those between large eddies of w and T (or q), which causes decreased CO 2 uptake as well as decreased H and LE. Section 2 introduces the study site and experimental data. Section 3 briefly describes the methods used in this study, including ensemble empirical mode decomposition (EEMD) and Hilbert transform (HT). Section 4 presents the results, section 5 discusses these results, and section 6 presents some concluding remarks of this study.

Study site
The study site (AmeriFlux site US-Hn1; 46°24′32″ N, 119°16′30″ W) is located within the Hanford Area in central Washington, United States (figure S1 is available online at stacks.iop.org/ERL/14/044016/ mmedia). This area has a semi-arid climate, with a 30 year mean (1986-2016) annual total precipitation of 197 mm (Missik et al 2018). The site is flat and covered with scattered sagebrush and short grasses (figure S1). During the period from February to April when soil water availability is adequate, the vegetation is active with large carbon uptake (Missik et al 2018). The soil texture at the site is sandy with small rocks and gravel interspersed (Gao et al 2017b).

Experimental data
A 10 m EC tower was erected in the field in December 2015. One three-dimensional sonic anemometer (Model CSAT3, Campbell Scientific Inc.) and one open-path gas analyzer (Model LI-7500A, LI-COR Inc.) were installed at 5 m above ground level (a.g.l.) measuring longitudinal, lateral, and vertical wind velocity components (u, v, and w), T, q, and c. These data were logged with a datalogger (Model CR5000, Campbell Scientific Inc.) at a sampling frequency of 10 Hz, and were then processed to obtain half-hourly turbulent fluxes. The procedures of data processing, including despiking, double rotation for the sonic anemometer velocities, and corrections of sonic temperature and density applied to the fluxes, are detailed elsewhere (Liu et al 2009, Gao et al 2016, 2017a. Here, each 30 min averaging interval is referred to as a run. Soil temperature and water content were measured at five depths below the surface (2.5, 5, 10, 15, and 20 cm) using soil temperature probes (Model 109SS, Campbell Scientific Inc.) and volumetric water content reflectometers (Model CS616, Campbell Scientific Inc.). An infrared radiometer (Model SI-111, Apogee Instruments Inc.) was installed on the EC tower at 3.8 m a.g.l. aiming at the soil sensors to measure the soil surface temperature. Soil heat flux was measured at a depth of 5 cm using a self-calibrating soil heat flux plate (Model HFP01SC, Huskeflux Thermal Sensors). These soil data were sampled at 1 Hz and stored as 30 min averages. G 0 was then determined by the calorimetric method (Liebethal et al 2005, Russell et al 2015, Gao et al 2017b combining the measured soil heat flux at 5 cm depth and the change in heat storage in the layer above 5 cm depth. We used two sublayers (0-2.5 cm and 2.5-5 cm) between the surface and 5 cm depth to calculate the change in heat storage. The mean soil temperature and water content for each sublayer were calculated by averaging the corresponding top and bottom measurements for the sublayer. Soil thermal properties (e.g. soil thermal conductivity and diffusivity) were determined by laboratory analysis of soil samples collected from the experimental field (Gao et al 2017b).
R n was measured using a two-component net radiometer (Model CNR2, Kipp & Zonen) located at 8 m a.g.l. on the EC tower, which was replaced on 18 March 2016 with a four-component net radiometer (Model CNR4, Kipp & Zonen). Other instruments relevant to this work included a temperature and relative humidity probe (Model HMP45C, Vaisala) colocated with the EC system, and a cup anemometer and wind vane (Model 03002 Wind Sentry, R M Young) at 10 m a.g.l. on the EC tower. Data from these instruments were also stored as 30 min averages with a sampling frequency of 1 Hz.
The daytime data collected during the months from February to April in 2016 and 2017 with large carbon uptake signals (i.e. Fc<-0.5 μmol m −2 s −1 ) were chosen, including a total of 1948 half-hour runs before the following data selection criteria were applied. Runs were then selected based on the following criteria: (1) Turbulence intensity, I u (I u, u u s = / where u and u s represent the mean and standard deviation of u), was less than 0.5 to assure Taylor's frozen turbulence hypothesis (Stull, 1988). Approximately 11% of the 30 min runs were excluded by this step.
(2) Net radiation (R n ) was larger than 200 W m −2 to ensure that the trends in Fc were not due to changes in photosynthetically active radiation (PAR). Although PAR measurements were not available at this site during these periods in 2016 and 2017, the incoming short-wave radiation exhibited a linear relationship with R n (not shown here). Approximately 28% of the runs left from the previous step were excluded.
(3) Wind direction was between 300°and 60°, under which the influence of the tower structure on the EC system was minimal and the measurements had a large uniform fetch to minimize effects of unsteady flows and horizontal advection on turbulent transport. Approximately 43% of the runs left from the previous step were excluded.
Overall, 344 half-hour runs accounting for about 18% of the total runs were selected and analyzed.

Methodology
We used the EEMD tool (Huang et al 1998, Huang and Wu 2008, Wu and Huang 2009) to detect and extract low-frequency signals (i.e. large eddies) from the 10 Hz turbulence data, as shown in section 3.1. In section 3.2, we introduced the HT (Huang et al 1998, Wu 2008, Wu andHuang 2009) to calculate both the magnitudes and phase differences between the large eddies of w and c, T, and q, respectively.

Ensemble empirical mode decomposition (EEMD)
Turbulence is highly nonlinear and non-stationary in nature, with eddies of different scales super-imposed on each other; these characteristics of turbulence may pose challenges in using Fourier-and wavelet-based methods to extract large eddies. EEMD has two advantages that make it capable of handling such demands (Huang et al 1998, Huang and Wu 2008, Wu and Huang 2009. First, EEMD is adaptive and has high locality without using a pre-determined basis function so that it can handle the nonlinear and nonstationary occurrences in turbulence data. Second, using a sifting process, EEMD decomposes the turbulence data into a limited number of oscillatory components, each with its own distinguishing frequency Wu 2008, Wang et al 2013). Thus, large eddies can be simply extracted as the sum of certain oscillatory components (Gao et al 2017a).
The details of EEMD can be found in Huang et al (1998) and Wu and Huang (2009). Briefly, for a 30 min turbulence time series x (x=w′, c′, T′ and q′, where the prime denotes the turbulent fluctuations relative to the 30 min block average), 13 oscillatory components C k x , (k=1, 2, K, 13) are extracted along with a residual marked as the fourteenth component C , Numerous studies separate large eddies from small eddies by identifying a gap frequency in the velocity spectra or flux cospectra Mahrt 2003, Cava andKatul 2008). However, the spectral gap is difficult to identify on some occasions, especially under unstable atmospheric conditions when convective motions tend to overlap the scales between large eddies and small eddies (Mahrt et al 2001). Based on previous studies and our tests, large eddies, here defined as x C ,

Hilbert transform (HT)
For a time series x t , where P is the Cauchy principle value integral. From the original series and its HT, an analytic signal can be formed using x F ( ) the instantaneous amplitude and phase functions, respectively, for the time series x are given by The Hilbert cross spectrum of two time series x t ( ) and y t ( ) is where * denotes complex conjugation. The instantaneous phase difference (from -90°to 90°) between two time series is calculated using Im HCS Re HCS , 6 xy xy xy where Im and Re are the imaginary and real parts of the Hilbert cross spectrum, respectively. Here, we used the time-averaged absolute phase difference xy F | | to represent phase differences between two times series; i.e. l wc F | | denote the mean absolute phase difference between large eddies of w and c. Separate tests indicate that the calculated phase difference using HT is close to the pre-defined phase difference when sinusoidal signals are processed (figure S5). Given that the sign of the CO 2 flux is opposite to that of H and LE during the daytime, wc F | | refers to the mean absolute phase difference between w and negative c to simplify analysis. As the time-averaged amplitude A x ( ) and standard deviation of a time series x are well correlated

Variations in Fc with surface energy balance closure
For the selected runs, Fc varied from -0.5 to -7.5 μmol m −2 s −1 with negative values denoting carbon uptake by plants ( figure 1). To investigate the impacts of different environmental drivers on variations in Fc, we plotted the distributions of Fc against R , n air temperature, and vapor pressure deficit (VPD) (figures 1(a)-(c)). While Fc had no clear dependence on R n when R n >200 W m −2 , as the air temperature or VPD increased the magnitude of Fc first increased and then decreased when air temperature was higher than 290 K or VPD>10 hPa. As the sum of turbulent heat fluxes H LE + ( )increased, the magnitude of Fc also increased ( figure 1(d)). This pattern was similar to what was observed in the diurnal variations of Fc and turbulent heat fluxes (figure S3). The run-to-run variations in Fc were well correlated with the variations in H and LE (i.e. large carbon uptake and large heat fluxes occurred concurrently) and thus were negatively correlated with variations in the residual of the energy balance closure. Therefore, the magnitude of Fc increased as the surface energy balance closure (CR) increased ( figure 1(e)). Also, within a certain range of air temperature (as indicated by red circles in figure 1(e)), the magnitude of Fc still increased with increasing CR. This result indicates that the increased carbon uptake as CR increased was not likely caused by changes in air temperature or VPD. We found that when Fc was decomposed into fluxes caused by small eddies and large eddies (see section 3.1, figure S4), a large part (∼60%) of the trend of increasing magnitude of Fc with increasing CR was caused by large eddies (figure 1(f)).

Changes in large-eddy flux contributions with phase differences
As large eddies explained a large portion of the variation in Fc with CR, we now investigate mechanisms leading to variations in large-eddy flux contributions for CO 2 and turbulent heat fluxes. Figure 2 and (c)). In addition, a reduction in l wc F | | corresponded to an increase in the CR, while there was no relationship between A l wc and CR ( figure 2(d)). Overall, as the phase differences between large eddies of w and the corresponding scalar decreased, the magnitudes of both the CO 2 and turbulent heat fluxes caused by large eddies increased and thus the energy balance closure was improved. In contrast to flux contributions by large eddies, CO 2 and turbulent heat fluxes contributed by small eddies did not exhibit an obvious trend with the phase difference between small eddies of w and the corresponding scalar (figure S7).

Similarity of large-eddy phase differences between w and different scalars
To explore the similarity in phase differences between large eddies of w and different scalars, we compared the phase differences between large eddies ( figure 3). Given that the phase differences between large eddies of w and scalars were similar, to examine the causes of the varying phase differences we plotted the variations of phase differences between large eddies with the stability parameter ( z L, z = --/ where z and L are the measurement height and Obukhov length, respectively) in figure 4. As the atmospheric instability increased, the phase differences between the different scalars did not exhibit obvious change (albeit with some degree of scatter) ( figure 4(a)), whereas the phase differences between large eddies of w and different scalars decreased ( figure 4(b)).

Discussion
Our results indicate that large eddies which influenced turbulent heat fluxes (i.e. H and LE) also affected Fc in a similar way. The variations in flux contributions by large eddies were highly related to their phase differences (figure 2), consistent with the results in the previous study (Gao et al 2017a). Since underestimation of turbulent heat fluxes is believed to be linked to the influence of large eddies (Gao et al 2016, 2017a, Zhou et al 2018, these large eddies would cause an underestimation of CO 2 fluxes (i.e. decreased magnitudes in CO 2 fluxes). This finding was also supported by the increasing CR (i.e. increasing turbulent heat fluxes) with decreasing phase difference between large eddies of w and c ( figure 2(d)). In summary, large eddies with smaller phase differences have larger flux contributions and thus larger CR, whereas large eddies with larger phase differences have smaller flux contributions and thus smaller CR. Therefore, enlarged phase differences between large eddies linked the surface energy balance non-closure to underestimated CO 2 fluxes.
The similarity of phase differences between large eddies of w and different scalars also reflects two implications for turbulent transport processes at this site. First, the small phase differences among these three scalars (figure 4(a)) indicate that sources and sinks for these three scalars were coherently linked, resulting in their synchronized changes in response to dynamic processes of large eddies. Second, the dynamic processes controlling the turbulent transport of large eddies of T, q, and c were similar for a given phase difference between large eddies of w and scalars. As a consequence, large eddies causing smaller H and LE (thus smaller CR) would also concurrently cause decreased Fc magnitude by enlarging the phase differences between large eddies of w and all three scalars simultaneously. Therefore, at the EC sites where the sources and sinks for T, q, and c are similar, large eddies could cause underestimation of Fc if the turbulent heat fluxes are underestimated due to the influence of large eddies on turbulent processes.
Our results also show that the phase differences between large eddies of w and different scalars decreased as instability increased. One possible explanation is that the variations in phase differences were affected by changes in the structures of large eddies with instability. Li and Bou-Zeid (2011) suggested that large eddies of hairpin vortices and thermal plumes are dominant under near-neutral and free-convection conditions, respectively. Using a quadrant analysis (text S1, figure S8, and Katul et al 1997aKatul et al , 1997b, we computed the flux contribution fraction and time fraction in each quadrant (figure S9). Increasing instability led to an increased imbalance between ejecting and sweeping eddy contributions, consistent with previous results in Li and Bou-Zeid (2011) and Li et al (2018). Also, as instability increased, the absolute flux contributions and time fractions of inward and outward interactions decreased, leading to a decrease in the phase difference between large eddies of w and different scalars. Therefore, the variations in phase differences between large eddies of w and different scalars are most likely caused by changes in the structures of large eddies under increasing instability.

Conclusions
Using EC data from the 2016 and 2017 growing seasons over a sagebrush ecosystem, we analyzed the link between the surface energy balance closure and CO 2 fluxes. Our results show that large eddies explained about 60% of the observed trend of increasing CO 2 flux magnitude with improved energy balance closure. The synchronized variations and small phase difference among large eddies of T, q, and c suggest that the sources and sinks of the scalars were similar at this site. The varying phase differences between large eddies of vertical velocity and all three scalars (T, q, c) were most likely caused by variations in the structures of large eddies with varying instability. Therefore, an enlarged phase difference between large eddies of w and scalars resulted in small flux contributions by large eddies. Our results demonstrate that reduced CO 2 flux magnitudes are mechanistically linked to decreased closure of the surface energy balance.
Improving our understanding of the linkage between energy balance non-closure and CO 2 fluxes is crucial for better interpretation of EC measurements of turbulent heat and carbon exchanges. Although non-closure has been reported at most EC sites, whether or not large eddies would simultaneously cause underestimated turbulent heat and CO 2 fluxes  depends on whether the sources and sinks for T, q, and c are similar. Given widespread observations in nonclosure in the flux community, the processes identified here may bias CO 2 fluxes at many sites and cause upscaled regional and global budgets to be underestimated. Our study presents direct evidence that more studies are needed to investigate how landscape heterogeneity, which may result in dissimilarity in scalar sources and sinks, influences CO 2 fluxes through the influence of associated large eddies on flux exchange.