Non-closure of the surface energy balance explained by phase difference between vertical velocity and scalars of large atmospheric eddies

It is now accepted that large-scale turbulent eddies impact the widely reported non-closure of the surface energy balance when latent and sensible heat fluxes are measured using the eddy covariance method in the atmospheric surface layer (ASL). However, a mechanistic link between large eddies and non-closure of the surface energy balance remains a subject of inquiry. Here, measured 10 Hz time series of vertical velocity, air temperature, and water vapor density collected in the ASL are analyzed for conditions where entrainment and/or horizontal advection separately predominate. The series are decomposed into small- and large- eddies based on a frequency cutoff and their contributions to turbulent fluxes are analyzed. Phase difference between vertical velocity and water vapor density associated with large eddies reduces latent heat fluxes, especially in conditions where advection prevails. Enlarged phase difference of large eddies linked to entrainment or advection occurrence leads to increased residuals of the surface energy balance.


Introduction
The non-closure of the surface energy balance, i.e. the failure of the sum of the turbulent fluxes of sensible and latent heat (H þ LE) to account for the available energy (R n À G 0 , where R n and G 0 are the net radiation and the ground heat flux, respectively), remains an unsolved problem in micrometeorology , Foken 2008, Leuning et al 2012. The lack of surface energy balance closure raises concerns when utilizing eddy covariance flux data to evaluate or calibrate land surface flux schemes in weather and climate models (Williams et al 2009, Blyth et al 2010. A number of causes have been proposed to explain the surface energy balance nonclosure, including mismatch in instrument footprint, measurement errors in all components of the surface energy balance, advective flux divergence, and inadequate sampling of large-scale, low-frequency turbulent motions (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, Russell et al 2015. Any violations in the assumptions (e.g. steady state and horizontally homogeneous landscapes) and corrections (e.g. density effects) made for eddy convariance measurements also contribute to the non-closure. However, even taking into account many of these causes adequately in specially designed field experiments or accounting for all terms rigorously, a surface energy imbalance was still reported (Mauder et al 2007, Foken 2008, Leuning et al 2012. A recent panel discussion concluded that large eddies are likely to be one of the primary contributors to the underestimated turbulent fluxes and thus the non-closure of the surface energy balance . This view is supported by Large Eddy Simulation (LES) studies (Kanda et al 2004) that demonstrate lack of closure is attributed to finite time averaging of large-scale convective structures. However, mechanistic understanding of how large eddies affect turbulence structures and what factors determine the contribution of large eddies to fluxes remains a subject of inquiry.
Variations in large eddies can be spectrally characterized by phase difference between large structures for vertical wind speed (w, m s À1 ) and scalars such as temperature (T, K) or water vapor density (q, g m À3 ) along with their magnitudes (Mahrt 1991, Li andBou-Zeid 2011). As shown here, the phase difference ∅ between large-scale w and q is far more critical in explaining the lack of energy balance closure than their spectral amplitude counterparts. Consider two low frequency, sinusoidal time series x t ð Þ ¼ sin vt ð Þ and y t ð Þ ¼ sin vt þ ∅ ð Þwith a certain frequency v but the same unitary amplitudes separated by ∅. This phase difference between x t ð Þ and y t ð Þ also leads to a time-lag t ¼ ∅=v between the two series. Hence, low frequency (or small v) produces large lags for finite ∅ and conversely, small lags are associated with high frequencies for the same finite ∅. It can be readily shown by time elimination that y ¼ xcos ∅ ð Þ ∓ ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À x 2 p sin ∅ ð Þ: For ∅ ¼ 0; the two series are perfectly correlated (y ¼ xÞ with maximum flux contribution as indicated by their covariance irrespective of whether v is large or small. As ∅ increases, the correlation between these two series is weakened and their covariance is reduced, which is caused by two mechanisms: a direct slope reduction due to finite cos ∅ ð Þ and a hysteretic term generating loop-like behavior ( ∓ ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À x 2 p Þ modulated by sin ∅ ð Þ. Hence, any shifts in phase between the large structures of w and those of q cause reductions in their covariance magnitude and contribute to the energy balance non-closure. This conjecture and potential causes for the phase difference at small v are explored using high frequency atmospheric surface layer (ASL) data collected during the Energy Balance Experiment (EBEX-2000) described elsewhere (Gao et al 2016), where advection predictably impacts the energy balance non-closure during certain conditions.

Experimental data
The data used were measured during EBEX over an irrigated cotton field in the San Joaquin Valley, California, USA (Pacific Daylight Time, PDT = UTCÀ7 h; 36°06 0 N, 119°56 0 W; 67 m a.s.l.) from 20 July to 24 August 2000 (Mauder et al 2007. Although the topography of the experimental field was flat with canopy height of about 1 m, patch-bypatch irrigation created heterogeneous soil moisture conditions across the field . Ten tower sites were erected during the field campaign. The eddy covariance data analyzed here were collected at Site 7 by one three-dimensional sonic anemometer (CSAT3, Campbell Scientific, Inc.) and one krypton hygrometer (KH20, Campbell Scientific, Inc.) measuring longitudinal, lateral and vertical wind speed components (u, v and w), T, and q at 8.7 m. The data were collected at 10 Hz using a datalogger (CR23X, Campbell Scientific, Inc.) and then stored for future post-processing. All sensors were compared before the experiment (Mauder et al 2007) and were corrected accordingly. Routine maintenance was performed to keep the source and detector windows of KH20 to minimize the scaling, and the calibration was also performed during the experiment (Foken and Falke 2012). Other ancillary data included 30 min averaged net radiation (CNR1, Kipp & Zonen, Inc.), air pressure, air temperature and relative humidity (HMP45C, Vaisala, Inc.), soil temperature, soil water content (Q soil ; CS615, Campbell Scientific, Inc.), and soil heat flux (G(z p ); HFT-3.1, Campbell Scientific, Inc.; z p = 80 mm below the soil surface). Only the daytime data (i.e. 08:00À18:00 PDT) were used to ensure that measured R n À G 0 far exceeds instrument noise and some uncertainties associated with the determination of G 0 explained later .

Post-field data processing
The post-field data processing procedures are documented elsewhere (Zhang et al 2010, Zhang and Liu 2014, Gao et al 2016. Briefly, these procedures include: 1) removing physically impossible values and spikes from the time series; 2) double rotation for the sonic anemometer measured w (Kaimal and Finnigan 1994); 3) calculation of averages, variances and covariances using a 30 min block average; 4) sonic temperature correction (Schotanus et al 1983, Liu et al 2001, oxygen cross sensitivity correction for KH20 (Tanner et al 1993), density correction applied to LE (Webb et al 1980), and correction for flux attenuation due to spatial separation of CSAT3 and KH20 ; and 5) quality check for stationary and developed turbulent conditions (Foken et al 2004).

Ground heat flux (G 0 )
Using the method in Liebethal et al (2005) and Oncley et al (2007), G 0 was determined from measurements by heat flux plates at a depth z p and the heat stored in the soil layer between the surface and z p , where T soil is the mean soil temperature between the surface and z p , and c soil is the volumetric heat capacity of moist soil calculated by: where rc ð Þ water and ðrcÞ soil;dry are the volumetric heat capacities of water (4.2 Â 10 6 J m À3 K À1 ) and soil minerals (2.5 Â 10 6 J m À3 K À1 ), respectively, and Q soil is the layer-averaged volumetric water content. Unlike soil heat conductivity that can vary with the particular configuration of soil pore structure, soil particle contact areas and water bridge configurations between them, specific heat capacities vary linearly in proportion to the constituent material. Previous studies suggested that the method described above is one of the highly recommended methods depending on data availability (Liebethal et al 2005, Russell et al 2015.

Ensemble empirical mode decomposition (EEMD)
To quantify the phase difference between lowfrequency series of w, q, and T caused by large eddies, we first detected and extracted large eddies (i.e. low-frequency structures) from the 10 Hz time series data. Traditionally, Fourier-based filters and wavelet decompositions were deployed to separate large eddies from an oscillatory signal using pre-determined basis functions. Non-linear and non-stationary nature of turbulence could prevent successful implementation of the traditional methods on turbulence decomposition because both Fourier and wavelet transforms cannot readily accommodate non-linear and nonstationary series without 'extra' ad-hoc post-processing (Huang et al 1998, Huang andWu 2008). Compared to Fourier-based filters and wavelet decompositions, the empirical mode decomposition (EMD) is adaptive and has high locality, and therefore capable of handling some of the non-linear and non-stationary occurrences that are ubiquitous in turbulence data (Huang et al 1998, Huang andWu 2008). EMD has been favored in several atmospheric boundary layer (ABL) studies (Hong et al 2010, Barnhart et al 2012, Wang et al 2013, Gao et al 2016. The ensemble empirical mode decomposition (EEMD) (Wu and Huang 2009), a method based on EMD (Huang et al 1998, Huang andWu 2008) for time series analysis, was adopted here. Using a sifting process, EEMD decomposes a time series x t ð Þ into a finite set of amplitude-frequency modulated oscillatory components C j,x (j = 1, 2, . . . , n) and a residual r n (t). This residual is a time series that is either monotonic or containing only one extremum from which no more oscillatory components can be further extracted. Hence, To ensure the equivalence of variances and covariances calculated from all extracted components (i.e. C j,x and r n ) and the post-field data processing, the EEMD was applied to each 30 min time series data of w 0 , q 0 and T 0 . The prime denotes the turbulent fluctuations relative to the 30 min block average after despiking and double rotation (i.e. steps 1 and 2 in section 3.1). EEMD can only sift out oscillatory components that differ in period by more than a factor of two (Huang et al 1998, Huang and Wu 2008, Barnhart et al 2012. For each 30 min time series run, thirteen oscillatory components (i.e. n ¼ 13) were extracted with one overall residual r 13 , marked as the fourteenth component (C 14;x ). The residual was included because it still contributes to the total variance and possibly covariance as compared to the traditional post-field data processing (Barnhart et al 2012). Each oscillatory component has its own characteristic frequency (Wang et al 2013). Since all time-domain components are additive, the sum of certain oscillatory components can be interpreted as large eddies if a critical frequency is identified and assigned as a delineator of large scales (Wang et al 2013).
Several methods for identifying a critical frequency have been proposed to separate large eddies from small eddies, including a spectral gap in the v spectra (Cava and Katul 2008), a copsectral gap in uw and wT cospectra (Vickers and Mahrt 2003), and the first clearly identified maximum in wavelet variance spectrum (Thomas and Foken 2005, Barthlott et al 2007, Zhang et al 2011. In the EBEX-2000 dataset, a spectral gap at about 0.02 Hz did exist for some cases in their v spectra (Zhang et al 2010, Gao et al 2016. Zhang et al (2011) calculated the characteristic scales of coherent structures (i.e. large eddies) under different atmospheric stability conditions using the same dataset, and found that the averaged frequency for dominant large eddies was also about 0.02 Hz during the daytime for this site. Based on these previous studies, oscillatory components with the mean frequencies smaller than 0.02 Hz are labeled as large eddies (i.e. x 0 l ¼ P 14 j¼9 C j;x , where x = w, q, and t, and l refers to large eddies for x), and larger than 0.02 Hz as small eddies (i.e. x 0 s ¼ P 8 j¼1 C j;x , where s refers to small eddies for x). Different tests indicated that the selection of the critical frequency only has minor effects on quantifying phase difference and flux contribution. Turbulent fluxes of large and small eddies were corrected following step 4 in the post-field data processing (section 3.1).
There are three methods for estimating phase difference between different quantities: Fourier transform, wavelet transform, and Hilbert transform (Stull 1988, Quiroga et al 2002, Grinsted et al 2004. Separate tests (not shown here) indicated that these three methods resulted in similar phase difference when sinusoidal signals with known phase difference were processed. Nevertheless, taking into account the concerns in previous studies that both Fourier and wavelet transforms may induce uncertainty when nonlinear and non-stationary data are processed, the Hilbert transform was deemed better suited for this purpose (Huang et al 1998, Quiroga et al 2002.

Hilbert transform (HT)
For a signal xðtÞ, its Hilbert transform,x t ð Þ, is given as,x Environ. Res. Lett. 12 (2017) 034025 where PV denotes a Cauchy principle value integral. Mathematically, the HT is a phase shift of the original signal by À90°. Thus, an analytic signal can be constructed from the signal and its Hilbert transform as , A is the instantaneous amplitude, and F is the instantaneous phase function. Because the Hilbert transform is by construct sensitive to instantaneous phase shifts encoded as t À t in equation (4), it offers advantages over wavelet or Fourier methods.
Similar to Fourier and wavelet cross spectrum, the Hilbert cross spectrum of two signals xðtÞ and yðtÞ is defined as where Ã denotes complex conjugation. The instantaneous Hilbert phase difference between two signals is obtained as where Im and Re denote the imaginary and real parts of the Hilbert cross spectrum, respectively. Here, the time-averaged F dif f t ð Þ associated with large eddies is used to represent the phase difference between w 0 l and q 0 l for each run.  1(a)). The variations suggested that other external forcing beyond R n affected the evapotranspiration, because R n followed a nearly ideal sine shape, with a maximum value of about 700 W m À2 at 13:00 PDT (local noon; figure 1(a)). When partitioning LE into flux contributions by large eddies (i.e. LE l caused by w 0 l and q 0 l ) and small eddies (i.e. LE s caused by w 0 s and q 0 s ), more pronounced point-to-point variations in LE l was noted compared to LE s during the late morning to early afternoon. In the late afternoon after 14:00 PDT, LE l showed a large decreasing trend and approached Environ. Res. Lett. 12 (2017) 034025 zero by about 18:00 PDT, whereas LE s exhibited a small decreasing trend during the entire afternoon ( figure 1(b)). The decreased LE l was largely attributed to horizontal advection under warm (dry)-to-cool (wet) transition from the upstream drier patches to the wet cotton field , Zhang et al 2010, leading to a reduced surface energy balance closure in the afternoon than in the morning at the study site. As shown in table 1, a case study is featured illustrating that despite the comparable magnitudes between low-frequency signals of w and those of q caused by large eddies (i.e. w 0 l and q 0 l ), the phase difference between them led to large reduction in the contribution of large eddies to LE. Would this conclusion be generally applicable to explain the influence of large eddies on ASL turbulence and energy balance non-closure? How do simultaneous changes in both magnitudes and phase difference between two quantities (i.e. w 0 l and q 0 l ) modulate the influence of large eddies on turbulent fluxes? To make our calculation manageable, the daytime data in EBEX were divided into five groups according to the ranges of the magnitudes of w 0 l (A w l ) and q 0 l (A q l ). Each group consists of 30 min data points with A w l A q l being located within a magnitude range of 0-0.1, 0.1-0.15, 0.15-0.2, 0.2-0.25, and > 0.25 g m À2 s À1 (= m s À1 g m À3 ). Figure 2(a) shows that for all five groups, LE l generally decreased with the increased phase difference between w 0 l and q 0 l , though the decrease trends of LE l with the increasing phase difference were different for the groups. The larger the magnitude ranges, the greater the slopes, suggesting that the flux contribution from more energetic large eddies was more sensitive to changes in phase shifts. Given the same phase difference, LE l was greater for the larger magnitude ranges, reflecting that the w 0 l and q 0 l magnitude affected flux contributions of large eddies. Figure 2(a) also indicates that LE l varied across larger ranges under smaller phase difference (e.g. jF dif f j ¼ 40°) than larger phase difference (e.g. jF dif f j ¼ 80°). LE l from all groups converged to smaller values when phase differences became larger, and all groups approached zero when jF diff j became 90°. When w 0 l and q 0 l became out of phase, the large eddies contributed little to LE no matter how energetic large eddies are (as quantified by their w 0 l and q 0 l magnitudes). For the study site, H was relatively small and even negative during the daytime, while it is still true that H contributed by large eddies (i.e. w 0 l and T 0 l ) decreased with the increasing phase difference given the comparable magnitudes of A w l A T l (figure S1 available at stacks.iop.org/ERL/12/034025/mmedia). We also computed the flux contribution and phase difference between small eddies of w and q (i.e. w 0 s and q 0 s ). The mean absolute phase difference between w 0 s and q 0 s only varied within a small range of 60°À78°b ecause small eddies are becoming locally isotropic and de-correlated because of the action of pressurevelocity and pressure-scalar interactions (Li and Bou-Zeid 2011). As a result, the decrease in LE s with the increase in phase difference between w 0 s and q 0 s was not as obvious as that for the large eddies (figure S2).

Phase difference largely responsible for the nonclosure of surface energy balance
Since enlarged phase difference between w 0 l and q 0 l led to the reduction in LE l , it is expected that the surface energy balance closure was degraded with enlarged phase difference. This statement is confirmed by figure  2(b) in which the ratio of the sum of the turbulent fluxes to the available energy (i.e. (H þ LE)/(R n À G 0 )) decreased with the increase in phase difference between w 0 l and q 0 l . As stated above, LE l increased with A w l A q l at fixed jF dif f j. Therefore, (H þ LE)/(R n À G 0 ) also increased with A w l A q l given the same phase difference (figure 2(c)). To demonstrate whether the reduced phase difference alone is able to improve the closure for the dataset in EBEX-2000 presented in figure 2(d), we artificially modified the phase difference between w 0 l and q 0 l (T 0 l ) with their amplitudes intact by using Fourier transform, as the inverse Fourier transform of a signal equals to the original signal. For instance, the phase of q 0 l was artificially set to the phase of w 0 l but amplitudes were not altered. This artificial adjustment would cause an increase of about 27% in LE (figure S3, table S1). For w 0 l and T 0 l , under unstable conditions, the phase of T 0 l was artificially set to the phase of w 0 l but, under stable conditions, the phase of T 0 l was artificially set to be out of phase with w 0 l so that w 0 l and T 0 l would have a negative flux contribution. The linear regression between the original H and the modified H has a slope of 1.25 ± 0.02 and an intercept of (2.02 ± 0.71) W m À2 with an absolute value for R 2 of 0.94 (figure S3, Table 1. Characteristics of large eddies for the three cases A (11:00-11:30 PDT), B (11:30À12:00 PDT), and C (12:00-12:30 PDT), respectively. LE l is the latent heat flux contributed by large eddies. A w l ± SD and A q l ± SD are mean and standard deviation of the instantaneous amplitudes for w 0 l and q 0 l , respectively. jF dif f j ± SD is mean and standard deviation of the absolute instantaneous phase difference between w 0 l and q 0 l in units of degree. DLE l is the increase of the latent heat flux contributed by large eddies by artificially modifying the phase difference between w 0 l and q 0 l with their amplitudes intact by using Fourier transform and its inverse transform. table S1). After the modification of phase difference for all data points, the surface energy balance closure improved significantly (P < 0.001) from 0.77 ± 0.18 to 0.97 ± 0.20 with similar values of R 2 as shown in table 2. The above tests suggest that phase differences between w 0 l and q 0 l were responsible for a significant portion of the energy balance non-closure.

Variations of phase difference
What causes these phase differences in large eddies is difficult to discern from single-point time series analysis. However, a number of conjectures may be offered. The daytime unstable ASL over moist landscapes is statistically predominated by upward warm/moist (e.g. w 0 s > 0 and q 0 s > 0) and downward cool/dry (e.g. w 0 s < 0 and q 0 s < 0) turbulent eddies when their size is commensurate with the measurement height z m , leading to positive H and LE The 'post-phase shift' means that the phase of large eddies for scalars (temperature and water vapor density) was forced to equal to the phase of large eddies for vertical velocity. This is the only difference between the calculation of 'post-phase shift' and 'pre-phase shift' turbulence fluxes. The dash lines are the linear best fits, and their linear regression coefficients are listed in table 2. The gray line is a 1:1 line. Table 2. Linear regression coefficients for energy balance closure for 'pre-phase shift' and 'post-phase shift' . EBR refers to mean and standard deviation of the ratio of (Hþ LE)/(R n ÀG 0 ). The difference between EBR for 'pre-phase shift' and 'post-phase shift' is statistically significant (P < 0.001).

Intercept
Slope R 2 EBR upward warm/moist (e.g. w 0 l > 0 and q 0 l > 0) and downward cool/dry (e.g. w 0 l < 0 and q 0 l < 0) large eddies would contribute positively to sensible and latent heat fluxes (e.g. þw 0 l À Á þq 0 l À Á > 0 and Àw 0 l À Á Àq 0 l À Á > 0). Under this circumstance, largescale changes in time-series of vertical wind velocity and scalars (e.g. water vapor density) are considered to be in phase (i.e. phase difference jF dif f j ¼ 0°). If downward moving large eddies carry warm/moist air masses (e.g. w 0 l < 0 and q 0 l > 0), these large eddies would contribute negatively to sensible and latent heat fluxes (e.g. Àw 0 l À Á þq 0 l À Á < 0). In this case, large-scale changes in the time-series of vertical wind velocity and scalars (e.g. water vapor density) for large eddies are considered to be perfectly out-of-phase (i.e. jF dif f j ¼ 180°). In reality, large-scale phase between vertical wind velocity and scalars varies from 0°to 180°. Hence, in the absence of advection and for the case of entrainment and storage (leading to a time lag) within the much studied convective ABL, it follows that F diff will be finite.
The problem of horizontal advection of dry air onto a wet surface is of particular interest here because its qualitative effects on the energy balance nonclosure is predictable and serves as another test case for understanding the genesis of F diff . In the presence of advection, the mean continuity equation for water vapor density at z m above a canopy of height h reduces to where U is the mean longitudinal wind speed and x and z are the longitudinal and vertical directions, respectively. The overbar denotes the 30 min average. Because advection of dry air onto a wet surface generates a positive ∂q=∂x, the U ð∂q=∂xÞ is always positive (as is the case here) and w 0 q 0 z m ð Þ ¼ w 0 q 0 h ð Þ À Ð z m h U @q @x dz: Because jw 0 q 0 z m ð Þj < jw 0 q 0 h ð Þj, the proper energy balance residual ¼ R n À G o À L v w 0 q 0 h ð Þ is expected to be smaller than the measured R n À G o À L v w 0 q 0 z m ð Þ, where L v is the latent heat of vaporization of water. The jF diff j for these advective conditions is shown in figure 3. The stability parameter z is calculated using z ¼(z m À d)/L, where d is the zero-plane displacement and L is the Obukhov length. The analysis in figure 3 demonstrates that advective conditions did result in negative sensible heat flux (i.e. z > 0 for daytime values), increased energy balance residual, and a correlation coefficient between air temperature and water vapor density that is sufficiently negative consistent with expectations of dry and warm air parcels advecting onto the sensor. In these cases, the measured q 0 2 appears larger than that predicted from surface fluxes represented by q Ã ¼ Àw where q Ã and u Ã are the scale for water vapor density and the friction velocity, respectively, and those large excursions are associated with increases in jF dif f j. Hence, it can be surmised that advection increases the phase between w 0 l and q 0 l presumably due to finite travel time of dry and warm air parcels to the sensor location. Both entrainment/storage or advective travel time lead to increases in jF dif f j with afternoon advection having a more pronounced impact on the energy balance non-closure (i.e. larger jF dif f j).

Conclusions
Using data measured over a flood-irrigated cotton field where the latent heat flux is roughly commensurate with net radiation, the effects of large eddies on fluxes and the surface energy balance closure are analyzed. Large eddies caused 30 min mean, point-topoint variations in daytime latent heat fluxes. The case study and the statistical analyses indicated that phase difference between large-scale vertical wind speed and water vapor density caused by large eddies were the primary cause for point-to-point variations in latent heat fluxes and thus the corresponding variations in the surface energy imbalance. Our tests suggest that when the phase difference was eliminated, latent heat fluxes were enlarged due to the increased contributions from large eddies and the closure was substantially improved. Therefore, the phase difference and not the co-spectral amplitude between vertical velocity and scalars of large eddies explains the non-closure. The phase difference at low frequencies is associated with largest time lags between vertical velocity and water vapor density. Mean advection and entrainment effects can introduce such time lags. In the case of advection, the locally sensed large-scale water vapor density may be lagging the locally produced vertical velocity turbulent velocity excursion due to new water vapor sources contributing from farfield horizontal transport by the mean longitudinal velocity. Likewise, entrainment fluxes introduce lags between water vapor density and locally generated vertical velocity at low frequencies. These lags are mainly due to changes in q 0 as large eddies in contact with a source (land) and a sink (ABL top) experience long turn-over times (analogous to changes in storage). Our tests also suggest that eliminating phase difference increased sensible heat flux, leading to an improved closure under daytime unstable conditions. Thus, we speculate that, for the sites with different Bowen ratios from our site, phase difference between large eddies of wind speed and scalars is also likely to cause non-closure. Last, we emphasize that the work here is not suggesting a phase shift correction be applied to post-field data processing. Such phase shifts in large eddies is not a measurement error to be corrected; it is inherent to the canonical structure of ASL turbulence. Environ. Res. Lett. 12 (2017) 034025