Linkage between surface energy balance non-closure and horizontal asymmetric turbulent transport

A number of studies have reported that the traditional eddy covariance (EC) method generally underestimated vertical turbulent fluxes, leading to an out-standing non-closure problem of the surface energy balance (SEB). Although it is recognized that the enlarged surface energy imbalance frequently coincides with the increasing wind shear, the role of large eddies in affecting the SEB remains unclear. On analyzing data collected by an EC array, considerable horizontal inhomogeneity of kinematic heat flux is


INTRODUCTION
As a fundamental concept in the theoretical description of the global climate system, the surface energy balance (SEB), which depicts the energy exchange processes at the Earth's surface, has been extensively investigated for decades.The eddy covariance (EC) methodology is widely employed to directly measure turbulent heat fluxes, soil heat flux, and net radiation, supporting the assessment of the SEB on the local ecosystem scale (Baldocchi et al., 2001;Mauder et al., 2007;Oncley et al., 2007;Liu et al., 2021;Metzger et al., 2021).The closure of the SEB is expected under the idealized circumstance that fluxes are perpendicular to the horizontally uniform two-dimensional exchange surface without a canopy (Mauder et al., 2020).However, it has been recognized that, at most EC sites, the sum of sensible (H) and latent heat fluxes (LE) is, to some degree, lower than the available energy (R n − G 0 , where G 0 and R n represent soil heat flux and net radiation respectively), raising an outstanding non-closure problem of the SEB (Foken, 2008;Mauder et al., 2013Mauder et al., , 2020;;Stoy et al., 2013;Gao et al., 2017;Reed et al., 2018;Liu et al., 2021).Such observed SEB non-closure poses challenges to the calibration and validation of land-surface models that fulfill the aforementioned ideal assumptions (Jung et al., 2009;Williams et al., 2009).Therefore, how to address the discrepancies in SEB closure between observations and models has been a subject of active research.
Several reasons have been identified to explain the observed SEB non-closure.Mauder et al. (2020) classify these reasons into four main aspects: instrumentation errors (e.g., transducer shadowing; Horst et al., 2015), data-processing errors (e.g., inadequate averaging and choices of flux correction methods; Finnigan et al., 2003;Fratini and Mauder, 2014), underestimated or ignored contributions of heat storage (e.g., underestimation of heat storage underground and in the canopy; Leuning et al., 2012;Higgins, 2012), and sub-mesoscale transport processes (e.g., secondary circulations and large-scale organized structures; Kanda et al., 2004;De Roo and Mauder, 2018;Morrison et al., 2022).Given the critical role in transporting energy and scalars in the atmospheric boundary layer (ABL), the influence of large-scale coherent eddies on SEB non-closure has drawn extensive attention (Foken et al., 2011;Zhou et al., 2018).In the unstable ABL, large eddies that are either generated aloft or developed from the heated surface frequently disturb the ABL and modulate structures of local turbulent eddies, thus leading to enhanced or decreased fluxes (McNaughton, 2004;Zhang et al., 2010).Using measurements over a vegetative canopy with large evapotranspiration, it is revealed that enlarged phase differences between signals of vertical velocity and water vapor associated with large eddies account for the decrease in LE (Gao et al., 2017).Therefore, eliminating such phase difference increases the LE, leading to an improved SEB closure over a surface characterized by a low Bowen ratio.
Although the SEB closure can be improved by artificially synchronizing such a phase shift, the underlying mechanisms regulating the variation of the SEB closure ratio are worth further investigation (Stoy et al., 2013;McGloin et al., 2018;Gao et al., 2020;Liu et al., 2021).The convection process is proposed to explain the concurrence of decreased SEB closure ratio and strong surface heating.Spectral analysis indicates that turbulent eddies with vertical scales much larger than twice their observation height are energized due to enhanced thermal convection (Gao et al., 2020).The flux contribution of sweeps and ejections related to these large eddies becomes significantly imbalanced under strongly convective conditions, coinciding with the observed increase in SEB non-closure with enhancing instability (Liu et al., 2021).Apart from single-tower measurements, a spatial EC tower set-up has been employed to investigate the spatial variability of turbulent fluxes and its influence on SEB non-closure (Morrison et al., 2021(Morrison et al., , 2022)).Katul et al. (1999) observed that H is spatially relatively homogeneous compared with LE above a pine forest, akin to the results of Christen and Vogt (2004) that reveal large eddies contribute significantly to the momentum dispersive flux but less to H. On the contrary, both airborne and multi-tower EC measurements suggest that H determined by the spatial EC method is systematically higher than the counterpart derived solely from the time-average method, indicative of the spatial flux contribution to the SEB closure (Mauder et al., 2007(Mauder et al., , 2008;;Engelmann and Bernhofer, 2016).Although these findings provide direct evidence that the missing flux contributions related to large eddies behave as a significant reason for the SEB non-closure, it remains lacking how such large eddies modulate those fluxes, which is the compass of the work herein.
In this work, data collected using a spatial EC array set-up are analyzed, aiming to (a) evaluate how SEB closure can be improved by utilizing a combined EC method proposed by Engelmann and Bernhofer (2016), and (b) explore the role that large eddies have in modulating the spatial flux fluctuation, and thus influencing the SEB closure.The paper is organized as follows: Section 2 introduces information on the observation site, measurement set-up, and post-field-data processing procedures; results and conclusions are presented in Section 3 and Section 4 respectively.

2
EXPERIMENT, DATA, AND METHODOLOGIES

Experimental site, measurement set-up, and post-field-data processing
The observation site, as a part of the Guangdong Provincial Observation and Research Station, is located in a subtropical suburban environment (23.20 • N, 113.49• E, 30 m above sea level).The measurement, conducted during an intensive observation period from September 5 to October 31 in 2018, is a part of the long-term land-atmosphere exchange observation project in the Pearl River Estuary area (Lan et al., 2022;Wang et al., 2022).The experimental site is surrounded by a diversity of land types (e.g., rural roads in the east, forested hills in the west, sparse short buildings to the north, and a small parking lot to the south) that provide a complex underlying surface condition (Figure 1a).Within such a forested canopy, an EC array was set up over a horizontal area that is generally flat and covered by continuous grass cover with a height of ∼20 cm.The area of 15 × 21 m 2 was covered by nine three-dimensional sonic anemometers (WindMaster; Gill Instruments, Lymington, UK) that were horizontally separated (Figure 1b,c; the distance between two sonic anemometers was 10.3 m and 7.3 m in the north-south and west-east directions respectively).Each sonic anemometer was mounted on a 2-m tower and was aligned carefully vertically and to the north.Before the intensive observation, all sonic anemometers were factory calibrated and inter-compared.The climatological flux footprint is similar across the nine sonic anemometers, indicative of 90% of the source area, which covers the fetch of order of 50 m (Supporting Information Figure S1).Considering that the EC array is within the forest canopy (Figure 1b), flux contributions from the nearby tall trees via the transport associated with large eddies might be captured by the EC array.Along with our EC array, a long-term operated EC tower (EC main tower) is located to the northeast, providing energy flux measurements (e.g., H, LE, R n , and G 0 ) in such a typical subtropical monsoon climate area.The water vapor density  v measured by the EC main tower was used to derive actual air temperature from sonic temperature measured by the EC array.Detailed instrumentation information is provided in Table 1.The measured time-series from the EC array and the EC main tower were sampled at 10 Hz using Campbell Sci CR3000 dataloggers that were automatically synchronized every 24 hr.
The 10 Hz raw data were processed by a post-field-data program to determine 30-min time-averaged statistics   (Lan et al., 2018(Lan et al., , 2022)).Briefly, the data-processing procedures include the following steps: (a) time-series de-spiking based on the median absolute deviation technique (Mauder et al., 2013) et al., 1980;Schotanus et al., 1983;Liu et al., 2001); and (f) determining quality flags for turbulent fluxes (Foken et al., 2005).In the following sections, 30-min statistics with data quality flag of 0 are used for exploring the general time variations (2,644 data segments of length 30 min are used).Among these, data segments with turbulence intensity (i.e., I u =  u ∕u, where u and  u are respectively the mean and standard deviation of the horizontal velocity aligned to the mainstream) larger than 0.5 or H <10 W⋅m −2 are further excluded to ensure the validity of Taylor's hypothesis and well-resolved energy fluxes by the instruments.Hence, 401 data segments of length 30 min (16% of the entire dataset) are used for spectral analysis and scaling argument.It is worth noting that the planar fit method was employed to perform tilt correction and transform the wind vector into a streamline-parallel coordinate system (i.e., u is aligned to the streamline; v and w are perpendicular to the streamline and the ground surface respectively).Compared with the double rotation used in other studies (McMillen, 1988;Engelmann and Bernhofer, 2016), the planar fit method is not only able to correct the vertical velocity influenced by the inevitable inclination of the sonic anemometers but also reveal the actual vertical wind velocity (Mauder et al., 2008).Supporting Information Table S1 provides the rotation angles ( = 0.0100,  = 0.1060) and the offset of the vertical velocity (−0.0208) derived from the entire observation period.Since the velocity field measured by each sonic anemometer was rotated using the same tilt coefficients, it is reasonable to cross-compare the velocity field and fluxes within the EC array.

Temporal, spatial, and combined EC method
As the cornerstone for the EC technique and flux calculation, the Reynolds decomposition can be performed upon either temporal or spatial averaging.Hereafter, for a time-series of an arbitrary variable a(t), the temporal and spatial averaging are expressed as a(t) and ⟨a(t)⟩ respectively.Correspondingly, the Reynolds decompositions based on temporal and spatial averaging are denoted as a ′ (t) (i.e., a ′ (t) = a(t) − a(t)) and a ′′ (t) (i.e., a ′′ (t) = a(t) − ⟨a(t)⟩) respectively.Hence, the kinematic heat fluxes derived from temporal and spatial EC methods are where M and N represent the number of data records within an averaging interval (i.e., M = 18,000 for a 30-min window and 10 Hz sampling rate) and the number of spatial measurement points (i.e., N = 9 for the EC array) respectively.
As proposed by Engelmann and Bernhofer (2016), the combined EC method integrates the temporal and spatial EC methods, yielding fluxes that represent the vertical exchange over the spatial planar covered by the EC array: where . Note that the first term on the right-hand side of Equation (3) represents the space-time averaged kinematic heat flux (i.e., ⟨w ′′ T ′′ ⟩), matching the conventional EC from the temporal perspective.Compared with the single-tower measurements, the second term on the right-hand side of Equation (3) (i.e., ⟨w⟩ ′ ⟨T⟩ ′ ) represents the temporal EC of the spatial mean values over a larger area.Although it signifies the spatial fluctuations induced by large eddies, the validity of Taylor's hypothesis is the prerequisite for applying temporal averaging to this term.

Wavelet spectral analysis
To illustrate the turbulence structures observed by the EC array, wavelet transform is employed to determine the averaged w − T cospectrum.For an arbitrary time-series x(t), the continuous wavelet transform is determined as: where  is the mother wavelet, N and n represent the length of the time-series and the localized time index respectively, and √ t∕s is a normalization factor to ensure unit energy at each scale s.Similar to the traditional fast Fourier transform, the cross-wavelet spectrum of two time-series is calculated as where the superscript asterisk denotes the complex conjugate.Since the cross-wavelet spectrum is complex, the wavelet cospectrum can be expressed as | W xy n (s) |.Therefore, summing up the wavelet cospectrum over all scales equals the covariance between x(t) and y(t).Compared with the traditional fast Fourier transform, the global wavelet cospectrum has the advantage of determining the unbiased and consistent estimation of the true power cospectrum (Torrence and Compo, 1998).

Theoretical scaling argument about the governing factors of spatial variations of w ′ T ′
To explore the dominant effect of large eddies on the spatial variation of w ′ T ′ , and thus the influence on SEB closure, a scaling argument based on the budget equation is performed.For stationary atmospheric surface layer flow with ignorable molecular destruction effect (Katul et al., 2013(Katul et al., , 2014;;Li et al., 2018;Liu et al., 2021), the budget equation of w ′ T ′ is expressed as: where g is the gravitational acceleration.The terms on the right-hand-side of Equation ( 6) represent the production associated with the temperature gradient, the horizontal advection related to the horizontal inhomogeneity of w ′ T ′ , the third-order vertical flux transport, the pressure decorrelation due to the interaction between pressure and temperature, and the buoyancy effect associated with thermal stratification.The pressure decorrelation term can be parameterized by the Rotta (1951) model: where C is a proportionality coefficient and  T is the relaxation time scale, delineating how fast a turbulent eddy loses its coherency.To simplify the calculation, C and  T are considered as constants (Lan et al., 2019).Substituting Equation ( 7) into Equation ( 6) and rearranging to give w ′ T ′ yields Hence, the right-hand-side of Equation ( 8) directly delineates the contribution of gradient diffusion production (i.e., w ′ 2 (T∕z)), horizontal advection (i.e., u ) ), vertical flux transport (i.e., w ′ w ′ T ′ ∕z), and buoyancy destruction (i.e., −(g∕T)T ′ 2 ) to the local kinematic heat flux.
In the constant flux layer (i.e., w ′ T ′ ∕z = 0), the third-order vertical flux transport can be parametrized by the cumulant expansion method, which is based on the joint probability density function of w ′ and T ′ (Raupach, 1981; Nagano and Tagawa, 1995): The influence of large eddies on vertical flux transport is quantified by f , which is a measure of the asymmetry between sweeps and ejections.Furthermore, f is generally parameterized by the incomplete cumulant expansion method (Cava et al., 2006;Katul et al., 2018;Li et al., 2018): where ΔS 0 = S sweep − S ejection represents the imbalance of flux contribution caused by sweeps and ejections (i.e., referred to as the turbulent transport asymmetry), ranging from −1 to +1,  ww =  w ∕u * ,  indicates the ratio of the dimensionless turbulent transport of scalar flux and the dimensionless turbulent transport of the scalar variance, determined by  = M 21 ∕(M 12 − 1), where , where  w and  T are the standard deviations of vertical velocity and temperature respectively.Detailed information about how flux contribution is calculated by quadrant analysis can be found in the Supporting Information and prior studies (Lan et al., 2022;Wang et al., 2022).Taking the partial derivative of Equation ( 9) with respect to x and taking the absolute value for both sides establishes the link between the turbulent transport asymmetric and horizontal inhomogeneity of w ′ T ′ : The terms on the right-hand-side of Equation ( 11) illustrate that the observed horizontal variations in w ′ T ′ are mainly attributed to (a) horizontal inhomogeneity of third-order flux transport (i.e., w ′ T ′ T ′ ∕x), (b) disproportionate changes (e.g., the asymmetry in flux contribution between sweeps and ejections is horizontally inhomogeneous) of turbulent transport asymmetry in the spatial ensemble (i.e., f −1 ∕x), and (c) horizontal inhomogeneity of shear intensity (i.e., u −1 * ∕x).The diurnal variations of air temperature showed a perfect sine shape, and net radiation showed a single peak occurring around noon (Figure 2b,c, Supporting Information Figure S2).Although the diurnal variation of net radiation did show some fluctuations attributed to the cumulus, these observed features indicate generally clear sky conditions as well as the absence of high-impact weather.The minimum temperature occurred at night with a value of 15 • C, whereas the maximum temperature occurred during the day with values up to 36 • C. The R n reached its peak at noon, with average values of 550 W⋅m −2 , and the averaged Bowen ratio during daytime is 0.35, suggesting that the values of LE are equal to 80% of R n .This phenomenon is attributed to the high evapotranspiration in the Pearl River Delta region, consistent with results in the climatological study on the surface energy partitioning in this area (Qian et al., 2017).The average diurnal circle shows that near-neutral and stable conditions occurred either in the early morning or nighttime period, whereas the unstable conditions were mostly observed during the daytime accompanied by peak values of R n (Figure 3).It also captured a salient feature that large SEB non-closure coincided with strong surface heating and wind shear, in compliance with previous studies (Stoy et al., 2013;Gao et al., 2017Gao et al., , 2020;;Liu et al., 2021).

Measurement conditions and intercomparison within EC array
To confirm that each sonic anemometer in the EC array was under the influence of the same flow field, statistical variables (e.g., mean, variance, and covariance) determined by the temporal EC method -that is, Equation ( 1) -from all nine sonic anemometers are intercompared (Table 2).Using EC05 located in the center of the EC array as a reference, average horizontal and vertical wind speeds not only showed high correlation, but their magnitudes were also in good agreement (e.g., S was close to unity and R 2 > 0.95).This finding suggests that all sonic anemometers in the EC array measured the same flow field, further verifying the applicability of the planar-fit method that treats the EC array as a well-defined and static plane (Supporting Information Figures S3 and S4).For the variances of wind velocity components, although the regression slopes were not as close to unity as the  mean wind velocity, good correlations were still observed (Supporting Information Figures S5 and S6).Compared with the flow field, the temperature field showed relatively large variations, as evidenced by the comparatively small correlations (Supporting Information Figures S7 and S8).This finding raises the point that the temperature might vary slightly across patches even in a small area with continuous grass cover.It is worth noting that since  v measured by the EC main tower was employed to perform the Schotanus-Nieuwstadt-Debruin correction throughout all sonic anemometers, the inevitable spatial difference in the air humidity would also contribute to the large variances in the temperature field.As a result, the kinematic heat flux varied substantially within the EC array, challenging the interpretation of flux measurement from a single EC tower (Supporting Information Figure S9).It further motivates the investigation of the contribution of such spatial flux variation to the SEB non-closure.were always larger than w ′ T ′ .Moreover, relatively large discrepancies generally coexisted with peak values of kinematic heat fluxes, implying that the spatial flux contribution was more prominent around noon (Figure 4b).This finding supports previous studies that spatial covariance significantly contributes to vertical flux under strong radiative forcing conditions (e.g., large value of R n and unstable stratification) (Mauder et al., 2008;Stoy et al., 2013;Engelmann and Bernhofer, 2016).Since the first term on the right-hand-side of Equation ( 3) represents the 30-min bin-average value of ⟨w ′′ T ′′ ⟩, it is reasonable to compare its contribution to the vertical flux in the temporal ensemble under unstable conditions.Although there is no systematic difference between ⟨w ′′ T ′′ ⟩ and w ′ T ′ , ⟨w ′′ T ′′ ⟩ only accounts for 89% of the vertical fluxes measured by a single tower.In comparison, aggregating the spatial flux contribution associated with large eddies bring forth the increment in both regression slope and correlation (e.g., S and R 2 increase from 0.89 to 1.04 and from 0.85 to 0.93 respectively).On average, employing the combined EC method considerably increases the vertical kinematic heat flux by 21% compared with the traditional temporal EC method.

Combined EC method yields improved SEB closure
To further demonstrate how the SEB closure can be improved by employing the combined EC method under unstable conditions, Figure 5a shows the comparison between the sum of turbulent heat fluxes (i.e., H + LE) and  available energy (R n − G 0 ).Note that H combined is derived from all nine sonic anemometers' measurements using the combined EC method -that is, kinematic heat flux is determined using Equation ( 3) and then converted to the sensible flux -whereas others (i.e., H temporal , LE, R n , and G 0 ) are derived from the EC main tower measurements using the traditional temporal EC method.As expected, the combined EC method considerably improves the SEB closure, supported by the enhanced correlation and regression slope that approach unity (Figure 5a).Another noticeable feature shown in Figure 5a is that data points become more scattered when with enlarged available energy.Furthermore, the magnitude of H combined not only increases but also varies substantially in coincidence with enhancing u * (Figure 5b).Although such phenomena have been brought to light by prior studies (e.g., Balogun et al., 2009;Gao et al., 2017Gao et al., , 2020;;Kutikoff et al., 2019), the role of large eddies in affecting the SEB non-closure is not fully understood and will be further investigated by a scaling argument provided in the following sections.

Horizontal variation of w ′ T ′ attributed to large eddies
Figure 6 shows the averaged normalized w − T wavelet cospectra as a function of non-dimensional frequency (n = fz∕⟨U⟩, where f is natural frequency, z is measurement height, and ⟨U⟩ is the mean wind speed determined by temporal-spatial averaging).It is worth noting that the individual cospectrum for each 30-min data segment is normalized by the variances of vertical velocity and temperature before averaging.It can be seen that the non-dimensional frequencies corresponding to the   cospectral peak are comparable within the EC array.Moreover, regardless of the slight difference in cospectral energy, both the non-dimensional frequency corresponding to the cospectral peak and the energy cascade pattern in the inertial subrange can well depicted by the model cospectrum (Kaimal et al., 1972).It further corroborates the validity of employing the planar-fit method for tilt correction, providing confidence that each individual sonic anemometer in the EC array measured the same flow field, which is the prerequisite for the following theoretical argument.Nevertheless, the scattered data points shown in Supporting Information Figure S9 and the improved SEB closure resulting from the combined EC method (Figure 5) highlight the underestimation of turbulent fluxes derived from the traditional EC method, implying the crucial role of large eddies in modulating the spatial flux, and thus affecting the SEB closure (Kanda et al., 2004;Foken, 2008;Mauder et al., 2008).
To further delineate the role of large eddies in influencing the horizontal variation of w ′ T ′ , the scale dependence of difference in normalized kinematic heat flux is evaluated.Since the global wavelet cospectra represent the energy density as a function of non-dimensional frequency, the area under any portion of the cospectra curve continues to be proportional to the covariance.Therefore, the scale-dependent cospectral difference between each individual sonic anemometer and the reference is estimated by It is worth noting that the influence of shear intensity is emphasized by dividing Δf CO into five groups with different ranges of u * (Figure 7).Distinct features with respect to different turbulent scales are observed in Figure 7, interpreting explicit contributions of large eddies  and small eddies to the horizontal inhomogeneity of w ′ T ′ .For Δf CO in the high-frequency range (i.e., n > 10 −2 ), the cospectral difference drops dramatically as the scale of turbulent eddies reduces and the Δf CO is independent of u * , indicating that the background small turbulent eddies are isotropic, rarely contributing to the horizontal variation of w ′ T ′ .Despite the sudden increase of cospectral difference in the low-frequency end, Δf CO increases then drop with the enlarged turbulent scales, showing an energy peak at the mid-frequency range (5 × 10 −3 < n < 10 −2 ).In this frequency range, the conspicuous feature is that Δf CO for the larger u * groups have larger magnitudes than the counterparts with smaller u * .Plainly, such large eddies with frequencies ranging from5 × 10 −3 to 10 −2 are mainly responsible for the horizontal inhomogeneity of w ′ T ′ , thereby explaining the observed feature that large SEB non-closure concurred with increased u * .Nevertheless, this observed correlation does not stand for the direct mechanism of how large eddies regulate the SEB closure.To further explore such causation, diagnosing secondary circulations and quantifying the associated flux transports should be conducted using numerical simulations (e.g., large eddies simulation and land surface model) over surfaces with a broad range of heterogeneity.

Inhomogeneous asymmetric turbulent flux transport
Based on Equation ( 8), the influence of vertical flux transport on local turbulent flux has been investigated using multi-level EC measurements (Li et al., 2018;Lan et al., 2019;Liu et al., 2021).The results found that the disproportionate changes in turbulent transport asymmetry with instability regulate the varying flux contributions and coincide with the widely observed increase in the non-closure with increasing instability.However, the link between horizontal inhomogeneity of w ′ T ′ attributed to large eddies and the variation of SEB closure with u * is rarely explored and will be discussed in this section.increases with | w ′ T ′ ∕x |, such a horizontal difference in kinematic heat flux cannot be captured by a single-tower measurement.On the one hand, it explains the observed feature in Figure 3, that large non-closure with strong wind shear.On the other hand, it confirms that the combined EC method, which integrates the temporal and spatial fluxes, yields larger magnitudes in the kinematic heat flux and thus improves the SEB closure (e.g., Figures 4 and 5).To further determine the dominant factor that contributes to the horizontal inhomogeneity of w ′ T ′ , the variation of | w ′ T ′ ∕x | is presented as a function of each individual term on the right-hand-side of Equation ( 11).Both calculations (Table 3) and Figure 10 confirm that the horizontal inhomogeneity of w ′ T ′ is largely attributed to the disproportionate changes of turbulent transport asymmetry in the spatial ensemble.For instance, the average value of the second term on the right-hand-side of Equation ( 11) has the largest value, compared with the other terms:

Sonic anemometer
01.On the one hand, as evidenced by the increased flux fraction related to both sweeps and ejections, flux contribution related to large eddies increases with enhancing wind shear.On the other hand, the prominent increment of flux contribution is observed in sweeps compared with the ejections (Figure 11).This further supports that the enlarged turbulent transport asymmetry is primarily attributed to the disproportionate increment of flux contribution associated with sweeps and ejections.Hence, the frequently observed dispersive fluxes are explained by such enlarged turbulent transport asymmetry associated with large eddies under strong wind shear, which cannot be captured by the single EC tower and the traditional temporal EC method, and thus degrade the SEB closure.

CONCLUSIONS
Using a high-density spatial EC set-up over an area with continuous grass cover, the SEB is evaluated using both temporal and combined EC methods.The calculations indicate that ⟨w ′′ T ′′ ⟩ accounted for ∼89% of the kinematic heat flux measured by a single tower, whereas the kinematic heat flux determined by the combined EC method increased to 104%, indicating that the vertical flux cannot be completely captured by a single tower.The EC array was able, to some degree, to resolve such a spatial flux contribution related to large eddies, thereby increasing the SEB closure.Moreover, it is observed that, under unstable conditions, the enlarged energy imbalance coincided with increasing u * consistent with many other SEB experiments.The relation between surface energy imbalance and horizontal variation of flux, as well as the spectra analysis, suggests the predominant role of large eddies in contributing to the dispersive flux, explaining the concurrence of enlarged surface energy imbalance and increasing u * .Based on the turbulent flux budget equation, the scaling analysis confirms that changes in the horizontal variation of flux are primarily explained by the disproportionate changes of turbulent transport asymmetry in the spatial ensemble.As u * increases, the flux contribution related to sweeps exhibits greater increment compared with the counterpart associated with ejections.Consequently, disproportionate changes in such a turbulent transport asymmetry with u * regulate the horizontal variation of flux and coincide with the widely observed enlarged SEB non-closure with increasing wind shear.

F
I G U R E 1 Site information of the Guangdong Provincial Observation and Research Station: (a) the Google Earth image of the observation site -adapted from Lan et al. (2022); (b) a photograph of the spatial eddy covariance (EC) array set-up taken by the drone; (c) the schematic illustration of the spatial EC array set-up; (d) wind rose diagram during the intensive observation period from September 5 to October 31, 2018.[Colour figure can be viewed at wileyonlinelibrary.com]

Figure 2
Figure2portrays the measurement conditions during the entire intensive observation period.Wind speed showed relatively high fluctuations (Figure2a), and the wind direction was mostly southwest (Figure1d), ensuring the consistency of wind pattern during the entire intensive observation period and minimizing the flow perturbation due to the trees to the northeast of the EC array.Regardless of a few data gaps induced by power outages and poor data quality, the measurement conditions were fairly excellent, Average diurnal variation of the (a) friction velocity and non-dimensional stability parameter ( = z∕L, where z and L are the observation height and Obukhov length respectively); (b) energy fluxes and the residuals of the surface energy balance closure (Imb).All variables shown here are determined by the traditional temporal eddy covariance (EC) method using EC main-tower measurements.LST: local sidereal time; UTC+8.[Colour figure can be viewed at wileyonlinelibrary.com]

Figure
Figure 4a provides a visual illustration of the diurnal courses of kinematic heat flux derived by the temporal (i.e., w ′ T ′ determined by measurements from EC05) and combined EC (i.e., [ w ′′′ T ′′′ ] ) methods.On the one hand, diurnal variations of w ′ T ′ and [ w ′′′ T ′′′ ] exhibited high correlation, characterized by the single peak pattern a, b) The diurnal course of (a) the kinematic heat flux derived from the temporal (blue line) and the combined eddy covariance (EC) method (red line) and (b) the difference of kinematic heat fluxes derived from the combined and temporal EC method.(c, d) The comparison of the kinematic heat flux determined by (c) the space-time averaged method -first term on the right-hand-side of Equation (3) -and (d) the combined and temporal EC method.[Colour figure can be viewed at wileyonlinelibrary.com] a) Comparison of turbulent fluxes (H + LE) and the available energy (R n − G 0 ) for the 30-min segments that meet the data selection criteria; (b) variations of H derived from the temporal and combined eddy covariance (EC) method as a function of friction velocity.Light blue and magenta markers respectively indicate the H derived from the temporal and combined EC method.The dashed black line denotes the 1-1 line.[Colour figure can be viewed at wileyonlinelibrary.com]Averaged normalized w − T cospectra of the 30-min segments that meet the data selection criteria.The shaded area and black dashed line respectively represent the standard deviation of the averaged cospectra and the model cospectra based on the Kansas experiment (Kaimal et al., 1972).[Colour figure can be viewed at wileyonlinelibrary.com] Scale dependence of the cospectral difference between each individual sonic anemometer and the reference (EC05) for different ranges of u * as listed in the legend in (a).The shaded area represents the standard deviation of the averaged scale-dependent cospectral difference for the 30-min segments that meet the data selection criteria.[Colour figure can be viewed at wileyonlinelibrary.com] ) F I G U R E 8 Variations of the residuals of the surface energy balance closure (Imb) with the horizontal difference of kinematic heat flux between each individual sonic anemometer and the reference (EC05).The continuous lines are markers in subplots are the curves determined by the unweighted bin-averaged method.[Colour figure can be viewed at wileyonlinelibrary.com]

Figure 8
Figure 8 demonstrates the variation of the surface energy imbalance (i.e., Imb = R n − G 0 − H − LE) as a function of | w ′ T ′ ∕x |.It can be seen that the surface energy imbalance enlarges with increasing | w ′ T ′ ∕x |.The dependence of | w ′ T ′ ∕x | on u * confirms that strong fluctuations in spatial covariance always coincide with Figure 9 shows that | w ′ T ′ ∕x | increases with enhancing u * , indicating that horizontal difference in kinematic heat flux becomes more prominent when wind shear enhances.Although Figure 8 implies that the magnitude of surface energy imbalance F I G U R E 9 Variations of the horizontal difference of kinematic heat flux between each individual sonic anemometer and the reference (EC05) with u * .The continuous lines in subplots are the curves determined by the unweighted bin-averaged method.[Colour figure can be viewed at wileyonlinelibrary.com]TA B L E 3 Ensemble of each individual term in Equation (10) in units of K⋅s −1 .

F
I G U R E 10 Relations between changes in the horizontal difference of kinematic heat flux between each individual sonic anemometer and the reference (EC05) and the contributions of disproportionate changes of turbulent transport asymmetry in the spatial ensemble (| u −1 * w ′ T ′ T ′ ( f −1 ∕x ) |).The continuous lines in subplots are the curves determined by the unweighted bin-averaged method.[Colour figure can be viewed at wileyonlinelibrary.com]

1477870x, 0 ,
Downloaded from https://rmets.onlinelibrary.wiley.com/doi/10.1002/qj.4562by Karlsruher Institution F. Technologie, Wiley Online Library on [12/09/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Variations of flux contributions (left y-axis) from ejections (red markers) and sweeps (blue markers) and the turbulent transport asymmetric (right y-axis) associated with large eddies with u * for kinematic heat flux.Red and blue continuous lines indicate the bin-averaged curves related to flux contributions from ejection and sweeps respectively.Black curves are the bin-averaged curves of the turbulent transport asymmetry.[Colour figure can be viewed at wileyonlinelibrary.com]

TA B L E 1
Instrumentation of the continuously operating micrometeorological tower (eddy covariance main tower).
Intercomparison of statistic variables measured by the eddy covariance array.
TA B L E 2Note: S and R 2 respectively represent the slope of linear regression and the coefficient of determination related to the comparison of the listed variables between each individual sonic anemometer and the reference.