Detailed descriptions of the study site and long-term flux measurement system can be found in Bi et al. (2007). Only a brief account is provided here.

Preliminary measurements were carried out in May 2004 over a flat, homogeneous, grassland site (area: 300 m × 400 m) in southern China (Bi et al., 2007; Qi et al., 2015). The site was 12.5 m above sea level and was located in the tropical monsoon region (22.43°N, 113.25°E). The air temperature and wind velocity components ( u , v, and w, respectively) were measured using a 3D sonic anemometer (CSAT3, Campbell Scientific Inc., Logan UT, USA). Water vapor density and CO₂ were measured using a LI-7500 system (LiCor Biosciences, Lincoln, NE, USA). Fast response sensors were installed 3.9 m above the ground, with a sampling frequency of 10 Hz. Other supporting data, such as soil heat flux, soil temperature, and radiation, were measured in the experiment (Bi et al., 2007); however, these measurements were not included in our analysis.

Following Foken et al. (2004), post-field data quality control processes were applied. In particular, noise and various kinds of interference from 30-min measurements of turbulence using a criterion of X t < X ¯ − 4 σ or X t < X ¯ + 4 σ were eliminated, where X(t) denotes the measurement (i.e., wind speed components and temperature), X ¯ is the mean over the interval, and σ the standard deviation. Data during and after rain events were removed because the sonic anemometer would be in error in these cases. After quality controlling of the 40 days’ data, 1733 out of 1920 half-hour data remained.

The study period was 40 days from 1 June 2004 to 10 July 2004 to re-estimate the sensible heat flux. Three 30-min runs were conducted daily to obtain data under different stability conditions: 3:00–3:30 for stable conditions, 11:00–11:30 for unstable conditions, and 19:00–19:30 for neutral conditions. The mean meteorological conditions for the three runs are shown in Table 1, where the friction velocity u * , sensible heat flux H t (calculated from the conventional EC method), and stability parameter ζ are estimated in Eq. 1: u * = w ′ u ′ ¯ 2 + w ′ v ′ ¯ 2 1 4 , H t = ρ C p w ′ T ′ ¯ , ζ = z / L = − κ z g w ′ T ′ ¯ + 0.61 T w ′ q ′ ¯ / T ¯ u * 3 , (1) where ρ , C p , and g represent air density (kg m⁻³), the specific heat under constant pressure (J kg⁻¹ K⁻¹), and the gravitational acceleration, respectively. L = − T ¯ u * 3 / κ g w ′ T ′ ¯ + 0.61 T w ′ q ′ ¯ is the Obukhov length, which includes the buoyancy correction due to water vapor; κ is the von Karman constant, and L is the height at which the buoyant force begins to predominate over the shear force. Notably, all scalar fluctuations in (1) are derived from the time moving-averaging operation.

The heat transferred by eddies cross a horizontal surface in the atmosphere, at any height z, is given by Priestley and Swinbank (1947): H = ρ ¯ C p w T − T 0 , ¯ (2) where ρ ¯ and C p are the mean air density (kg m⁻³) and the specific heat capacity of dry air at a constant pressure (J kg⁻¹ K⁻¹), respectively. Here, w is the vertical wind velocity and T − T 0 is the temperature change with respect to the environmental temperature (Priestley and Swinbank, 1947) or the base temperature T 0 (Webb et al., 1980). Webb et al. (1980) describe the base temperature T 0 as “…taken as constant at any given height, representing roughly an assumed initial base temperature from which each element of air is warmed (or cooled) during the vertical transfer of heat supplied (or removed) at the underlying surface.” Henceforth, we use the term “environmental temperature” for T 0 . If we substitute the mean temperature T ¯ into Eq. 2, it becomes H = ρ ¯ C p w ¯ · T ¯ − T 0 + ρ ¯ C p w ′ T ′ ¯ = ∆ H + H t , (3) where T ′ represents the temperature fluctuation with respect to the mean temperature T ¯ . Additionally, we replace T ¯ − T 0 with ∆ T .

In the conventional approach, the first term on the right side of Eq. 3 is usually neglected as small values of w ¯ exit in the ASL (Webb, 1982) and the difference between the mean and environmental temperature is considered minimal (Webb et al., 1980). However, we demonstrate in the following that this is not the case, as indicated by Mauder et al. (2008, 2020).

First, consider the two terms on the right side of Eq. 3 denoted as ∆ H and H t . The additional flux ∆ H is ignored by the conventional EC method that considers only H t , which is determined by the product of the mean vertical velocity and the difference between the mean and environmental temperatures; it does not account for the additional heat transport contributed by anisotropic thermal structures in the ASL. The decomposition of the sensible heat flux in Eq. 3 suggests that the heat flux H is determined by two temperature scales and two velocity scales. It should be noted that both ∆ H and H t are generated by atmospheric eddy motion, regardless of the mean flow and velocity fluctuations (Shang et al., 2004).

To more accurately illustrate ∆ T for the measurements, detailed information about the three runs with different stability conditions is described. Table 1 lists the mean meteorological conditions. The stable condition run had a relatively small horizontal wind speed (1.18 ms⁻¹), a small u * (0.12 ms⁻¹), and small σ T ′ (0.11 ℃ ) and σ w ′ (0.18 ms⁻¹). Under unstable conditions, a moderate horizontal wind speed (3.23 ms⁻¹) and a large u * (0.38 ms⁻¹), σ T ′ (0.59 ℃ ), and σ w ′ (0.47 ms⁻¹) were observed. For the neutral condition, the horizontal wind speed was moderate at 3.53 ms⁻¹; large values were recorded for u * (0.36 ms⁻¹) and σ w ′ (0.50 ms⁻¹), yet σ T ′ was small (0.06 ℃ ). All incoming winds of the three runs had nearly the same direction.

Figure 3 shows the PDFs for the temperature fluctuations of the three runs under the different stability conditions. To contrast them more fully, all maximums of the PDFs were normalized to the value of 1. Distinct peaks were found in the PDF distributions for all stability conditions. For the three runs, the differences between the mean and environmental temperatures ( ∆ T ) were −0.07 °C, 0.29 °C, and 0.01 °C for stable, unstable, and neutral conditions (Table 2), respectively. To understand the dynamic mechanism for ∆ T , detailed information including the skewness of temperature ( S T = T − T ¯ 3 / T − T ¯ 2 3 / 2 ) and vertical velocity ( S w = w − w ¯ 3 / w − w ¯ 2 3 / 2 ) and the skewness of the temperature derivative ( S T ′ = ∂ T / ∂ t 3 / ∂ T / ∂ t 2 3 / 2 ) is also reported in Table 2. Under stable conditions ( ζ = 0.147 ), both S T and S w were greater than 0. However, under unstable conditions ( ζ = − 0.109 ), negative values were found for both S T and S w . It is known that the skewed distributions of temperature and vertical velocity indicate the nature of the upward thermal structures in the ASL (Chu et al., 1997). Here, the thermal structures under stable and unstable conditions were generated by the instabilities of the velocity layer for S T ∙ S T ′ < 0 (Belmonte and Libchaber, 1996). The vertical movement of thermal structures led to the existence of a non-zero ∆ T . For example, upward warm and cold structures create temperature differences of ∆ T > 0 and ∆ T < 0 , respectively. Good agreement was attained between ∆ T and the stability conditions such that ∆ T < 0 for stable conditions and ∆ T > 0 for unstable conditions. Neutral conditions resulted in a ∆ T value of approximately 0.

Figure 4 shows the variations in ∆ T for the period 01 June to 10 July 2004. ∆ T had a maximum value on 04, 12, and 27 June. Generally, positive values of ∆ T were observed in the daytime hours; at night, ∆ T was smaller and negative on average (Figure 5). Over the observational period, under the condition of strong solar forcing on sunny days, large ∆ T ( > 0.8°C) was observed on 4 of the 40 days: 04, 12, 14, and 27 June. On 6 June it was rainy, so ∆ T remained small ( 0.19 ℃ ). The diurnal variation of average ∆ T was consistent with the diurnal variation in the solar radiation: both exhibited strong solar forcing and positive ∆ T during the day, reached a maximum at noon, but then declined to near 0 or below with negative ∆ T values at night. The relationships between the stability parameter ζ and ∆ T for 40 days are shown in Figures 5B, C. Because the unstable boundary layer was more conducive to the development of thermal structures, large positive values of ∆ T were usually found under unstable conditions in the daytime. In contrast, the stratification of the boundary layer at night prevented heat and mass exchange between the Earth’s surface and upper air layers, such that small ∆ T dominated.

Wind data from sonic anemometers are usually transformed into a mean streamline parallel coordinate system to correct for non-zero w ¯ effects (Wilczak et al., 2001; Finnigan et al., 2003; Dellwik et al., 2010). Following Wilczak et al. (2001), the planar-fit technique was adopted in this study. The corrected vertical velocity can be calculated using Eq. 8: w r = w − b 0 − b 1 u ¯ − b 2 v ¯ , (8) where: b 0 is a mean zero offset, possibly due to electronic problems and flow perturbation in the measured vertical velocity; b 1 and b 2 are wind direction-dependent coefficients; and u ¯ and v ¯ are the longitudinal and lateral velocity, respectively, as measured by a sonic anemometer.

The wind data for 27 to 30 June were selected to determine the linear regression coefficients ( b 0 , b 1 , and  b 2 when the sonic anemometer was functioning properly; a wide range of wind directions were found for these four days. The resulting b 0 , b 1 , and  b 2 were 0.0132 m s⁻¹, 0.1197, and 0.0156, respectively.

Figure 4 shows relatively small changes in the vertical velocity ( < 0.1 ms⁻¹); notably, the lack of diurnal variation in w ¯ dominated our measurements. The grasslands had a slight effect on the wind direction, depending on the vertical velocity. However, no distinct vertical component was found even at midday when vertical convection is more likely.

There was a linear relationship between the daily maximum of ∆ T and the corresponding local vertical velocity (Figure 6). In general, ∆ T increased with w ¯ . Since ∆ T was generated by the thermal structures, it can be deemed an indicator of the activity level of those structures in the ASL. The relationship between the maximum ∆ T and w ¯ suggests that an increase in the local w ¯ can result in more active thermal structures. In turn, the development of the thermal structure can be suppressed by terrain gradients, large eddies in the ASL, and/or by local thermal circulations, regardless of whether the local w ¯ is induced.

The measurement location was in the tropic monsoon region of southern China. Because measurement commenced in June, the air temperature was relatively high and the diurnal variation in the maximum temperature ranged between 26 ℃ and 38 ℃ (figure not shown). Figure 7 shows the sensible heat flux over the 40-day period. The maximal conventional EC method was used to calculate H t ; most values were larger than 80 Wm⁻² at noon, with 04, 14, 17, and 27 June showing a clear sky and strong solar forcing that resulted in flux H t maxima of 231.4, 232.4, 236.3, and 229.8 Wm⁻², respectively. Relatively frequent rainfall was found within the observation period, and the missing results of H t were largely induced by rainfall events. Notably, rainfall at noon on 06 June produced a low H t value of 32.6 Wm⁻².

Large differences (>50 Wm⁻²) between the new method used to estimate the total sensible heat flux H and the conventional EC method used to calculate H t were found for eight days: 03, 04, 13, 14, 15, 25, 26, and 27 June; ∆ H for those days were 78.3, 63.4, 66.6, 61.2, 93.5, 51.1, 57.4, and 59.9 Wm⁻², respectively, with the ratio of the ∆H to H t reaching an astonishing 50.7%, 29.1%, 42.2%, 30.6%, 43.3%, 44.5%, 29.6%, and 26.1%, respectively. Those large differences reached a maximum at noon when strong solar forcing occurred under a large ∆ T of 0.77, 1.02, 0.39, 0.54, 0.42, 0.48, 0.59, and 0.80 ℃ , respectively. At the same time, a relatively small vertical convection with positive local vertical velocities of 0.07, 0.04, 0.06, 0.08, 0.11, 0.08, 0.07, and 0.07 m s⁻² was observed on these eight days. Notably, on occurrences with large ∆ H , the wind speeds were not large, ranging 1.6 to 3.9 m s^-1, and winds were mostly from the south (Figure 8). In addition to those eight days, additional flux less than 50 Wm⁻² but larger than 30 Wm⁻² was observed for an additional 11 days: 01, 05, 11, 12, 16, 22, 23, and 24 June, and 01, 06, and 08 July.

In addition to the large positive additional flux observed, there were also relatively large negative additional flux values (<−30 Wm⁻²) for some of the other days: −44.0, −32.0, −30.1, −48.2, −33.4, −39.4, −54.9, and −34.3 Wm⁻² for 01, 02, 08, 19, 21, 22, and 28 June, and 08 July, respectively. As with the positive flux, all of the negative flux occurred at noon, with ∆ T ranging between 0.40 ℃ and 0.68 ℃ , but with small local subsidence in the vertical velocities ranging between −0.07 and −0.04 m s⁻¹.

Mauder and Foken (2006) indicated an uncertainty of 5% or 10 Wm⁻² from sonic anemometer measurements, using the conventional EC method to estimate the sensible heat flux. The large discrepancies described in sensible heat flux calculated from the conventional EC method and the new EC method at noon in our observations seem plausible as the thermal structures were active at this time of day. The presence of small local vertical velocities resolves these large discrepancies.

The conventional EC method was developed under the condition of homogeneous and isotropic turbulence flow, and it was thought able to capture the heat flux contributed by all signals of small-scale motions. The visible contradiction between the prerequisite of the conventional EC method and the widespread anisotropic turbulence in the ASL was the first indication that the conventional EC method may fail to estimate the correct heat flux. The identified thermal structures were responsible for the majority of the vertical heat, momentum, and mass fluxes in the ASL (Gao et al., 1989; Chen et al., 1997). Behind these thermal structures are anisotropic flows. The ability of an anisotropic flow to transport heat differs considerably from that of an isotropic flow.

We obtained the local environmental temperature from the ramp-like structures in our model (Figure 2). If warm (cool) structures occupy the space, then the time mean temperature will be higher (lower) than the environmental temperature. The time scales of the thermal structures are much less than 30 min; thus, they can be captured by the 30-min average operator. Additionally, their spatial scales are not very large, such that the results represent the local heat flux. To measure the additional heat flux contributed by large-scale organized structures, Mauder et al. (2008) used the time–space-averaged temperature as the environmental temperature. However, this may be a limitation of their theory in practice.

When using the most probable temperature within a 30-min period as the environmental temperature, an additional flux due to the transportation of thermal structures could be detected. The additional vertical sensible eddy heat flux ∆ H is partly decided by the difference between the mean and environmental temperatures, ∆ T . The non-zero ∆ T presented here was determined by similar existing thermal structures in Figure 2, meaning that both ∆ H and H t are a result of the vertical transportation of eddies (Shang et al., 2004). Moreover, even ∆ H has an inseparable relationship with local vertical convection. A good linear relationship between heat flux w ′ T ′ ¯ and ∆ T was found (Figure 9B), such that a larger sensible heat flux H t ( H t = ρ C p w ′ T ′ ¯ ) occurred with larger ∆ T , while smaller ∆ T corresponded to a smaller H t . This linear relationship is also evident when comparing Figures 5A, 9A as both ∆ T and H t have the same diurnal variation. The identified thermal structures are responsible for the majority of vertical heat flux in the ASL (Gao et al., 1989; Chen et al., 1997). Large H t means that the thermal structures are active and an apparent ∆ T emerges. In these active thermal events, additional fluxes are visible if non-zero w ¯ occurs. On the other hand, no additional flux exists in the absence of thermal structures, even with strong local vertical convection. Unlike the similar approaches proposed by Lee (1998) and Mauder et al. (2008), in which the former quantifies the convective sensible heat flux by measuring the transported local vertical gradient of the scalar while the latter undertake this by measuring the difference between the time-averaged scalar and the time–space-averaged scalar, ∆ H and H t cannot be studied separately. They are both part of the turbulent flux but can be captured by a single-point tower measurement.

The linear regression relationship from Figure 9 provides a modeled estimation of ∆ T from the conventional EC method w ′ T ′ ¯ ( H t ): ∆ T = 0.0061 + α w ′ T ′ ¯ , (9) where α is equal to 3.55. The coefficient of determination ( R 2 ) of the linear regression is 0.77. Here, w ′ T ′ ¯ is the dimensionless operator, and the coefficient of 3.55 is determined by the thermal structures. The first term on the right side of Eq. 9 will be discarded because zero heat flux exchange means that no thermal structures appeared. We should note that α varies at different underlying surfaces, but the linear relationship in Eq. 9 is reliable.

Substituting Eq. 9 into Eq. 3 leads to the following relationship: H m o d e l = 1 + α w ¯ H t , (10)

Based on Eq. 10, it is clear that, once the coefficient is known, the additional flux coinciding with local vertical convection can be resolved. For example, an additional flux of 35.5% of H t should be included if a slight subsidence of w ¯ = 0.1   m s − 1 existed for our measurements.

The additional vertical sensible eddy heat flux ∆ H is not only decided by the difference between the mean and environmental temperatures ∆ T but is also sensitive to the local mean vertical wind speed w ¯ . ∆ T is usually greater than zero in the daytime (Figure 5A); thus, the local subsidence or the local uplift would lead to an underestimation or overestimation of the total sensible heat flux (Figure 10). However, because only small ∆ T occurs at night, there will be no large underestimation or overestimation. This case somewhat differs from the “cool down” events described by Williams and Hacker (1992) or McNaughton (2004) that underestimate the sensible heat flux for downdrafts as cooler than the ambient air (Mauder et al., 2008). Variations in the vertical velocity perturbation of more than one order larger than w ¯ lead to a slight subsidence that cannot prevent warm structures from propagating upward due to the instability of the velocity layer as opposed to buoyancy. Closer examination of Figure 5 shows that, although there was slight subsidence at noon on more than 10 of the 40 days, all ∆ T were positive.

It is known that the majority of vertical heat, momentum, and mass fluxes in the ASL are due to the movements of the identified thermal or coherent structures (Qiu et al., 1995; Chen et al., 1997). The creation of these thermal structures is dynamically dominated by the vertical shear of wall-bounded turbulent flows (Gao et al., 1989), and their thermodynamics are dominated by the temperature difference between the upper and bottom boundaries of Earth’s surface or the Rayleigh–Bénard convection (Shang et al., 2004). The structures are usually more pronounced under unstable conditions (Figure 1). Similar ramp patterns of humidity and CO₂ can also be found, and they are both shaped by the turbulence structures in the ASL. Therefore, additional heat flux in the form of latent heat and additional CO₂ transport due to this turbulence should also be considered.

In this study, both the conventional EC estimate and the additional contribution are fluxes transported by the turbulence in the surface layer; that is, the additional contribution is a component of heat transported by the thermal structures but cannot be identified by the conventional EC method. A more intuitive understanding comes from a laboratory experiment of Shang et al. (2004), who compared the heat transports calculated by the EC method and those injected by heater in a turbulent Rayleigh–Bénard convection system. They concluded that the conventional EC estimates can significantly underestimate the injected heat flux under convective conditions, and the underestimated heat transport was captured by the new method (Eq. (3)). The results revealed by Shang et al. (2004) may be inspiration for addressing the energy imbalance problem in the terrestrial surface. The energy imbalance problem was widespread at FLUXNET sites and was most apparent during the day (Wilson et al., 2002). According to Wilson et al. (2002), the difference between the available energy and turbulent energy fluxes reaches its maximum value at noon. Foken et al. (2010) also reported an imbalance of 20%–30% during the daytime in the LITFASS-2003 experiment. Based on the results of our study, the additional contribution ∆ H , which is usually ignored and also usually reaches its maximum value at noon, may be one of the main reasons for the energy imbalance.

The new estimate considers heat transport by thermal structures. These are not unique to the underlying surface of grassland types but are rather a common feature shared by all types of underlying surfaces, indicating the universality of this research. The contribution of thermal structures generated by plant canopies such as grasslands, forests, and crops to heat transport and their related boundary layer energy imbalance problem have received widespread attention (Chen et al., 1997; Castellvi et al., 2008; Holwerda et al., 2021; Cely-Toro et al., 2023). As they are all based on the conventional EC results to study the contribution of thermal structures to heat flux, there is also the inevitable issue of the energy imbalance problem.

Using the conventional EC method, the large-eddy simulation study of the energy imbalance problem by Steinfeld et al. (2007) showed that the imbalances are relatively small (∼5%) at a height of 20 m. However, reports from field experiments indicated a larger imbalance (10%–30%) at lower measurement heights closer to the surface (Wilson et al., 2002; Oncley et al., 2007; Foken et al., 2010). No viable reasons were offered to explain the imbalance difference between the different measured heights.

We describe the presence of thermal structures in the ASL that determine the existing ∆ T . Those structures, without mixing with the surrounding air, create the ramp-like patterns in the model, leading to temperature skewness (Table 2). When the thermal structures propagate into the upper air and mix with the surrounding environment, the thermal structures gradually fade, and a corresponding decrease in ∆ T occurs. Thus, the value of ∆ H would also decrease with height while having the same w ¯ at different heights. This means that, if the ignored ∆ H is included in the total vertical sensible heat flux, the question of the reduction in the imbalance rate with height may be partly answered. Of course, further study is needed to confirm this premise.