1 Introduction

1.1 Context

First published as an assumption to facilitate a discussion on the spectrum of turbulence, Taylor’s ‘frozen turbulence’ hypothesis (Taylor 1938) (hereafter TH) has since become a fundamental theoretical assumption behind turbulence measurement and theory. In his own formulation, Taylor (1938) argues that if the turbulent velocity is small compared to the averaged wind velocity that transports that turbulence, then:

$$\begin{aligned} u=\phi (t)=\phi \left( \frac{x}{U} \right) , \end{aligned}$$
(1)

where u is a turbulent velocity, x and t are distance and time, and U is the Eulerian wind velocity. I.e., when the turbulence is considered ‘frozen’ as it moves past a stationary sensor, the time domain can be exchanged for the spatial domain via the mean wind speed. This equation thus allows the spatial structure of turbulence to be determined from time-resolving instruments. TH is employed wherever a turbulent time series is used to infer spatial characteristics of turbulence, such as when computing surface turbulent fluxes of heat, moisture, and trace gases via eddy covariance (Foken et al. 2012). Should TH not strictly apply to these measurements, then this may introduce errors in to surface flux calculations, which may contribute to the long-standing problem of surface energy balance closure (Foken 2008; Mauder et al. 2020). Because of its importance and utility, TH has been investigated since the 1950s theoretically (e.g., Lin 1953; Lumley 1965; Wyngaard and Clifford 1977), in laboratory conditions (e.g., Wills 1964; Krogstad et al. 1998; Buxton et al. 2013), in the field (e.g., Panofsky et al. 1974; Mizuno and Panofsky 1975), and more recently in simulation (e.g., Horst et al. 2004; del Álamo and Jiménez 2009).

Field investigations of TH are traditionally difficult and expensive as they require measurement arrays that are both spatially and temporally dense. Generally these investigations have been conducted with sparse arrays of anemometers arranged in one direction (Panofsky et al. 1974; Powell and Elderkin 1974; Horst et al. 2004) or in two axes (Han et al. 2019). Recently, advances in sensor technology have created instruments capable of measuring turbulent variables densely enough in both space and time to be used to investigate TH in the field, such as Raman lidar for humidity fluctuations (Higgins et al. 2012) and fibre-optic (FO) distributed temperature sensors (DTS) for temperature fluctuations (Cheng et al. 2017; Han and Zhang 2022). Thomas et al. (2012) first demonstrated the ability of FO sensors to resolve atmospheric turbulent structures at sub-metre and \(\approx \)0.5 Hz resolution. Since then, DTS has been employed to study many turbulent flows where both the spatial and temporal domains of atmospheric flows are of interest (Zeeman et al. 2015; Peltola et al. 2021; Zeller et al. 2021; Thomas and Selker 2021). The advent of such sensors has re-ignited discussion of TH by allowing field studies at finer scales.

Explicit in Eq. 1 are many assumptions, among them

  1. 1.

    Eddies travel, on average, at the mean Eulerian wind velocity.

  2. 2.

    Eddies of all sizes travel, on average, at the same velocity.

  3. 3.

    There is no dependence on height of assumptions 1 and 2.

  4. 4.

    There is no dependence on dynamic stability of assumptions 1 and 2.

We will use these assumptions as a framework to guide this work. Taken together, assumptions 1 and 2 create the broadest possible interpretation of TH. Namely that for Eq. 1 to be true, eddies must all travel at the same velocity, and that velocity must be equal to the mean Eulerian wind velocity. It would be possible to satisfy one of these assumptions yet still fail the general assumption of TH, e.g., if all eddies travel at the same velocity but this velocity is not equal to the mean Eulerian wind. It is pointed out by Taylor himself that these assumptions obviously cannot be true at all scales and conditions, and much work has been done to investigate the spatial and temporal scales and turbulent conditions over which TH is valid. It is also of course true that any individual eddy may be travelling faster or slower than the mean Eulerian wind as it moves past a sensor, and that the question is whether there exists a systematic bias upon ensemble averaging of many realisations. Panofsky et al. (1974) and Mizuno and Panofsky (1975) found that TH was applicable to smaller eddy sizes but that larger eddies travelled faster than the local mean wind, and Davison (1974) proposed structures advect according to the velocity at the centre of gravity of the eddy. This has remained the general ‘textbook’ understanding of TH and conditions have been broadly acknowledged in which TH is valid, such as turbulence intensity \(I < 0.5\) or \(I^2 < 0.1\) (Willis and Deardorff 1976; Stull 1988). Some theoretical work, such as the anisotropic cascade models described by Schertzer et al. (1997) has argued that a single scale-independent velocity is insufficient to describe turbulent flows. In a laboratory setting Krogstad et al. (1998) found, via a fixed and a moving probe, that instead it is smaller eddies that travel slower than the local mean and to which TH does not apply. This finding has been supported by recent fieldwork in the atmospheric surface layer (Cheng et al. 2017; Han et al. 2019). Cheng et al. (2017) note that the apparent decreased velocity of smaller eddy sizes is attributed to the loss of their coherence due to turbulent diffusion and hence violate of the ‘frozen’ assumption of TH. To fulfill this requirement eddies must remain coherent for the time it takes to advect past a sensor, i.e., for a given eddy size the ratio of a dissipation timescale to an advective timescale must be much greater than 1 (Lin 1953; Everard et al. 2021). It is clear therefore that assumption 2 is incorrect, though the task of defining the relationship between the spatial scale of eddies and the local wind is incomplete. Additionally, for spatial scales over which eddies do seem to travel at a uniform velocity, the local wind is likely not the best estimate of this (Buxton et al. 2013; Han et al. 2019), violating Assumption 1.

Assumptions three and four have been somewhat investigated in the atmospheric surface layer. Mizuno and Panofsky (1975) used two measurement heights of 1.8 and 6.7 m above a lake and found that agreement between eddy velocity and Eulerian wind velocity improved with height in unstable conditions. Davison (1974) found that in a convective layer the advective velocity at 3.5 m represented the Eulerian wind at 100 m agl. Han et al. (2019) used four measurement heights at or below 5 m agl and found no significant changes with height in near-neutral stabilities. Higgins et al. (2012) found that measurement height itself was important to the scale of structures that are measured, and thus will affect the applicability of TH. Cheng et al. (2017) find no change in behaviour between unstable and near-neutral cases. Han and Zhang (2022) specifically investigate the applicability of TH in stable, near-neutral, and unstable cases, and conclude that TH fails in highly turbulent unstable cases. In laboratory studies some variation is found between the ratio of advective to Eulerian velocity with distance from the wall, but as Han et al. (2019) point out, the variation observed in laboratory conditions is less than 5% over the full boundary layer depth (Atkinson et al. 2015; Krogstad et al. 1998; del Álamo and Jiménez 2009). There seems to be an interaction therefore between assumptions 3 and 4: atmospheric stability and the applicability of TH changing with measurement height, but this interaction remains somewhat unclear.

1.2 NamTEX

The Namib Turbulence Experiment (NamTEX) campaign of March 2020 was conducted specifically to investigate some of the theoretical assumptions behind turbulent flow and (sub-)surface—atmosphere heat transfer in the surface layer Hilland et al. (2022a). A field site in the Namib desert (23.516\(^\circ \) S, 15.089\(^\circ \) E) was selected that fulfills the requirements that themselves underly the assumptions behind textbook descriptions of boundary-layer turbulent flow: homogenous and flat fetch on the order of many kilometres and a flat, homogenous surface. The desert site had little to no moisture, no vegetation, and provided a wide range of atmospheric stabilities. The consistency and repeatability of the diurnal meteorological conditions at the site (see Hilland et al. 2022a, Fig. 3) allows for the observation and ensemble averaging of many cases for any given atmospheric condition(s). The NamTEX measurement domain was a flat area \(300 \times 300\,\textrm{m}^2\). The surface cover was densely packed fine sand (Fig. 1A) and the surface aerodynamic roughness length was determined empirically from a sonic anemometer profile to be 0.88 ± 0.36 mm (Hilland et al. 2022b). Using this flat, unvegetated, homogenous field site in combination with spatial- and temporal-resolving instruments, we systematically investigated the applicability of TH with respect to the four guiding assumptions by comparing spatially integrated and spatially scale-dependent eddy velocities to traditional single-point wind speeds.

2 Methods

Fig. 1
figure 1

Photos from the NamTEX field site. A Aerial view showing central mast and X footprint of the DTS system. B One of the dense section poles showing 16 dense transect heights and the cable routing. C One of the long section stakes showing cable routing around the support disk

2.1 Field Methods

2.1.1 Fibre-Optic Temperature Measurements

A Silixa Ultima-S (Silixa, Elstree, UK) fibre-optic (FO) distributed temperature sensor (DTS) was used to sample the turbulent air temperature field. A 3000-m FO cable (AFL, Duncan SC, USA) with 9 \(\upmu \)m core and 1.6 mm outer diameter (teal-coloured PVC sheath) was arranged in two vertical cross-sectional planes in SE–NW and SW–NE orientations between opposing corners of the domain. The DTS provided temperature measurements at 0.76 Hz and each 0.254 m along the cable. Each spatially-integrated 0.254 m portion of FO cable is hereafter referred to as a segment.

Fig. 2
figure 2

Schematic of the DTS array (not to scale) showing the SW-NE axis in red and the SE-NW axis in blue. Information on the calibration baths is provided in Appendix 1. A more thorough description of the cable routing is available in Hilland et al. (2022b). The ‘dense’ section extends between ±13.5 m from the origin, and the ‘long’ section between ± 163.5 m from the origin

Two cable arrangements are distinguished within each cross-section: a ‘dense’ section centred within the domain of 16 measurement transects, each 27 m long horizontally (\(\approx \)106 measurement segments per transect) and separated vertically by 0.17 m from 0.28 m above the surface to 2.85 m, and a ‘long’ section of 327 m horizontal distance (\(\approx \)1287 measurement segments per transect) with cables at two measurement heights: 0.21 and 0.38 m agl (Fig. 2). The dense section provides a ‘grid’ of measurements enabling a small vertical cross-section of the temperature field to be observed, while the long section samples a more spatially extensive area.

The FO cable was suspended on white plastic disks (diameter 0.17 m) bolted either to aluminum poles spaced 9 m apart in the dense section (Fig. 1B), or to iron stakes spaced 10 m apart in the long section (Fig. 1C) with the cable looped completely around each support disk on the latter. The tension of the cable was controlled by hand at least once per day. During the night of 9–10 March the lower SW long transect was cut by an animal chewing through the cable. This break was spliced in the field, however the field splice failed to hold the necessary tension and the tension was largely removed from this transect. The lower height of the long transect is therefore not used for the present analysis.

The temperature spectrum captured by the DTS in a representative day-time 30-min period (16 March 1330–1400 UTC) is shown in Fig. 3 along with the sonic temperature spectrum of one of the sonic anemometers (see Sect. 2.1.2). Figure 3 shows two different measurement locations in the array which demonstrate the increasing noise level along the length of the FO cable. The dense measurement location (earlier in the cable) captures more of the turbulent temperature signal and has lower contribution from noise than the measurement location of the long section (later in the cable). Also shown is the average spectrum within the baths before the measurement array, demonstrating that only noise is measured within the baths. The same spectrum is measured in the baths after the measurement array, albeit with significantly higher noise (not shown).

To shield the measurement unit it was placed inside a large white environmental enclosure and actively ventilated with two small fans and the enclosure itself was shielded from direct sunlight. The system clock of the measurement unit was actively synced using a GPS mouse (NaviLock NL-8012U, NaviLock, Berlin, Germany) and the NMEATime2 software (v. 2.2.1, VisualGPS).

Fig. 3
figure 3

Temperature spectra for a representative 30-min period. \(T_{SA}\)—sonic air temperature from the 3-m sonic anemometer; \(T_{DTS}\) (dense)—average spectrum of measurement segments at 2.85 m agl; \(T_{DTS}\) (long)—average spectrum of measurement segments at 0.38 m agl (only central dense-overlap section); \(T_{Baths}\)—mean spectrum of both warm and cold baths before the measurement array

2.1.2 Sonic Anemometer Measurements

The centre of the domain was defined as the intersection of the two DTS axes. 2 m east of this intersection was a 21-m mast on which a vertical profile of four CSAT-3 sonic anemometer-thermometers (Campbell Scientific, Logan UT, USA) were mounted at 0.5, 1.0, 2.0, and 3.0 m agl. An IRGASON combined sonic anemometer-thermometer and infrared gas analyser (Campbell Scientific, Logan UT, USA) was mounted at 12.0 m. The azimuth of all instruments was 222\(^\circ \) and the measurement volume was < 1 m from the domain centre. All 5 instruments sampled a 3d wind vector (m s\(^{-1}\)) and acoustic air temperature (\(^\circ \)C) at 60 Hz and output at 20 Hz. All instruments were logged on a Campbell Scientific CR1000X datalogger. The logger clock was synced to GPS time via a Garmin GPS16X-HVS (Garmin, Olathe KA, USA).

2.2 Processing

2.2.1 Calibration and Mapping

A warm (ambient) and cold (ice-water mixture) bath were maintained at the field site to provide reference fibre temperatures for calibration. The setup of the baths is elaborated in Appendix 1. The raw fibre temperatures were calibrated as per des Tombe et al. (2020) using the dtscalibration python package (des Tombe and Schilperoort 2020).

The resultant calibrated cable temperatures are not air temperatures due to the significant convective and radiative errors on the cable. Instead we will refer to calibrated fibre temperatures \(T_{\textrm{f}}\). An energy balance model is not applied in the current analysis as the focus is on correlations of turbulent temperature fluctuations rather than accurate calculation of air temperatures. An energy balance model adapted from Sigmund et al. (2017) has previously been applied to these data (Hilland et al. 2022a), but the low-frequency signal imposed on the time series due to small changes in wind direction when near-parallel with the fibre was found to impair, rather than improve, the subsequent correlation analysis.

Converting the length along fibre output of the DTS to a useful coordinate system was performed by indexing the measurements based on average day- and night-time temperatures. The thermal imprint of the supporting structure was clearly visible as regularly spaced abnormally warm (cool) spikes in averaged day (night) temperatures along the fibre. The centre of the support disks was taken as the local peak of these spikes. In the long section where the cable loops completely around a support disk two measurement segments were removed at each disk: the local temperature peak itself was removed and the adjacent points averaged together. The small (2.6 cm) error incurred due to the difference in instrument spatial resolution and disk circumference was accounted for by shifting subsequent points until the next support. Points were mapped to a convenient local cartesian coordinate system in which the SW–NE transects lie along the x axis and the SE–NW transects lie along the y axis with the origin at the central cable intersection. After the spatial mapping, any remaining measurement segments at each support structure that showed strong temperature influence from the disks were discarded from the analysis, totalling 88 segments on the SE-NW axis and 93 segments on the SE–NW axis. A temperature field that interpolates these discarded segments is shown in Fig. 4. More detail on the mapping procedure can be found in Hilland et al. (2022b).

Fig. 4
figure 4

Calibrated DTS cable temperature at one measurement height (0.38 m) over a 10-min sampling period (16 March 1330–1340 UTC) for the SW–NE axis (top) and SE-NW axis (bottom). Each column is a spatial sample, each row a time series from a single measurement segment. At length along cable 0 the transects overlap. Temperatures that were influenced by the support structure have been removed and interpolated for illustrative purposes, causing some persistent horizontal artifacts. Mean Eulerian wind direction \(\overline{\Omega }\) is 272\(^\circ \) ± 3\(^\circ \) at 2.0 m agl and wind speed \(\overline{U}\) is 4.4 m s\(^{-1}\) ± 0.9 m s\(^{-1}\). Lines A-D show segment pairs in the along-wind direction (A and B; C and D). The inset figure shows how segments A-B and C-D are correlated using Eq. 3 and produce \(C_T\) values at \(r'\approx 97\) m in Fig. 6. This same case is shown in Figs. 3, 5, and 6

2.2.2 Case Selection

Continuous DTS measurements began on 8 March at 0700 UTC and ended on 17 March at 0440 UTC, producing 1283 10-min measurement periods. Ten minutes was chosen to reach a compromise between capturing the full spectrum of turbulence and having more measurement periods and was suitable to capture the spectral peak at the site in day-time unstable conditions such as those in Fig. 3. These periods were filtered using a variety of selection criteria removing any 10-min period in which: the cable was broken; an incomplete period was monitored (due to, e.g., a power supply interruption, a hard-drive swap, etc.); the periods during and immediately after the replenishing of ice in the cold calibration bath; the turbulent temperature fluctuations as measured by the fibre-optic cable were below the noise level of the instrument (i.e., the space-time standard deviation of temperature within the measurement array was at or below the space-time standard deviation within the calibration baths of 0.4 K); the wind field was non-stationary, following the method described by Foken and Wichura (1996) using a threshold of variance of 30%. Following this initial filtering, 597 10-min periods remained. After constructing the correlation fields described in the following sections, an additional filter was applied to remove periods where the standard deviation of the space-time correlation field was less than 0.035; a noise-equivalent threshold determined empirically. Following this, 335 periods remained (26% of the full campaign).

Some characteristics of these periods are summarised in Table 1, where \(\overline{U}\)—mean Eulerian wind speed, I—turbulence intensity, \({T\!K\!E}\)—turbulent kinetic energy, \(\zeta \)—Obukhov’s stability parameter, and \(Q_H\)—turbulent sensible heat flux. Though a few nocturnal periods are included, the strong temperature forcings required to produce a signal in the FO cable bias the dataset towards day-time strongly unstable conditions. A wide range of wind speeds and degrees of instability are represented, and the turbulence characteristics generally conform with textbook ‘rules-of-thumb’ for the applicability of TH.

2.3 Analysis

2.3.1 Eulerian Wind

The advection velocities are compared against a mean Eulerian wind calculated from the sonic anemometer profile. For each analysis period the Eulerian wind at each of the DTS measurement heights was calculated by logarithmic interpolation to the mean winds of the lower four sonic anemometers (i.e., linear best fit of \(\overline{U}\) with ln(z)). The DTS heights (28 cm to 285 cm) are well represented by four co-located sonic anemometers (50 to 300 cm).

2.3.2 3-Dimensional Space-Time Correlation

The common one-dimensional normalised space-time correlation function used to determine advective velocities can be written as:

$$\begin{aligned} C_s (r, \tau ) = \frac{\overline{s'(x+r, t+\tau )\;s'(x,t)}}{\sqrt{\overline{s'(x+r,t+\tau )^2}\;\overline{s'(x,t)^2}}}, \end{aligned}$$
(2)

where x and t are position and time, r is a spatial separation, and \(\tau \) is a time lag. The variable \(s'\) is any turbulent flow variable. In field studies where the flow can not be experimentally controlled and repeated, a horizontal line of sensors is typically used and analyses restricted to flows whose mean direction is approximately along the line of sensors (e.g., Cheng et al. 2017; Everard et al. 2021). At the NamTEX field site the prevailing wind direction rotated consistently counter-clockwise through 360\(^\circ \) each day. With the 2-axis array correlations can be determined independent of wind direction by considering the two-dimensional horizontal separation of each pair of measurements. Replacing \(s'\) with fluctuations of calibrated FO cable temperatures \(T'_{\textrm{f}}\):

$$\begin{aligned} C_T (r, \upsilon , \tau ) = \frac{\overline{ T'_{\textrm{f}} (x+r, y+\upsilon , t+\tau ) \; T'_{\textrm{f}} (x,y,t) }}{\sqrt{\overline{T'_{\textrm{f}}(x+r,y+\upsilon ,t+\tau )^2}\;\overline{T'_{\textrm{f}}(x,y,t)^2}}}, \end{aligned}$$
(3)

where x is position along the SW-NE axis and y the position along the SE-NW axis, r and \(\upsilon \) are spatial separations along the respective axes, and \(C_T\) is the normalised two-point correlation coefficient bound between -1 and 1. For the current analysis only horizontal spatial separations are considered. Iterating over all segment pairs not removed during mapping (i.e., those not influenced by the support structure) between the two axes yields in almost all cases two segments with identical (to the cm) r and \(\upsilon \) separations. Due to imperfect spacing between support structures, some (r, \(\upsilon \)) combinations only have a single unique match and such segments are filtered out to preserve symmetry about the origin. All other (overlapping within 1 cm) measurement pairs are averaged together.

For a single time lag \(\tau \), Eq. 3 produces a 2d spatial field of \(C_T\) consisting of \(\approx \)1.45 million correlation pairs. The gaps in the field left by the removal of measurements distorted by the support structure are bilinearly interpolated to complete a spatially continuous \(C_T\) field (Fig. 5A). This process is repeated for a series of time lags from \(\tau =-70\) s to \(\tau =70\) s, producing a ‘correlation cube’ (Fig. 5B).

For each analysis period the mean Eulerian wind direction \(\Omega \) as calculated from the central sonic anemometer at 2.0 m was used to slice the \(C_T\) cube. The result is a 2d space-time \(C_T\) field aligned along the mean Eulerian wind direction (Fig. 5C) in which the ordinate is a \(\tau \) from \(-70\) to 70 s and the abscissa a spatial separation \(r'\) with variable spatial resolution and maximum distance, depending on the wind direction. The locus of points formed by following the local maxima of \(C_T\) through space-time indicates the integral advective velocity \(\overline{U_a}\) (Sect. 2.3.3) of eddies passing through the array (Dennis and Nickels 2008; Atkinson et al. 2015).

Fig. 5
figure 5

Illustration of the ‘space-time correlation cube’ construction and slicing using the long transect at 1330–1340 UTC 16 March. A The \(C_T\) field at \(\tau =0\). B Stacking \(C_T\) fields from \(\tau =-70\) to \(\tau =70\) s. C The cube slice showing \(C_T\) in the direction of \(\Omega \) (only positive space lags \(r'\) shown). The dashed white line qualitatively indicates the \(C_T\) peak through space-time, whose slope is proportional to the integral eddy velocity \(\overline{U_{a}}\)

Table 1 Summary conditions for the 335 filtered periods

2.3.3 Integral and Scale-Dependent Advective Velocity

An integral advective velocity was determined for each 10-minute period and each height by tracing the \(C_T\) peak through (\(\tau \), \(r'\)) space (Fig. 5C). The peak of each \(r'\) row was found by taking the highest \(C_T\) value in the row and fitting a polynomial around the maximum value and adjacent 6 values. The peak of this polynomial was used to find the \(\tau \) of maximum correlation at sub-second precision. Points with \(C_T>0.1\) were kept, and points with \(C_T<0.1\) were discarded. A linear regression was performed on these points with a forced-0 y-intercept (\(y=mx\) model), as \(\tau \) must equal 0 at \(r'=0\). The slope is the integral advective velocity \(\overline{U_{a}}\) in m s\(^{-1}\). Goodness of fit was determined via the mean absolute error (MAE) in s.

Fig. 6
figure 6

Illustration of the integral velocity method applied to the 10-minute segment 1330–1340 UTC 16 March. Left – the \(C_T\) field with \(r'\) on the ordinate and \(\tau \) on the abscissa. Right – the axes swapped. The purple line indicates the average Eulerian wind speed (4.4 m s\(^{-1}\)). The green dots show the peak values of \(C_T(\tau ,r')\) for each \(r'\). Black dots indicate peaks which pass the \(C_T>0.1\) threshold. The green line shows the linear best fit through the black dots, the reciprocal slope of which is the integral DTS velocity \(\overline{U_a}\) (4.7 m s\(^{-1}\)). In this case of good agreement between the two, \(\overline{U_a}/\overline{U}=1.1\) and the MAE of the linear fit is 1.9 s

An example of the integral method is shown in Fig. 6, in which \(\overline{U_a}\) and \(\overline{U}\) are similar and MAE of the linear fit is small. This example demonstrates the difficulty of the NamTEX measurements: a \(C_T\) cutoff must be defined to avoid the noisy correlation field at larger \(r'\), but in doing so many observations that are qualitatively valid are eliminated, and decrease the overall quality of the measurement.

A scale-dependent advective velocity was determined by calculating a series of slopes between all pairs of row peaks \(r'=0\) and \(r'=i\) for each row of \(r'\) in which the peak of the row was \(C_T>0.1\). I.e., the slopes between \((\tau ,r')=(0,0)\) and all black points in Fig. 6. The green points are discarded and no slopes are calculated to these. Fluctuations below the temporal resolution of the instrument were filtered by removing measurements where \(r'<(1/0.76)\cdot \overline{U}\). The FO cable acts as a low-pass filter with a delayed time response as shown in Fig. 3 and discussed in, e.g., Peltola et al. (2021) or Sayde et al. (2015). However, a wind speed derived from correlations between FO measurement segments requires only that both cable segments respond equally, on average, to the wind field, and therefore separations up to the instrument temporal resolution are used (see discussion in Appendix 2). If eddies of all sizes travel at the same velocity, this collection of slopes should remain equal. This approach does not isolate single wavenumbers or eddy sizes, but does nevertheless incorporate the effect of scale: by averaging together the velocities calculated at increasing distances of \(r'\), only structures of sufficient size will remain correlated while smaller structures are filtered out (Krogstad et al. 1998). Each 2-point velocity is normalised against the average Eulerian wind velocity, as interpolated at the DTS transect measurement height. Scale-dependent velocities were calculated for all heights and both dense and long sections.

3 Results

3.1 Correlations

The absolute values of \(C_T\) calculated along the integral best-fit lines are generally low. In Fig. 7 the ensemble median \(C_T\) values along the linear best-fit path are shown separated by stability classes. The median case for all stabilities reaches the cutoff of 0.1 by about 80 m spatial separation and in all cases \(C_T\) decreases approximately exponentially with increasing \(r'\). Most cases fall in to the very unstable class (\(n=1872\) cases) or the unstable class (\(n=3085\) cases) and there is significant variation in \(C_T\) within these classes. \(C_T\) is generally higher in the very unstable and unstable classes. \(C_T\) decreases most rapidly in the near-neutral class (\(n=276\) cases), suggesting that in cases which pass all the exclusion criteria, near-neutral cases lose spatial coherence most rapidly. The number of cases described in Fig. 7 refers to unique measurement heights and 10-minute periods, i.e., 16 measurement heights for one 10-minute period comprise 16 cases. Additionally, because there is only one measurement height that extends the full distance of the array, there are only four stable cases at \(r'>15\) m which produces the large noise in the stable class.

Fig. 7
figure 7

Median \(C_T\) values along the integral velocity line of best fit for four different stability classes: \(-\infty <\zeta \le -1\) (very unstable); \(-1<\zeta \le -0.1\) (unstable); \(-0.1<\zeta \le 0.1\) (near-neutral); \(0.1<\zeta <\infty \) (stable). For very unstable and unstable classes the 10th and 90th percentiles are also shown. \(C_T\) is also shown for the cold bath for (\(r'>0\))

3.2 Integral Advection Velocity

The integral advective velocities of the DTS (\(\overline{U_a}\)) are compared against the mean Eulerian wind velocities measured from the central sonic anememoter profile (\(\overline{U}\)). In Fig. 8\(\overline{U_a}\) and \(\overline{U}\) are displayed for a single measurement height (285 cm, dense section). The median \(\overline{U_a}\)/\(\overline{U}\) is 1.1 and inter-quartile range is 0.9 to 1.4. The distribution of the ratios has a right skew which biases the mean of the ratios higher to 1.3. The mean absolute error of the fits is overwhelmingly below 1 s, with a median MAE of 0.5 s.

Despite low MAE, and suitably filtered conditions, there remain a significant minority of cases (22) in which \(\overline{U_a}\)/\(\overline{U}\) exceeds 2 at this measurement height. These cases tend to be caused by a very strong correlation field centred around \(\tau \approx 0\) that spans the full range of \(r'\). This causes the locus of peaks through the \(C_T\) field to remain centred around \(\tau \approx 0\) and the advective velocity therefore to appear quite high. Due to the explicit spatial filtering caused by the physical array these may be an artifact of much larger waves or eddies that the array can not capture.

Fig. 8
figure 8

Left - integral-scale advective velocity of \(T_f\) (\(\overline{U_a}\)) compared with the mean Eulerian wind \(\overline{U}\) for one measurement height (285 cm agl). Right - histogram of the ratio of the left panel. For clarity the axis of the right panel is clipped, removing three outliers at \(\overline{U_a}\)/\(\overline{U}>4\). 283 10-minute cases are displayed. TH is supported when \(\overline{U_a}=\overline{U}\)

The \(\overline{U_a}\)/\(\overline{U}\) distribution shown in Fig. 8 repeats for all measurement heights (not shown): a median of 1.05 ±0.1; inter-quartile range of 0.85 ±0.1 to 1.35 ±0.1; and a right skew caused by a small number of very high \(\overline{U_a}\) values. This pattern also emerges when applied to the larger spatial extent of the long DTS section, though in this case the median \(\overline{U_a}\)/\(\overline{U}\) is slightly higher at 1.2 (IQR 1.05 to 1.50) and the MAE also increases to 1.6 (IQR 1.05 to 3.5) s. At all measurement heights the \(\overline{U_a}\) and \(\overline{U}\) distributions were tested with a paired two-sample t-test which found the distributions to be significantly different at \(p<0.001\) except for 215 cm agl, where \(p=0.014\).

3.3 Scale-Dependent Advection Velocity

The 2-point velocities from Sect. 2.3.3 are each individually normalised against \(\overline{U}\) at the transect height as described above. The collections of individual velocities are binned by logarithmically increasing intervals of \(r'\). In Fig. 9 the binned distributions are shown as boxplots where the box extent is \(Q_1\) to \(Q_3\), the median shown within the box, and whiskers extend to \(\pm 1.5 \times \)IQR or the maximum/minimum value in the distribution.

Fig. 9
figure 9

Distributions of 2-point velocity ratios binned at logarithmically increasing intervals of \(r'\). Box height is from \(Q_1\) to \(Q_3\) and the bold horizontal line the median. Whiskers extend to \(\pm 1.5 \times \)IQR. The horizontal lines at 1.0 indicate TH. Top—grouped by stability classes as per Fig. 7. Bottom—grouped by I quartiles at 2.0 m agl. Note that outliers (exceeding ±1.5 times the IQR) are disabled for these plots for clarity and in each panel comprise <2% of measurements

In Fig. 9 (top panel) the boxplots of velocity ratios are divided by stability using Obukhov’s stability parameter \(\zeta =z/L\) and four stability classes. Distributions were tested against one another for statistical similarity using the non-parametric two-sample Kolmogorov-Smirnov (KS) test. Due to the higher wind speeds associated with near-neutral stabilities, there are no valid velocity ratios at \(r'<5\) m for this class. For all other stability classes, the velocity ratios are generally below 1 at smaller \(r'\), reach a median ratio of 1 around 5–9 m, and exceed 1 at \(r'>28\) m. The exception is \(r'\) bin 16–28 m in which the very unstable class median is slightly below 1. There are variations in the relative paths traced along \(r'\) for each stability class, but they all follow the general trend that at smaller \(r'\) the medians are below 1, increase with increasing \(r'\), and eventually cross 1 and remain above 1 for larger distances. The stable class in the largest \(r'\) bin is much higher in median velocity ratio (\(>2\)) than the other classes. This is likely biased by the small number of stable cases and even fewer 2-point velocities that pass the \(C_T>0.1\) threshold to produce the distribution. The behaviour of the stable class at \(r'>\) 28 m can be safely assumed to be a failure of the method rather than reflective of extremely rapid eddies.

The very unstable and unstable classes are similar at \(r'<5\) m (KS \(p>0.01\) for the first three distance bins) and maintain similar median ratios up to \(r'<16\) m, after which the unstable class stabilises while the very unstable class shows greater fluctuation in median and range ratio. In the 85 m \(<r'<\)150 m bin the median ratio is 1.36 for the very unstable class, 1.22 for the unstable class, and 1.26 for the near-neutral class. The distributions in the 85 m \(<r'<\) 150 m bin are all significantly different from each other at \(p<0.001\).

In Fig. 9 (bottom panel) the velocity ratios are shown grouped by I quartile. As before, median \(\overline{U_a}/\overline{U}<1\) at small \(r'\) for all classes before reaching and exceeding 1 around 5 to 16 m. Perhaps surprisingly \(Q_4\), the highest turbulence intensities, remains closer to a ratio of 1 at greater distances, with a median ratio of 1.02 at 16 m \(<r'<\) 28 m compared to 1.15 and 1.16 for \(Q_2\) and \(Q_3\), respectively. The three distributions \(Q_2\) to \(Q_4\) are significantly different from each other at \(p<0.01\) in this bin. The large IQR for the lowest I class (\(I<0.1\)) at larger \(r'\) is again likely caused by the method failing in calm conditions. At lower separation distance 1.5 m \(<r'<\) 3 m there is statistical similarity between \(Q_1\) and \(Q_2\) and between \(Q_3\) and \(Q_4\) (\(p>0.05\) for both). In the 85 m \(<r'<\) 150 m bin the median ratio is 1.20 for I \(Q_2\), 1.24 for \(Q_3\), and 1.24 for \(Q_4\). Despite similar medians, the distributions in the 85 m \(<r'<\) 150 m bin are significantly different from each other at \(p<0.001\).

In Fig. 10 the 2-point velocity ratios are examined for each of the 16 heights within the dense section of the DTS for near-neutral cases (left) and all unstable (including very unstable) cases (right). Rather than presenting the full distributions, the fraction of \(\overline{U_a}\)/\(\overline{U}>1\) per bin is shown. I.e., at 0.5 the median \(\overline{U_a}\)/\(\overline{U}\approx 1\). Only 10-minute periods for which all 16 measurement heights passed the exclusion criteria are shown: 30 cases in near-neutral conditions and 198 cases in unstable conditions.

The near-neutral cases are disordered and there is no obvious relationship between measurement height and median velocity ratio. In the unstable cases there is a clear relationship between transect height and the median velocity ratio. At \(r'>2\) m the lowest transect heights are consistently 10 to 15% below the highest transect heights. The higher transect heights reach and exceed a median \(\overline{U_a}\)/\(\overline{U}>1\) at smaller \(r'\) separations.

Fig. 10
figure 10

Velocity ratios binned at 1-m \(r'\) intervals by height within the dense DTS section. For each bin, the fraction of \(\overline{U_a}\)/\(\overline{U}>1\) is shown, i.e., whether the median velocity ratio in the bin is above or below 1. Sixteen transect heights are evenly spaced between 28 and 285 cm. Left – near-neutral cases (\(-0.1\le \zeta <0.1\)). Right – unstable cases (\(\zeta \le -0.1\)). Only periods in which all 16 measurement heights pass the exclusion critera are considered

4 Discussion

Assumption 1

States that eddies travel at the mean wind velocity. Alternatively: the integral-scale advection velocity determined from the DTS should equal the Eulerian wind velocity of the sonic anemometer. The results shown in Fig. 8 and the results of the \(\overline{U_a}\)/\(\overline{U}\) distributions for the other measurement heights as well as the full extent of the long section do not support A1. At all heights except for one (215 cm agl) where the median was 0.95 the median \(\overline{U_a}\) exceeds 1. The median \(\overline{U_a}\)/\(\overline{U}\) of 1.05 suggests that eddies travel on average 5% faster than the Eulerian wind. However, when the full extent of the DTS array is considered, the median ratio increases to 1.2, i.e., eddies travel 20% faster than the Eulerian wind. This is similar to the value of 16% found by Han et al. (2019). The better agreement between \(\overline{U_a}\) and \(\overline{U}\) in the dense DTS section is likely the result of spatial filtering caused by the comparatively small spatial extent of this section in relation to the mean size of eddies. From Fig. 9 the median \(\overline{U_a}\) is very close to \(\overline{U}\) at spatial separations of 5–16 m in all conditions, and this is commensurate with the spatial extent captured by the dense DTS section. Furthermore as seen in Fig. 9, the IQR of 2-point velocity distributions decreases with increasing \(r'\) as the calculated velocities become less sensitive to minor variations, lending greater confidence to results obtained from the larger spatial extent.

Assumption 2

States that eddies of all sizes travel at the same velocity. The correlation-cube method used here does not allow isolation of individual wavenumbers as has been employed in previous studies that examine scale-dependent advection velocities (Cheng et al. 2017; Buxton et al. 2013). Instead, collecting distributions of 2-point velocities imposes a variable spatial filter on the array and acts as a proxy for spatial scale-dependence. While each individual velocity calculated at a given spatial separation along the mean-wind direction \(r'\) integrates all turbulent scales that the DTS can measure with the given separation, patterns emerge on ensemble averaging at different separation distances that indicate effects of spatial scale. If all spatial scales of turbulence travel at the same velocity, then the ratio of \(\overline{U_a}\) to \(\overline{U}\) should stay the same across a series of spatial separations.

A common pattern that emerges in this analysis is that at small separation distances \(\overline{U_a}<\overline{U}\), i.e., the advective velocity is less than the Eulerian velocity. This is observed in all stability classes except near-neutral (Fig. 9—top panel), all I classes (Fig. 9—bottom panel), and at all measurement heights (Fig. 10). This suggests that smaller turbulence scales are advecting slower than other turbulence scales. This pattern has been observed in laboratory and atmospheric flows (Krogstad et al. 1998; Han et al. 2019).

It is however also true that at larger separation distances \(\overline{U_a}>\overline{U}\), i.e., the advective velocity is greater than the Eulerian velocity. This is similarly observed in all stability classes, all I classes, and at all measurement heights. In some cases the median \(\overline{U_a}\)/\(\overline{U}\) remains essentially unchanging above some \(r'\): e.g., in I \(Q_2\) and \(Q_3\) the median ratio varies by less than 0.05 at \(r'>28\) m. In these conditions all but the smallest eddies appear to be travelling at the same velocity, but this velocity is not equal to the Eulerian velocity. In these cases there is a characteristic “global” \(\overline{U_a}\) that is a better estimation of the eddy velocity as others have proposed and which may be better suited to use as a transport velocity when applying TH (Buxton et al. 2013; Krogstad et al. 1998; Han et al. 2019; Dennis and Nickels 2008).

In other conditions however there is no such asymptotic approach of a characteristic velocity. In the very unstable class shown in Fig. 9 the median velocity ratio continues to increase over the spatial extent of the array, reaching 1.36 in the largest bin. This would be consistent with the centre-of-gravity hypothesis (Davison 1974) in which an increasingly unstable and large surface layer generates larger eddies which travel more rapidly. In this case the very unstable class should also asymptote towards a characteristic velocity but the NamTEX array was of insufficient size to capture this.

Assumption 3

States that there is no effect of height on the ratio of advective to Eulerian velocity. From Fig. 10 there is a subtle but clear gradient in the velocity ratios with height in unstable conditions. There is no such pattern in near-neutral conditions. Closer to the surface the median \(\overline{U_a}\) remains below \(\overline{U}\) until a higher spatial separation is reached compared with measurements farther from the surface. Note that the analysis presented in Fig. 10 does not describe by how much the velocities disagree, but only whether the median is greater or lower than \(\overline{U}\).

This is consistent with previous findings that there is no clear effect of height on the behaviour of TH in the near-neutral surface layer (Han et al. 2019) as well as findings that there is a change with height in the unstable surface layer (Panofsky et al. 1974). This supports the conclusion that in certain conditions measurement height must be considered when applying TH (Higgins et al. 2012), and that higher measurements are more influenced by larger-scale and therefore faster moving eddies.

Although the array represents only a tiny fraction of the full boundary layer height, reaching only 3 m above the surface, this is the region with the greatest velocity gradient with height, especially given the very smooth surface. Even within this small extent a change with height is noticeable. This suggests that in unstable conditions, the applicability of TH will depend on measurement height.

Assumption 4

States that there is no effect of dynamic stability on the ratio of advective to Eulerian velocity. Although each stability class presents the same general distribution of lower-velocity small scales and higher-velocity large scales, there is evidence of systemic differences between them.

The most important difference between the very unstable and unstable classes, as mentioned above, is that the unstable classes appear to asymptote to an approximately constant median velocity ratio. In contrast the very unstable class continues to increase through the spatial extent of the NamTEX array. Not just stability class, but also degree of instability affects the behaviour of TH: in very unstable conditions, quite large spatial separations may be needed to accurately determine the characteristic advective velocity of the flow.

While the stable class clearly deviates at larger separation distances, this can be attributed to the loss of coherence of the flow at these larger separations, after which even though some observations manage to pass the \(C_T>0.1\) threshold these are not representative of a well-correlated field. However this analysis suggests that where the correlations remain sufficiently large, stable surface layer flows behave similarly to other stability classes. Important to note of course is that the sample of stable flows that passed the various thresholding criteria is very small and represents a very small subset of well-correlated stable flows, rather than all stable flows during the observation period.

Ultimately the NamTEX DTS array was not sufficiently large in either horizontal or vertical extent to fully investigate TH and a larger array in both dimensions could provide more insight. Furthermore, it must be reiterated that most stable cases from the campaign were discarded and conclusions drawn from the small subset of remaining stable cases should be limited.

This analysis finds that assumptions 1 through 4 do not strictly apply. It begs instead the question: when, if ever, does TH apply? As stated by Everard et al. (2021) the utility of TH for both theory and practice, past and present, is not contested. This analysis suggests however, that there is only a small window of spatial separation, which is likely site-specific, for which TH is strictly valid. At the NamTEX site this is \(5<r'<16\) m, a separation over which \(\overline{U_a}\approx \overline{U}\) in almost all cases. Coincidentally this is the spatial extent covered by the dense section of the DTS, which, with some variance, produces an integral advective velocity ratio of 1.05. Outside of this window TH does not strictly describe the observed flows.

Any theoretical or practical application in which the spatial domain of turbulence is inferred from the temporal will be distorted to some extent by the failure of TH. An application of particular interest is the underestimation of turbulent heat fluxes calculated from the eddy-covariance method, which relies on TH to infer the spatial structure of turbulence (Foken et al. 2012). It is generally known that the eddy-covariance method underestimates turbulent heat fluxes, and a failure of TH has long been speculated to contribute to this (Mahrt 2010): If larger eddies, which usually contribute the most to the turbulent sensible heat flux signal, travel systematically faster than smaller eddies that do not transport heat as efficiently, then their time-fraction in a time-series measured by a single-point measurement is underestimated in comparison to their true spatial extent, which defines the landscape-averaged flux. Indeed, evidence suggests significant underestimation of fluxes where TH fails (Cheng et al. 2017), and correcting this error may contribute to the long-standing problem of surface energy-balance closure (Foken 2008; Mauder et al. 2020). Attempting to improve energy-balance closure by applying a correction for the space-time distortion caused by TH remains an important task and subject for future work, though this analysis suggests a single global correction may not be possible, as any correction would depend on stability, height, and potentially other controls on co-spectral distributions such as roughness and boundary layer heights.

5 Conclusions

This work is an attempt at a systematic investigation of Taylor’s Hypothesis (TH) at an ideal ‘textbook’ desert field site using temperature as a passive scalar. A two-axis array of fibre-optic cable was installed within the surface layer providing temporally and spatially dense temperature measurements: up to 300 m horizontal extent and a smaller section of 16 measurement heights and 27 m horizontal extent between 0.28 and 3 m agl. The two fibre axes were correlated with each other in space and time, creating 3-dimensional correlation fields which were sliced along the mean wind direction. Both integral and scale-dependent eddy advection velocities were then calculated and compared against the mean Eulerian wind at the same height.

At the integral scale, eddies travelled faster than the local wind, with the median velocity ratio 1.05 when considering spatial separations up to 14 m. At larger spatial separations the median integral velocity ratio increased to 1.2, indicating eddies travelled on average 20% faster than the mean Eulerian wind. Scale-dependent advective velocities were calculated using 2-point correlations. In all cases eddies were found to travel slower than the mean Eulerian wind at smaller separation distances, and faster than the mean Eulerian wind at larger separation distances. Depending on stability, the largest eddies travelled on average 22% to 36% faster than the mean Eulerian wind.

Four guiding assumptions implied by Taylor’s original equation were used to structure the work, and none of the four was supported by the analysis. Namely eddies of all sizes do not, on ensemble averaging, travel at the mean Eulerian wind velocity, nor at the same velocity. In unstable conditions there is a distinct influence of measurement height on the applicability of TH, and it does not behave uniformly across different atmospheric stabilities. Furthermore, although our array was able to capture and demonstrate a characteristic advection velocity that may be used instead of the mean Eulerian wind for the application of TH, in very unstable cases no such velocity was observed.

Understanding of when, and to what extent, TH is either applicable or not, has significant implications for any application in which the spatial domain is to be inferred from the temporal. Of special interest is turbulent heat flux calculations and impact on surface energy-balance closure. This analysis suggests that a global correction for the failure of TH is not possible, and a correction to TH must consider additional variables such as atmospheric stability and height above the surface.