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

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).

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, Oncley et al 2007, 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) co-located 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⁻² s⁻¹) 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}={\sigma }_{u}/\bar{u},$ where $\bar{u}$ and ${\sigma }_{u}$ 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⁻² 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.

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 non-stationary 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 (Huang and 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, ..., 13) are extracted along with a residual marked as the fourteenth component $\left({C}_{14,x}\right),$ and hence $x=\displaystyle {\sum }_{k=1}^{14}{C}_{k,x}$ (Wang et al 2013, Gao et al 2017a). The residual is included to ensure that the accumulated fluxes from the extracted components are comparable to the calculated fluxes from the traditional data processing (figure S2).

Numerous studies separate large eddies from small eddies by identifying a gap frequency in the velocity spectra or flux cospectra (Vickers and Mahrt 2003, Cava and Katul 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}_{l}=\displaystyle {\sum }_{k=6}^{13}{C}_{k,x},$ are mostly responsible for the run-to-run variations in fluxes (figures S2–S4). Small eddies are defined as ${x}_{s}=\displaystyle {\sum }_{k=1}^{5}{C}_{k,x},$ and their flux contributions present a smooth diurnal variation (figure S4).

For a time series $x(t),$ its HT $\hat{x}(t)$ is

$$\begin{eqnarray}&&\hat{x}\left(t\right)=\displaystyle \frac{1}{\pi }P\displaystyle {\int }_{-\infty }^{\infty }\displaystyle \frac{x\left(\tau \right)}{t-\tau }d\tau ,\end{eqnarray} \tag{ 1 }$$

where P is the Cauchy principle value integral. From the original series and its HT, an analytic signal can be formed using

$$\begin{eqnarray}&&{z}^{x}\left(t\right)=x\left(t\right)+i\hat{x}\left(t\right)={A}^{x}\left(t\right){e}^{i{{\rm{\Phi }}}^{x}\left(t\right)},\end{eqnarray} \tag{ 2 }$$

where $i=\sqrt{-1},$ and ${A}^{x}\left(t\right)$ and ${{\rm{\Phi }}}^{x}\left(t\right),$ the instantaneous amplitude and phase functions, respectively, for the time series $x$ are given by

$$\begin{eqnarray}&&{A}^{x}\left(t\right)={\left({x}^{2}+{\hat{x}}^{2}\right)}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}^{x}\left(t\right)={\tan }^{-1}\left(\displaystyle \frac{\hat{x}}{x}\right).\end{eqnarray} \tag{ 4 }$$

The Hilbert cross spectrum of two time series $x(t)$ and $y(t)$ is

$$\begin{eqnarray}&&{\rm{H}}{\rm{C}}{{\rm{S}}}^{xy}\left(t\right)={z}^{x* }\left(t\right)\cdot {z}^{y}\left(t\right),\end{eqnarray} \tag{ 5 }$$

where * denotes complex conjugation. The instantaneous phase difference (from –90° to 90°) between two time series is calculated using

$$\begin{eqnarray}&&{{\rm{\Phi }}}^{xy}\left(t\right)={\tan }^{-1}\left(\displaystyle \frac{\text{Im}\left[{\rm{H}}{\rm{C}}{{\rm{S}}}^{xy}\left(t\right)\right]}{\mathrm{Re}\left[{\rm{H}}{\rm{C}}{{\rm{S}}}^{xy}\left(t\right)\right]}\right),\end{eqnarray} \tag{ 6 }$$

where ${\rm{I}}{\rm{m}}$ and ${\rm{R}}{\rm{e}}$ are the imaginary and real parts of the Hilbert cross spectrum, respectively. Here, we used the time-averaged absolute phase difference $\left|{{\rm{\Phi }}}^{xy}\right|$ to represent phase differences between two times series; i.e. $\left|{{\rm{\Phi }}}_{l}^{wc}\right|$ 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₂ flux is opposite to that of H and LE during the daytime, $\left|{{\rm{\Phi }}}^{wc}\right|$ refers to the mean absolute phase difference between w and negative c to simplify analysis. As the time-averaged amplitude $\left(\overline{{A}^{x}}\right)$ and standard deviation of a time series x are well correlated (figure S6), we used ${{\rm{A}}}^{xy}=\overline{{A}^{x}}\cdot \overline{{A}^{y}}$ to represent the magnitude of two time series.

For the selected runs, Fc varied from –0.5 to –7.5 μmol m⁻² s⁻¹ 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⁻², 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 $\left(H+LE\right)$ 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)).

Zoom In Zoom Out Reset image size

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₂ and turbulent heat fluxes. Figure 2(a) shows the variations of CO₂ fluxes contributed by large eddies $\left(F{c}_{l}\right)$ with the phase differences $\left(\left|{{\rm{\Phi }}}_{l}^{wc}\right|\right)$ between large eddies of $w\left({w}_{l}\right)$ and c (c_l). In general, the magnitude of $F{c}_{l}$ decreased with increasing $\left|{{\rm{\Phi }}}_{l}^{wc}\right|.$ To illustrate the influence of the magnitudes $\left({{\rm{A}}}_{l}^{wc}\right)$ of ${w}_{l}$ and ${c}_{l}$ on their flux contributions, we also divided the selected runs into six groups based on ranges of ${{\rm{A}}}_{l}^{wc}.$ Given the same $\left|{{\rm{\Phi }}}_{l}^{wc}\right|,$ the absolute values of $F{c}_{l}$ were larger for larger ranges of ${{\rm{A}}}_{l}^{wc},$ whereas for different groups, the absolute values of $F{c}_{l}$ were larger for smaller $\left|{{\rm{\Phi }}}_{l}^{wc}\right|.$ The variations in H and LE contributed by large eddies (i.e. ${H}_{l}$ and $L{E}_{l}$) with their corresponding phase differences (i.e. $\left|{{\rm{\Phi }}}_{l}^{wT}\right|$ and $\left|{{\rm{\Phi }}}_{l}^{wq}\right|$) followed the same patterns as those in $F{c}_{l}$ (figures 2(b) and (c)). In addition, a reduction in $\left|{{\rm{\Phi }}}_{l}^{wc}\right|$ corresponded to an increase in the CR, while there was no relationship between ${{\rm{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₂ 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₂ 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).

Zoom In Zoom Out Reset image size

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). Clearly, $\left|{{\rm{\Phi }}}_{l}^{wc}\right|,$ $\left|{{\rm{\Phi }}}_{l}^{wT}\right|,$ and $\left|{{\rm{\Phi }}}_{l}^{wq}\right|$ were quite similar to each other across different ranges of phase differences. $\left|{{\rm{\Phi }}}_{l}^{wc}\right|$ and $\left|{{\rm{\Phi }}}_{l}^{wT}\right|$ had a larger ${R}^{2}$ value than $\left|{{\rm{\Phi }}}_{l}^{wc}\right|$ and $\left|{{\rm{\Phi }}}_{l}^{wq}\right|$ (i.e. 0.79 versus 0.56).

Zoom In Zoom Out Reset image size

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 ($\mbox{--}\zeta =-z/L,$ 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)).

Zoom In Zoom Out Reset image size