1 Introduction

In the surface layer, the transport of momentum, temperature, and trace gases between the surface and the atmosphere is dominated by turbulent coherent structures (Finnigan 1979, 2000; Raupach 1981; Raupach et al. 1996). Therefore, understanding the behaviour of coherent structures can help improve our interpretation of surface fluxes and advance knowledge of fundamental transport physics. Conditional sampling broadly describes analyzing turbulent signals to identify patterns of coherent motion (Antonia 1981; Raupach et al. 1996). Some popular conditional sampling methods include quadrant analysis (Wallace et al. 1972), visual event detection (Chen and Blackwelder 1978; Gao et al. 1989; Paw et al. 1992), variable-interval time-averaging (VITA) (Blackwelder and Kaplan 1976), the window-averaged gradient approach (Bisset et al. 1990), and analyses via wavelet transforms (Farge 1992; Collineau and Brunet 1993a, b; Thomas and Foken 2005).

Quadrant analysis is one commonly applied method for understanding the nature of turbulent coherent structures (Wallace et al. 1972). It involves splitting the velocity field into turbulent perturbations via Reynold’s decomposition in the vertical and streamwise directions and displaying them as joint-probability distributions for a given averaging period. Wallace et al. (1972) named events in quadrant IV “sweeps” because they indicate the downward motion of higher momentum, positive streamwise velocity perturbations in the flow. Likewise, events in quadrant II are “ejections” because they represent the upward motion of lower momentum fluid. In near-surface flow, such as the atmospheric surface layer, sweeps and ejections represent down-gradient momentum fluxes, and quadrants I and III represent counter-gradient momentum fluxes. Willmarth and Lu (1972) and Lu and Willmarth (1973) introduced the concept of using a hyperbolic hole in quadrant analysis, where only events with Reynolds stress magnitudes greater than a selected threshold are counted for each quadrant.

It has been shown that sweeps dominate in the roughness sublayer, defined as the surface to approximately three times the canopy height (Kaimal and Finnigan 1994). Above the roughness sublayer, in the inertial sublayer, sweeps and ejections are of roughly equal importance (Finnigan 1979; Raupach 1981). Although quadrant analysis is a useful method for understanding the statistical frequency of coherent structures, it does not allow for structures to be connected in time or space and does not resolve individual coherent structures.

Other expansions of quadrant analysis include octant analysis (Suzuki et al. 1988) and a trajectory analysis technique (TRAT) (Nagano and Hishida 1989; Nagano and Tagawa 1995). Octant analysis involves including an additional variable (typically either the cross-stream velocity component or a scalar) to represent spanwise momentum or vertical heat fluxes for example. Using the TRAT approach, one can trace the motion of the turbulence as coherent structures cross between quadrants, and thus can be used to associate structures in time between quadrants. Like a hodograph, the TRAT approach preserves the history in defining the coherent structure. Wallace (2016) provides a thorough review of the expansions and applications of quadrant analysis in turbulence research.

Visual event detection is commonly used to check other conditional sampling methods, such as those based on scalar signals (usually temperature as in Chen and Blackwelder 1978) resulting from coherent structure velocity fields. A gradual rise followed by sharp drop in the temperature signal (a temperature ramp) may indicate an ejection of warmer air followed by a sweep of cold air associated with a coherent structure under unstable conditions (Gao et al. 1989; Paw et al. 1992). Hereafter, the term microfront is used to describe the sweep indicated by a temperature drop (Gao et al. 1989). Conversely under stable conditions, inverse temperature ramps have been observed (Shaw et al. 1989; Gao et al. 1992). Temperature ramps are not evident when there is no vertical temperature gradient or when there is insufficient turbulence.

The VITA (Blackwelder and Kaplan 1976) method is a jump-detection algorithm that compares the total variance in a time series to the variance of a windowed signal (e.g., velocity or temperature) (Blackwelder and Kaplan 1976). The VITA technique can be employed with a temperature signal, however, historically it was calculated using the wind speed signal. When used with the temperature signal, VITA identifies periods of high variance, which should correspond to a microfront. The signal produced from this technique is like a low-pass filter. When the value of VITA surpasses an arbitrary threshold, k, and its slope is negative, identification of a coherent structure is triggered.

Wavelet transforms are one of the most commonly used methods for identifying the presence and characteristics of coherent structures. Wavelet transforms provide information on the scale of the structure, in addition to the location in time of the events. Haar and Mexican hat wavelet transforms are commonly used in turbulence research (Farge 1992; Collineau and Brunet 1993a), because they roughly resemble the shape of turbulent signals and because the signals are localized in time. Nonetheless, turbulence signals rarely closely follow the shape of wavelets. Another limitation of wavelet analyses are leakage problems with orthogonal mother wavelets. These can occur either when the signal does not match the shape of the wavelet function or because of the translation of the signal during the transform (Qiu et al. 1995).

In this study, the continuous Mexican hat wavelet (MHAT) applied to the velocity and temperature signals is selected to indicate the presence of a coherent structure for integrated quadrant analysis (IQA). The method from Collineau and Brunet (1993b) for the MHAT wavelet is employed. The wavelet scale was selected based on the expected number of structures in the time period. The structure is defined when the transform crosses the x-axis. Like the VITA method, there is a slope criterion that is selected based on the expected direction of the signal jumps.

The final coherent structure identification method described here is the use of the rotated speed, \(U_r\), introduced by Shaw et al. (1989). This method has the benefit of identifying coherent structures when scalar tracers are not present. The rotated speed is calculated from a coordinate transformation based on the rotation angle (\(\theta \)) for an averaging period that is, in turn, based on \(u'\) and \(w'\)

$$\begin{aligned} \theta =0.5\tan ^{-1}{\left( \frac{2\overline{u' w'}}{\overline{u'^2}}-\ \overline{w'^2}\right) }. \end{aligned}$$
(1)

The rotated speed \(U_r\) is calculated from a coordinate transformation of the \(u'\) and \(w'\) velocity fields based on \(\theta \), and it represents sweeps and ejections in one term. The rotated speed has been used to identify coherent structures in two ways: (1) VITA was calculated with the rotated speed and (2) \(U_r\) can be integrated with respect to time to represent the strength of a sweep and ejection (Shaw et al. 1989).

There is no one accepted definition of turbulent coherent structures, nor is there one accepted method to identify them. Despite many options for using conditional sampling to identify the locations and physical characteristics of turbulent coherent structures, all the methods are limited by the assumption that temporal coherence implies spatial coherence. Calculating a coherence spectrum (Davenport 1961) provides the statistical spatial coherence, but it lacks the time information and cannot identify each individual structure’s physical shape. Other non-parametric methods including singular spectrum analysis (Vautard et al. 1992) and empirical orthogonal function analysis (Finnigan 2000), offer alternative tools for analyzing the statistical and periodic characteristics of turbulent coherent structures.

Because of the limitations in quadrant analysis and other conditional sampling methods, the IQA technique was developed to be used in addition to a new conditional sampling method to connect sweeps and ejections in time and space and better characterize turbulent coherent structures. The new conditional sampling method uses a measure that is the square root of twice the turbulence kinetic energy (\(uTKE= \sqrt{\overline{u^{\prime 2}}+ \overline{v^{\prime 2}}+ \overline{w^{\prime 2}}}\)) to identify the existence of a coherent structure in the ejection phase, before the temperature drop and sweep. The IQA method connects Eulerian point measurements to the Lagrangian motion of the air parcel by integrating the turbulent velocity field over time. Unlike the TRAT approach, the IQA technique connects coherent structures in physical space instead of by quadrant space. Here, the IQA method and the conditional sampling technique to identify the presence of turbulence coherent structures are described (Sects. 3 and 4). In Sect. 5.1, the conditional sampling method is compared to other conditional sampling and jump detection methods. In Sect. 5.2, the IQA technique is employed for two field experiments: one in a canopy and one in an open field above the roughness sublayer.

2 Field Experiment Sites

The IQA technique was tested using turbulence data from two field campaigns: the Eclipse Boundary Layer Experiment (EBLE; Higgins et al. 2019) and the Vertical Array Cherry Experiment (VACE). The EBLE campaign was conducted in a grass field (canopy height \(h \approx \) 150 mm, aerodynamic roughness length \(\approx \) 20 mm) with ultrasonic anemometers located at heights 10 to 30 times the roughness length, measuring the turbulence far above the roughness sublayer. Conversely, the VACE campaign was intentionally designed to measure turbulence at several heights within the canopy and within the roughness sublayer for a broader representation of IQA performance.

The EBLE campaign was conducted on 19–22 August 2017 before, during, and after the Great American Eclipse in Corvallis, Oregon (Higgins et al. 2019). Data from two ultrasonic anemometers on the south-central tower in a horizontal array (Tower A1 in Fig. 1 from Higgins et al. (2019)) were used; an ultrasonic anemometer (81000Re, RM Young, Traverse City, Michigan, USA) was mounted at 3.5 m above the ground and a gas analyzer/ultrasonic anemometer (IRGASON, Campbell Scientific Inc., Logan, Utah, USA) was mounted at 1.56 m. During the experimental campaign, the predominant wind direction was from the north-east where the surface roughness was uniform for the distance of more than 1 km, greater than the typical measurement fetch.

The VACE campaign was operational from November 2019 through July 2020 in Linden, California. Five ultrasonic anemometers, (CSAT3, Campbell Scientific Inc., Logan, Utah, USA) sampling at 20 Hz, were mounted on scaffolding between two cherry trees in an orchard with a north–south row orientation. The row spacing was 7 m and the tree spacing was about 6 m, with tree heights of 3.6 m. The anemometers were mounted facing \(310^{\circ }\) from north towards the predominant wind direction (north-west) adjacent to the crown of a cherry tree. They were located at heights of 1.45 m, 2.2 m, 3.0 m, 3.9 m, and 5.8 m above the ground. The anemometers located in the crown were placed between 0.20 m to approximately 1 m from the branches and leaves. The site also included an operational evapotranspiration station, so there were also sensors measuring net radiation, ground heat flux, soil moisture, and soil temperature.

The VACE campaign involved three intensive measurement periods: 8 November through 25 November 2019 for the decaying late fall leaf period; 9 January through 1 February 2020 for the no-leaf period; and 29 March through 10 April for the leaf-out period. In the November period, leaves remained on the trees for the entire period. Only the autumn study period was used to develop and verify the IQA method. Figure 1 shows the experimental set-up of the VACE campaign during the November 2019 period.

Fig. 1
figure 1

Experimental design for the November 2019 period of the VACE campaign (top panel). The left photo is the state of the canopy in November 2019 looking into the primary wind direction. The photo was taken between the top two ultrasonic anemometers. The right photo shows the position of the sonic anemometers

Full size image

3 Integrated Quadrant Analysis

The IQA method builds on traditional quadrant analysis by integrating the streamwise (u), cross-stream (v), and vertical velocity (w) perturbations with respect to time, after rotation into the mean wind direction for a 30-min averaging time. The recursive integration is described by

$$\begin{aligned}&X_i=\ X_{i-1}+\left( u_i^\prime \ \right) \delta t \end{aligned}$$
(2)
$$\begin{aligned}&Y_i=\ Y_{i-1}+\left( v_i^\prime \ \right) \delta t \end{aligned}$$
(3)
$$\begin{aligned}&Z_i=Z_{i-1}+\left( w_i^\prime \right) \delta t \end{aligned}$$
(4)

where \(X_i\) is the position for the streamwise component, \(Y_i\) is the position for the cross-stream component, \(Z_i\) is the integrated position for the vertical velocity component, and \(\delta t\) represents the timestep between samples of the measured data (e.g., 0.05 s for a 20 Hz sampling rate). The variables \(u_i^\prime \), \(v_i^\prime \), and \(w_i^\prime \) are the streamwise, cross-stream, and vertical velocity perturbations from the mean, respectively. The subscript i indicates the velocity field for one time. The integrated positions can be plotted on an X–Z plane or an X–Y–Z volume to visualize the trajectory of the measured air parcels.

Unlike traditional quadrant analysis where the quadrants represent sweeps, ejections and two interaction quadrants, with the IQA technique, flow trajectories are split into two regimes: bulk sweep and bulk ejection. In the IQA method, bulk sweeps are defined as periods when \(Z_i<\ 0\), and bulk ejections are defined as periods when \(Z_i>\ 0\). For this reason, bulk sweeps can consist of traditional quadrant analysis’s sweeps, ejections, and interaction periods, as the trajectory moves through three-dimensional space. The integrated trajectory therefore maintains the coherence of the turbulent structure, by describing its positional history as it passes the anemometer, and it yields a quantitative measure of the structure’s size. The values of \(X_i\) and \(Y_i\) do not contribute to the definition of bulk sweep or ejections. The position of \(X_i\) indicates the sign of the momentum transport with respect to the mean wind speed. The position of \(Y_i\) indicates a change in direction in the cross-stream component. The slope of the IQA line in Fig. 2 represents the trajectory.

Fig. 2
figure 2

An example of the IQA results plotted on an X–Z plane with artificial data. The structure is split into parts: the bulk sweep and the bulk ejection

Full size image

The IQA describes the motion of a parcel in both time and space (Fig. 2) given artificial data in two dimensions. One can imagine that the sensor is located at the origin on the plot. The structure starts at the origin as a bulk sweep (pink), and the parcels translate down with relatively fast wind. The green points having closer proximity indicate a decreasing wind speed. At times 2 and 3, the structure slows, and parcels move towards positive Z. Although the direction of the parcel indicates an ejection in traditional quadrant analysis, the parcel has not crossed the Z = 0 plane, and so it is still in the bulk sweep regime. At time 4 (blue points), the structure passes the height of the sensor and then is reclassified into the bulk ejection regime. This bulk ejection regime continues until the parcel moves below the height of the sensor, and the trajectory cycle represents one coherent structure.

Unlike quadrant analysis where the interaction quadrants are often interpreted as non-coherent turbulence, in the IQA method, the interaction times are incorporated in the definitions of the bulk events as they can account for up to \(40\%\) of the cycle. Moreover, the IQA technique provides the advantage of connecting sweeps and ejections in time. For example, one bulk ejection consists of a sweep followed by an ejection. The structures can be delineated; the sweep and ejection are components of one event. Although the IQA method is like the TRAT approach in that they both connect structures in time, the IQA method integrates the velocity perturbations so that the structure is resolvable in space through its trajectory.

To initialize the IQA technique, a coherent structure should be identified as having began, otherwise the solution drifts with an arbitrary starting point. The previously discussed conditional sampling methods use either the \(u'\) or \(\overline{u'w'}\) velocity signals or a scalar signal to identify the presence of a three-dimensional coherent structure. Using a detection method that evaluates all three velocity components should lead to more consistent identification of structures than two dimensional methods in the methods tested. This is because in some cases, the \(\overline{v'w'}\) signal is of equal or greater magnitude than the \(\overline{u'w'}\) signal.

4 Triggering the Start of a Coherent Structure

There is no consensus on a method for identifying turbulent coherent structures in a time series, nor do the existing methods for identifying turbulent coherent structures always identify the same coherent structures. A new conditional sampling method was developed to be used in conjunction with the IQA technique. This new conditional sampling method is compared to visual identification of temperature ramps, VITA, wavelet transforms, and a multi-level rotated velocity signal. These methods were described in detail in Sect. 1. Figure 3 is a schematic to help describe both the conditional sampling method and the IQA method that accompanies the explanation in this section. Hereafter, the word trigger is used to denote the identified start of the coherent structure at the measurement location for the purpose of the IQA technique.

Fig. 3
figure 3

Summary of IQA procedure. a Calculate \(\sqrt{u^{\prime 2}+v^{\prime 2}+w^{\prime 2}}\) for high-frequency data. Brown dashed lines are visually identified temperature ramps. b, c Apply low-pass filter to high-frequency \(\sqrt{u^{\prime 2}+v^{\prime 2}+w^{\prime 2}}\) d Apply MHAT wavelet to \(uTKE_{LP}\). e Trigger the start of the IQA procedure for coherent structures at minimum of the MHAT coefficients when the amplitude passes a threshold. f Calculate IQA starting with the integration constant re(set) to zero at each IQA trigger time g Display results of the IQA method in two or three dimensions

Full size image

The triggering method that is described here to initiate the IQA method employed the square root of the sum of the high-frequency velocity variances, which is also twice the square root of turbulence kinetic energy (uTKE). It also incorporates the turbulent cross-stream velocity component. Whereas the previously described methods use scalar signals that are not always present, or one- or two-dimensional velocity fields, large-eddy-simulation studies suggest that the cross-stream velocity component is important for maintaining a microfront (Fitzmaurice et al. 2004; Finnigan et al. 2009); therefore all three velocity components are crucial to coherent structures.

To calculate uTKE, the velocity data were first transformed via rotation into the mean wind direction (u) for a 30-min averaging time. A tilt correction was not appropriate for the VACE campaign because mean vertical velocities were probable; however, it was applied for the EBLE data. Because uTKE is axis invariant, the value of uTKE is not dependent on the selected coordinate rotation. After calculating the high-frequency components, the quantity uTKE was low-pass filtered by integrating the signal over a moving time window of size \(\Delta t\)

$$\begin{aligned} uTKE_{LP}= \frac{1}{\Delta t} \int _{t - \Delta t / 2}^{t + \Delta t / 2} \sqrt{(u'^2+v'^2+w'^2)}dt . \end{aligned}$$
(5)

Conceptually, this is analogous to a moving average multiplied by the window time interval. Integrating uTKE provides a measure of the distance of the structure in three directions. The length of time of the moving average might depend on the height of the canopy, but in this case, we found for both experiments that a 10-s integration time appears appropriate. This integration time was determined empirically by altering the integration time to match the peaks with the temperature ramps. The integrated uTKE signal trace resembled the MHAT wavelet (evident in Fig. 3c).

Next, the continuous MHAT wavelet (Torrence and Compo 1998) was applied to the integrated uTKE signal. The selected scale for the MHAT wavelet was determined empirically by comparing the signal of the coefficients at several scales to the temperature signal for several time periods with clearly identified ramps. For both the VACE and EBLE data, the selected scale of the wavelet was 10. This corresponds to a wave with a period of approximately 40 s. The period selected may be linked to the height of the canopy and the expected duration of temperature ramps.

The wavelet period was checked by calculating the temperature ramp amplitude, spacing, and duration using structure functions (van Atta 1978) with two lags (Paw U et al. 2005; Shapland et al. 2012). For the VACE campaign, the duration of the temperature ramps varied from 5 s to almost 1 min depending on the height above ground level and the time of day. Conversely in the EBLE campaign, the temperature ramps had shorter duration, so a scale that corresponds to a wave period of 30 s was selected. For both experiments, the duration and separation of the ramps were highly variable over time and space because of the multi-scale nature of turbulence.

At the determined scale, the amplitudes of the MHAT coefficients were used to trigger the detection of a coherent structure. The amplitudes of the wavelet coefficients depend on the measurement height and the mean uTKE for the period. The effect of sensor height on the amplitude was small and was accounted for using a linear regression to normalize amplitudes to the amplitudes at a reference height. The impact of the mean uTKE on the mean amplitudes was more important as it created a diurnal variation in the mean amplitudes. Therefore, this method needs to be calibrated for each dataset. The details of the calibration will follow in Sects. 4.1 and 4.2.

The IQA procedure is summarized in Fig. 3. After calculating the high-frequency uTKE, (a) it is low-pass filtered using an integration time of 10 s (b). The MHAT wavelet transform is performed on the quantity \(uTKE_{LP}\) (c). The amplitude of the MHAT coefficients is used to trigger the identification of the start of a coherent structure at a local minimum of the signal (d). The IQA is calculated with the \(u'\), \(v'\), and \(w'\) signals between the identified coherent structures (green lines on Fig. 3d) in the MHAT coefficient signal. Finally, IQA can be displayed in two or three-dimensions to visualize the trajectory of one coherent structure (g).

4.1 Trigger Calibration: Vertical Array Cherry Experiment

For the VACE campaign, approximately one week of data from 8–15 November (n = 324 30-min periods) was used to normalize the mean amplitudes to create a threshold for triggering coherent structure detection. For each 30-min period, the MHAT wavelet was applied to the high-frequency uTKE signal, and the amplitudes of the coefficients from minimum to maximum were calculated. The number of wavelets per period varied from 40 to 45. The amplitudes were averaged over each 30-min period to calculate mean wavelet coefficient amplitudes. The mean amplitudes were used to calibrate this coherent structure identification method. The general family of equations for the wavelet coefficient amplitude is

$$\begin{aligned} A_i=m_i A_{Z_i}+ b_i . \end{aligned}$$
(6)

The slope from Eq. 6 is a function of the relative height, \(z_i = (5.8 - z)\) where z is the measurement height (in metres) of the anemometer and 5.8 m is the height of the reference anemometer. Therefore, \(z_i\) is the distance below the reference height. The reference height was selected because it was approximately 1.5 times the canopy height. In the VACE campaign, the intercept is not a function of \(z_i\), however, in some cases, there may be an intercept effect, so it is checked for each dataset. The equation for the slopes and intercepts based on height are

(7)
$$\begin{aligned} b_i = \mu \times z_i + \beta . \end{aligned}$$
(8)

In Eq. 7, (Burmese letter “ma”) represents the slope, and (Burmese letter “ba”) represents the intercept of the \(m_i\). In Eq. 8, \(\mu \) and \(\beta \) represent the slope and intercept respectively of \(b_i\) for all of the sensor heights.

Equations 7 and 8 can be substituted into Eq. 6 to account for the impact of height on the amplitudes of the MHAT wavelet coefficients:

(9)

For the VACE dataset, there was no impact on the intercept with height, so \(\mu \times z_i + \beta \) was set to the average intercept (\(\approx \) 2.05). There is likely to be a relationship between the intercept with height in taller canopies.

In addition to the mean amplitude being a function of height, it is also a function of the turbulence. The results from the linear regression between the average amplitude at the reference height and the value of uTKE at the same height for each 30-min period can be described with

(10)

where \(A_{Z_{Ref}}\) is the average MHAT amplitude at the reference height. In Eq. 10, (Slavic Cyrillic letter “em”) represents the slopes of the regressions between the amplitude of the MHAT coefficients at reference height and any variable, and (Slavic Cyrillic letter “be”) represents the intercepts of the regression between the amplitude of the MHAT coefficients at the reference height and any variable.

By substituting Eq. 10 into Eq. 9 with \(A_{Z_{Ref}}\) for \(A_{Z_i}\), one can diagnose the mean MHAT amplitude at each level \(A_i\)

(11)

Table 1 has the coefficients to calculate \(A_i\) as a function of uTKE at the reference height in roughness sublayer (\(Z_{Ref}\) = 5.8 m/h = 3.6) for both the VACE and EBLE data.

Table 1 Calibration coefficients for VACE and EBLE experiments for the coherent structure identification technique
Full size table

Figure 4 shows the mean amplitude for each 30-min period against the calculated average of the mean amplitudes of the MHAT wavelet coefficients. The lowest heights are overestimated relative to the higher heights. There is also a visible deviation from the linear regression line at low amplitudes, where \(A_i\) is overestimated. This is due to periods of low values of uTKE (less than 0.3 m \(\mathrm{s}^{-1}\)). When these periods are removed, the slope is 0.9 instead of 1. Because of this bias, this method may not work well under low uTKE periods, which are less likely to have shear-induced turbulent coherent structures than high uTKE periods.

The amplitude of \(A_i\) represents the expected mean amplitude in each 30-min period, and it can be used to trigger the identified start of coherent structures. The amplitudes from the MHAT wavelet are not normally distributed over a 30-min. To optimize the triggering of identified individual events, the amplitude of the wavelet must surpass 1.25\(A_i\) for the VACE campaign. This threshold was determined empirically so that the identified coherent structures best matched the visually identified temperature ramps.

Fig. 4
figure 4

The MHAT wavelet amplitudes for the VACE campaign corrected by height and by uTKE at the reference height. Sonic 1 is the anemometer closest to the surface (z = 1.45 m), and Sonic 5 is the highest anemometer (z = 5.8 m)

Full size image

Table 2 displays the average number of times a coherent structure is identified per 30-min period separated by local time (UTC-8). There are more identifiable coherent structures triggered during the day than during the night. The analysis indicates that if the coherent structures were evenly spaced during the night-time, a structure would occur approximately once every 3 min. During the day, coherent structures occur closer to once every 2 min. In addition to a diurnal effect on the occurrence of coherent structures, the highest number of structures occur near the top of the canopy. This supports the results from Gao et al. (1989), who found that the strongest temperature ramps occur near \(80\%\) of the canopy height. This suggests that coherent structures are formed near the top of the canopy, and not all of them penetrate completely through the canopy roughness sublayer to the lower canopy layers near the surface (Paw U et al. 2005; Thomas and Foken 2005). The data from the lowest anemometer show relatively frequent coherent structures, which could be due to the overestimation of the value of Ai relative to the amplitude that is evident in Fig. 4.

Table 2 Average number of triggers per 30-min period using \(A > 1.25 A_i\) for the VACE campaign using the IQA coherent structure identification technique
Full size table

When the amplitude of a given MHAT wavelet passes the threshold of 1.25\(A_i\), a coherent structure is identified starting at the minimum of the wavelet coefficient’s wave. The minimum tends to occur preceding the microfront in the temperature signal, and the maximum tends to occur just following the microfront. The event identifications are triggered during the weak ejection phase of the coherent structure (blue points in Fig. 2). After the microfront, the uTKE signal peaks indicating the structure is in the stronger sweep phase (red and pink points in Fig. 2). The events are identified from one trigger until the next trigger, which means that the quiescent period between coherent structures is grouped with the proceeding structure. It would be possible to trigger the end of the structure as well as the beginning (possibly using a \(Z=0\) crossing in the IQA technique) so that one entire structure could be separated with the IQA method.

4.2 Trigger Calibration: Eclipse Boundary-Layer Experiment

Like triggering structures identification in the VACE data, the EBLE data were used to calibrate the linear regression model to correct the data for both height and amplitude. The height correction becomes negligible for sensors located above approximately 1.5 times the canopy height. In the EBLE data, the gas analyzers/sonic anemometer (\(z = 1.56\) m) and sonic anemometer (\(z = 3.5\) m) from Station A1 were used so that there were two heights. Because the measurement height of both sensors greatly exceeds the roughness height, the height dependence on the MHAT amplitudes is negligible in contrast to the VACE campaign.

In this case, the amplitude of the sonic anemometer at \(z = 1.56\) m was corrected by the sonic at \(z = 3.5\) m. Because the slope was approximately one and the intercept was approximately zero, there was little effect from using the height correction in this case. There is little scatter between the mean MHAT coefficient amplitudes at the heights because both sensors are located above the roughness sublayer. The height correction was used with this regression, not with the regression of the changes of slopes and amplitudes with height, as was done with VACE data.

After the height correction, the value of uTKE from the sonic at \(z = 3.5\) m was used to diagnose the mean amplitude at \(z = 1.56\) m. In this case, there is not the low magnitude uTKE cluster on the regression between the corrected and uncorrected MHAT amplitudes. This may be because the sensors in the EBLE campaign were located in the surface layer far from the roughness elements under higher wind speeds. There were fewer low wind speed periods in this experiment than the VACE experiment. As for the VACE campaign, the equation for \(A_i\) is corrected for the height (Eq. 6) and uTKE dependence (Eq. 10). Table 1 presents the coefficients used in the \(A_i\) correction. The threshold of \({A}/{A_i}>1.25\) was used to trigger coherent structure identification.

Table 3 has the average number of events triggered for each height during the entire experiment split by time. There tend to be more coherent structures per 30 min during the day than at night. The lower level has more structures during the day, and they are approximately the same at night.

Table 3 Average number of triggers per 30-min period using \(A > 1.25 A_i\) for EBLE data using the IQA identification technique
Full size table

5 Results and Discussion

The IQA technique was developed and tested with the data from the VACE campaign. The results from the triggering technique were compared with visual ramp identification from the temperature signal, VITA, wavelet identification, and the multi-level rotated velocity approach. To help automate the process and subsample the time series, periods with probable temperature ramps were determined by calculating the amplitude and duration of the temperature ramps using the two-lag method with van Atta’s structure functions (van Atta 1977; Paw U et al. 2005). Temperature ramps were considered to represent likely coherent structures if (1) they occurred at more than one level in the canopy and (2) had an amplitude of near \(1\,^{\circ } \mathrm{C}\) similar to the criteria used by Gao et al. (1989). The visually identified ramps met these two criteria. Periods with large temperature ramp amplitudes (calculated with van Atta’s method) did not necessarily correspond to periods with visually identifiable temperature ramps.

5.1 Triggering Turbulent Coherent Structures

To evaluate the performance of the uTKE trigger method, one 30-min period of VACE data on 14 November 2019 was selected to compare with other methods for identifying a coherent structure. This period was selected because it has strong temperature ramps at multiple levels from 1410 to 1425 LT, and no multi-level temperature ramps at the start of the period. Figure 5 shows the temperature perturbation (\(T'\)) signals for each level in this period. The dashed lines indicate the 13 temperature ramps that were visually identified using the aforementioned temperature ramp criteria.

Fig. 5
figure 5

Temperature perturbation signal for 14 November 2019 at 1400. The dashed lines indicate the visually identified temperature ramps that occur at more than one level and had an amplitude of about \(1\, ^{\circ } \mathrm{C}\)

Full size image

We applied VITA and the MHAT wavelet transformation to both the temperature and u signals for the same period. We determined VITA thresholds (k) by evaluating the results of the VITA function, so it best matched the visually identified temperature ramps. For this period \(k = 0.75\,^{\circ }\mathrm{{C}}^2\) for temperature and \(k = 1\,\mathrm{{m}}^2\, \mathrm{{s}}^{-2}\) for u. The scale of the MHAT transformation was selected so that it would approximately match the number of events detected with the aforementioned method.

Figures 6 and 7 present the \(T'\) and low-pass filtered \(u'\) signal (with the 5-s moving average) marked with the identification results of the times from VITA (blue dashed lines) and MHAT wavelet (orange dashed lines) marked on the graphs. We find that VITA works better for visually identifying the temperature ramps in this case than the wavelet. This is likely because the period contains a regime change after the first 10 min. Applying VITA to the temperature signal understandably replicates the results for the visual identification of temperature ramps. It identifies fewer structures than were visually selected at low levels and more structures at higher levels. This is likely because the ramps had to be identified at more than one level to count as a coherent structure in the visual method. This method successfully identified the strongest temperature ramps between 1414 and 1420 LT. It also signalled some events near 1410 that look like high-frequency turbulence instead of coherent structures related to ramps, but the largest difference between the VITA approach for temperature and the visually identified temperature ramps seems to be the acceptable range of the ramp, which is a function of the threshold. Because the wavelet transformation evenly spaces events, VITA performs better on both wind speed and temperature than the wavelet transformation as the coherent structures were irregularly spaced in this period. Applying VITA to wind speed triggers more events from 1415 to 1420 min than the rest of the time period, but it does indicate that there are some events that occur before 1410.

Fig. 6
figure 6

The temperature signal for 14 November 2019 from 1400 to 1430 LT. The blue dashed lines represent the coherent structures identified with the MHAT wavelet, and the orange dashed lines indicate the structures identified with the VITA temperature

Full size image
Fig. 7
figure 7

The low-pass filtered \(u'\) signal for 14 November 2019 from 1400 to 1430 LT. The blue dashed lines represent the coherent structures identified with the MHAT wavelet, and the orange dashed lines indicate the structures identified with the VITA \(u'\)

Full size image

Using the aforementioned wavelet transformation to identify coherent structures overestimates the number of coherent structures in the period (Fig. 7). The selected wavelet method is not calibrated, so every zero crossing of the transformed wavelet identifies a coherent structure. For this subset of VACE data, using this wavelet method over-diagnoses the number of events that occur in quiescent times and under-diagnoses the events when they are clustered in time. Another reason that the wavelet method does not work as well as VITA could be because the wavelet transformations using one scale identify evenly spaced events. The quiescent periods should be between events, not before the events start. By making the scale shorter, the wavelets trigger the events better in the middle of the period, but also add additional false events near the start of the period.

Because coherent structures are coherent in space as well as time, we compared one multi-level scheme with the single-level identification schemes. Averaging across all heights, this method should identify structures that are strong enough to penetrate the canopy. After calculating the values of the rotated velocity (\(U_r\)), all the levels were averaged (Shaw et al. 1989). Figure 8 shows the results for each method performed on height-averaged \(U_r\) values and layer-averaged \(T'\) values. Applying VITA to the layer-averaged value of \(U_r\) (k = 2) indicates that there are several events between 1405 to 1410 LT that are not evident in the temperature signal. Between 1415 LT and 1420 LT, VITA only indicates three temperature ramps. Integrated \(U_r\) values (k = 2.5 m) (see Sect. 4 for the method of integration) indicate several events near the start of the period, when there is no temperature signal, and it does identify several events between 1415 and 1420 LT that VITA also catches. Although there are more triggered events than the visual identification method, the \(U_r\) integration method does identify the events that are evident only through the velocity signal.

Fig. 8
figure 8

Multi-level rotated velocity analysis for 14 November 2019 at 1400 LT. The x-axis is displayed as minutes after 1400. a The temperature perturbation signal averaged over all sensor heights. b The multi-level rotated velocity with a 10-s moving average. c The VITA detection function marked with purple lines where the detection crosses the threshold k = 2, and d The results from the integrated rotated velocity. The orange lines indicate where the structure cross the distance threshold of k = 2.5 m

Full size image

Triggering the identification of coherent structures with the amplitude \(A_i\) is shown in Fig. 9. The results from the wavelet transform of the low-pass filtered uTKE are shown by the black line. The green lines indicate where the structure identification is triggered with this conditional sampling method, and the dashed brown lines are the visually identified temperature ramps. In this case, the structure identification is triggered just before the temperature ramps, which were visually identified. Like the other methods, there are several events not seen in the temperature signal in the beginning and end of the period, however, they tend to occur consistently across the levels. Unlike using a wavelet alone, this calibration helps to prevent a structure from being identified once per wavelength. Like VITA and unlike using a wavelet without calibration, the \(A_i\) method can handle long and irregular quiescent periods between structures, which should help it perform better in periods with few structures and intermittent turbulence.

Fig. 9
figure 9

The wavelet coefficients (scaled about zero) using the MHAT wavelet with a period of 4-s. The dashed lines indicate the times of temperature ramps and the green lines indicate where the amplitude of the coefficient triggers the start of a coherent structure

Full size image

For this period, there were 13 temperature ramps that were visually identified in more than one level with amplitude of at least \(1\,^{\circ }\mathrm{C}\). Table 4 compares the number of events found in all the previously discussed methods. “Unique Events” refers to the number of events that are triggered at one or more level; a structure that is identified at more than one level is only counted once. The VITA method tends to yield variable results depending on the height, which suggests that the threshold may have to change with height. It also identified few of the same events at different heights. Wavelet methods were more consistent with identifying the same events with height than the other methods; however, this is heavily a function of the wavelet scale selected. For example, the wavelet method cannot identify a coherent structure that is located closer in time to another than the selected wavelet scale. Although this wavelet method does not require an amplitude threshold, calibration is necessary to obtain results comparable to those of established methods. The height-averaged rotated velocity is useful because it considers both the vertical and streamwise velocitie components and ensures that an event occurs at all levels. Both methods that use height-averaged rotated velocity for identifying the structures require a threshold, which makes it more subjective to apply to longer time periods. The rotated velocity methods seem to miss several of the temperature ramps in the middle of the time period.

The method described in Sect. 4 has the advantage that once the equations for \(A_i\) have been calibrated, it can be reasonably applied to any time series from the same experiment. One drawback to the aforementioned method is the need to calibrate the method for each dataset; however, the procedure described in Sect. 4 should be universally applicable. Further research should be done to test the conditional sampling method in other canopies. Another advantage of this method is that it has several event identifications that are not visually identified, but in most cases, they occur across the heights. Moreover, there are few temperature ramps that are not identified.

Table 4 Average number of triggers per 30-min period using \(A > 1.25 A_i\)
Full size table

5.2 Results of Integrated Quadrant Analysis

To demonstrate the results of the IQA technique, one period each from VACE and EBLE campaigns are presented. Times were selected based on identifying temperature ramps that are associated with the coherent structures. One advantage of the \(A_i\) identification method is that it is not dependent on a scalar signal representing the velocity field, although having clear temperature ramps helps to independently confirm the timing of the phases of the sweep-ejection cycle.

For the VACE campaign, 14 November 2020 from 1415 to 1420 LT was selected. Figure 10 shows this 5-min period during which there were several temperature ramps that occurred at multiple levels. Near the end of the 5-min period, there are three temperature ramps with amplitudes of \(0.5^{\circ }\mathrm{C}\) to \(1^{\circ }\mathrm{C}\). In the middle of the time series, there are two temperature ramps that have comparable amplitudes, but do not occur at all levels. The brown dashed lines indicate that the visually determined temperature ramp is evident in two or more levels. One should note that there is likely a structure occurring between 1416 and 1417 LT based on the results of the \(A_i\) method and the velocity signals from Fig. 10; however, there is insufficient evidence of the start of a structure in the temperature signal.

Fig. 10
figure 10

Temperature perturbations at each height on 14 November 2019 from 1415 to 1420 LT (left panel) for the VACE campaign. The microfronts (identified by visual inspection of the \(T'\) signal at Z/h = 1.1) are marked in vertical dashed brown lines. In the right panel, the integrated distance from uTKE and its components with the \(w'\) for the same period of time. The visually identified microfronts are indicated by the vertical dashed brown lines labelled A through F

Full size image

The right panel in Fig. 10 presents the integrated quantities \(\sqrt{u^{\prime 2}+v^{\prime 2}+w^{\prime 2}}\), \(\sqrt{u^{\prime 2}+w^{\prime 2}}\), and \(\sqrt{v^{\prime 2}+w^{\prime 2}}\) for the same time. The peaks in these signals tend to occur after the microfront. This lag in the uTKE signal is not an artefact of the low-pass filter as it is a centred filter, and the data show the same pattern without the low-pass filter. The magnitude of the term with three wind components is highest compared to the components with only u and v, which indicates that the w contributions are important. In the microfronts labelled A, E, and F, \(\sqrt{u^{\prime 2}+w^{\prime 2}}\) accounts for most of the total distance. In the microfronts labelled B, C, and D, the \(\sqrt{v^{\prime 2}+w^{\prime 2}}\) component is of nearly equal magnitude of the instantaneous uTKE signal. In all cases, \(\sqrt{v^{\prime 2}+w^{\prime 2}}\) is of the same order of magnitude as \(\sqrt{u^{\prime 2}+w^{\prime 2}}\). This suggests that even though the data have been rotated into the mean wind direction, the cross-wind component is important in the detection of some of the coherent structures. Furthermore, between ramps A and B, there is an additional peak in the high-frequency uTKE signal that is not matched with a high amplitude temperature ramp that is evident in more than one height. The peak in uTKE suggests the presence of a coherent structure, but this is not supported by the scalar signal.

Figure 11 shows the scaled coefficients of the MHAT wavelet transform for the \(\sqrt{u^{\prime 2}+v^{\prime 2}+w^{\prime 2}}\) signal at the scale corresponding to 40-s for this time period. The brown dashed lines indicate the times of the microfronts from Fig. 5. The green solid lines indicate the times where IQA is triggered. The trigger occurs during the weak ejection phase of the coherent structures, and is followed by the sweep, where \(\sqrt{u^{\prime 2}+v^{\prime 2}+w^{\prime 2}}\) peaks and decreases. Triggering the structure identification at the minimum of the wavelet coefficient allows the ejection to be paired with the subsequent sweep. One should note that the ejection/sweep cycle is continuous, so it is somewhat arbitrary to chose the start of the structure during the ejection phase. However, starting with the ejection phase follows the concept of the weak ejection preceding the stronger sweep in canopies as introduced by Gao et al. (1989). Microfronts B and C are not triggered as coherent structures at all levels. Although an event could have a strong temperature signal, if it does not have a corresponding velocity signal, it may not be identified as a coherent structure.

Fig. 11
figure 11

The wavelet coefficients using the MHAT wavelet on \(\sqrt{u'^2 + v'^2 + w'^2}\) with a period of 40 s. The dashed lines indicate the times of microfront (at the temperature drop) and the green lines indicate where the amplitude of the coefficient triggers the start of a coherent structure event identification as defined in Fig. 10

Full size image

The IQA technique has the advantage of providing the trajectory of a coherent structure fluid elements for one or more specific events. If plotted in two dimensions or animated, it shows the motion of the air parcel relative to the mean advection. For example, when \(u'\) changes sign (and X changes direction), this indicates wind speed relative to the mean. Conversely, when \(v'\) and \(w'\) change sign (Y and Z change direction), it indicates a change in direction of the cross-stream and vertical velocity components, respectively. Figure 12 shows the X–Z plane of IQA for microfront F (Fig. 10), which occurs around 1419.45 LT for the levels at and below Z/h = 1.1. The colours indicate time in 40 s since the trigger. The trigger time shows a very slight lag in the coherent structure as it penetrates the canopy. At Z/h = 1.1, the structure begins 2 s before it is detected at the bottom of the canopy. The duration of this structure varies from about 55 s above the canopy to over 100 s at Z/h = 0.4. At the middle three heights, the time of the structure was approximately 70 s. The structure at Z/h = 0.4 includes more quiescent time than the other heights.

At Z/h = 1.1 and 0.8 (left panels of Fig. 12), the structure begins as a bulk ejection. In both cases, after the first ejection, the X signal moves back toward X = 0, which indicates the presence of relatively fast gusts that accompany the positive vertical velocity component. Before the sweep begins, Z is nearly constant, which indicates that there is a small separation in the vertical motions of sweeps and ejections in this case. Near 25 s after the structure begins, the sweep starts to bring faster moving air down towards the surface. The parcel returns to the origin as the sweep continues, with decreasing intensity, until the parcel crosses the Z = 0 threshold, where it appears that a quiescent period begins before the subsequent structure. At the lower levels, Z/h = 0.6 and Z/h = 0.4 (right panels of Fig. 12), the ejection at the start of the coherent structure is weaker, but the sweep remains strong in both cases, so the parcel stays in the bulk sweep regime for a longer period of time. The amplitudes of Z and X decrease as the coherent structure enters the canopy, indicating that the structures become weaker closer to the surface.

Fig. 12
figure 12

The IQA results for microfront F. The colour represents the time since the start of the structure in 40-s. The event ends when the next event is triggered. Please note that because the event is longest at Z/h = 0.4, the colour scale is different than for the other heights

Full size image

5.3 Integrated Quadrant Analysis in the Inertial Sublayer

Midday (1230 to 1300 LT) EBLE data from 20 August were used to evaluate IQA method above the roughness sublayer. Because measurements from the EBLE campaign are 200 times the canopy height, the temperature ramps are not always clear. For example, Fig. 13 shows the temperature signals and the integrated \(\sqrt{u^{\prime 2}+v^{\prime 2}+w^{\prime 2}}\) and its components. There are visible ramps at 3.5 m that have a mean duration and spacing of approximately 30 s. At 1.56 m, temperature ramps are still visible, however, there is also a lower frequency pattern with a period on the order of several minutes.

Although there is not a strong agreement between the temperature perturbation signals at the different heights, the velocity derived signals peak at the same times. Closer to the surface, \(\sqrt{v^{\prime 2}+w^{\prime 2}}\) is relatively more important than \(\sqrt{u^{\prime 2}+w^{\prime 2}}\) in the events near the middle and the end of the period, but they are on similar orders of magnitudes at both levels.

Fig. 13
figure 13

20 August 2017 from 1230 to 1300 LT was selected to test IQA results for the EBLE data. The top two panels display results for z = 3.5 m and the bottom two panels have results from z = 1.6 m. a and c show temperature perturbations for each layer. b and d show the results for uTKE and its components

Full size image

There were 13 and 12 events that this conditional sampling technique identified for z = 3.5 m and z = 1.56 m, respectively, during this 30-min period. Most of the events were concentrated in the last 10 min of the time period where there is a number of smaller peaks in \(\sqrt{u^{\prime 2}+v^{\prime 2}+w^{\prime 2}}\). When a structure is identified, the parcel starts in a strong prolonged ejection, before eventually returning to a sweep. Parcels remain in the bulk ejection regime for most of the time for these structures. This supports the idea that ejections are more important for momentum flux than sweeps outside of the roughness sublayer.

The IQA was calculated for all components of the velocity during this time period. The results for one event at 1252 LT are displayed in Fig. 14. At z = 3.5 m, the event begins in the bulk sweep regime at 1253 LT, and the event begins with relatively slow air being ejected. During the event, there are periods of higher momentum flux interspersed throughout it. After approximately 30 s, the parcel speeds up relative to the mean velocity and the vertical direction begins to change. The rest of the event consists of relatively fast air moving downward. This period would be labelled a sweep in traditional quadrant analysis. Conversely, with the IQA technique, the sweep is paired with the preceding ejection, so the IQA method identifies when the flow regime returns to a bulk ejection phase.

Closer to the surface, this event is identified about 1 s later. It begins in a relatively weak bulk sweep when high velocity air is brought to the surface. After about 15-s, the parcel enters the bulk ejection regime. The return of the structure after the ejection looks similar to that at 3.5 m: the parcel returns to the origin in the sweep phase. At this level, the parcel does not enter the bulk sweep regime before the next structure is identified. Above the roughness sublayer, ejections occur more frequently and are responsible for more momentum transfer than sweeps.

Fig. 14
figure 14

The IQA results for the event on 20 August 2017 at 1253 LT for both EBLE heights. The colours indicate time since the start of the structure in 40 s

Full size image

6 Summary and Conclusions

Integrated quadrant analysis is a new method that integrates the perturbations of the Reynolds-averaged velocity field with respect to time. It must be used in combination with a conditional sampling method like the one introduced to prevent drift in results. The conditional sampling method introduced to trigger the presence of coherent structures uses integrated uTKE values coupled with MHAT continuous wavelet analysis (Sect. 4). Structures should be detected using three dimensional velocity signals because the cross-wind component is important for maintaining microfronts, as shown in Fitzmaurice et al. (2004). Some microfronts are better described using the vertical and cross-wind velocity components than with the vertical and streamwise velocity components used in traditional quadrant analysis (Figs. 10, 13). The IQA method can be used in combination with the uTKE detection to identify and describe structures that do not have a corresponding scalar (e.g., temperature) ramp (Figs. 9, 11). The temperature signal cannot indicate the presence of a coherent structure when the temperature vertical gradient is weak whereas the uTKE method is not limited by the temperature gradient.

Although there are several methods to trigger identification of the start of a coherent structure in the ejection phase, the IQA technique may be useful for identifying the extent of the coherent structure. For example, the parcel trajectory crossing the Z = 0 axis twice might be appropriate for separating the coherent and quiescent times.

Conditional sampling methods, including the one we introduce, remain limited in light of the multi-scale nature of turbulence. Many methods cannot pick up smaller turbulent structures embedded in larger coherent structures. More work needs to be done to elucidate the dominant scales of turbulence. Non-parametric conditional sampling methods might help to solve this issue. Although the conditional sampling method selects one scale of turbulence, IQA can help to unravel the physical trajectories and size scales of turbulence. It can be used to display the smaller scale sweep/ejection cycles within larger coherent structures.

In summary, the IQA technique is a useful technique for connecting ejections and sweeps together in time. The IQA technique reveals two regimes (bulk sweep and bulk ejection) that are defined by the dominant motion of the structure, but the trajectory towards the Z = 0 represents the return component of the structure. For example, in sweep dominated flows, like the canopy turbulence in the VACE campaign, an event begins in a weak ejection and is followed by a stronger sweep. In ejection dominated flows, like above the roughness sublayer in the EBLE campaign, the event begins in a strong ejection and is followed by a weaker, longer sweep that does not necessarily return to the bulk sweep regime. From this analysis, one can approximate a size scale of the eddies by summing the distance the parcel travels during the coherent structure. The selected coherent structure in the EBLE campaign was larger in both distance travelled and time than the selected structure in the VACE campaign.

Benefits of applying IQA include connecting Eulerian point measurements to the Lagrangian motion of the air parcel, providing a size scale to the coherent structure measured at one point, linking the sweeps and ejection components of one coherent structure, and visualizing the relative motion of a turbulent air parcel in relation to the sensor. Future work should seek to characterize IQA for different canopy types. While we used the IQA method for two canopy types, taller and denser canopies may exhibit different characteristic structures.