Scaling Properties of Atmospheric Wind Speed in Mesoscale Range

: The scaling properties of turbulent ﬂows are well established in the inertial sub-range. However, those of the synoptic-scale motions are less known, also because of the difﬁcult analysis of data presenting nonstationary and periodic features. Extensive analysis of experimental wind speed data, collected at the Mauna Loa Observatory of Hawaii, is performed using different methods. Empirical Mode Decomposition, interoccurrence times statistics, and arbitrary-order Hilbert spectral analysis allow to eliminate effects of large-scale modulations, and provide scaling properties of the ﬁeld ﬂuctuations (Hurst exponent, interoccurrence distribution, and intermittency correction). The obtained results suggest that the mesoscale wind dynamics owns features which are typical of the inertial sub-range turbulence, thus extending the validity of the turbulent cascade phenomenology to scales larger than observed before.


Introduction
The atmospheric boundary layer (ABL), the lowest layer of the Earth's atmosphere varies in thickness, particularly over land, from d ≈ 100 m during the night to d ≈ 2000-3000 m or even more during the day. As the ABL is in contact with the Earth's surface it effectively connects the surface to the rest of the atmosphere and is dynamics are fundamental in the transport and exchange of moisture, heat and momentum with the underlying surface [1,2]. Clearly, the topography of the Earth's surface also plays an important role in the way in which the ABL interacts with the free troposphere. Within the ABL physical quantities such as flow velocity, temperature, moisture, density stratification (ABL stability), etc. are characterized by rapid or large amplitude turbulent fluctuations, and vertical mixing is strong. Moreover, numerical studies have shown that in case of strongly stratified flows a sort of intermittency can be found at large scales in a range of the Froude numbers between Fr ≈ 0.05, 0.3 [3,4].

Analysis of Mesoscale Wind Velocity Data
The data used in this work have been collected at the Mauna Loa Observatory, Hawaii (MLO, 19 • N, 155 • W, 3397 m above sea level). Before performing the analysis, the standard normalization procedure, was applied to the data (by subtracting the mean value of the data and dividing by the standard deviation) in order to obtain a zero mean and unit variance data set. MLO is best known for its long-term monitoring of greenhouse gases and the long range transport and chemical oxidation of other gaseous and pollutants species in the free troposphere (FT). The MLO is located on the northern slope of the Mauna Loa volcano and is mainly subject to two wind components: up-slope (radiative heating) and down-slope (cooling) radiation winds, characterized by a 24-h cycle. The second component, also known as barrier wind, is characterized by the interaction of FT flow with the island itself. The characteristics of the boundary layer and surface winds at the MLO have been studied intensively, mostly due to the necessity to determine the baseline conditions for atmospheric composition studies [19,40]. During day time, the MLO typically measures air masses coming from the the Marine Boundary Layer (MBL), while during the night time it measures air masses coming from the FT [19]. An exhaustive discussion of the wind velocity data acquisition process and instrumental description can be found in [41].
The complete time series analyzed in this work covers the time period between 1 January 2003 and 31 December 2017, with temporal resolution ∆t = 60 s, and with measurements recorded at 10 m and 38 m above ground level (AGL).
For this study, the data set was divided into 24 sub-intervals, each of the same duration, W L ≈ 90 days (W L ≈ 8 × 10 6 s), in order to avoid possible seasonal effects [42,43]. Figure 1 shows three samples of wind speed for three different sub-intervals. The autocorrelation function C( ) ≡ u(t)u(t + ) (where is a variable temporal lag) estimated for each 90-day interval, and shown in Figure 2 (left panel), provides a typical large eddy-turnover time t 0 ≈ 1 day, exposing the strong influence of the up-slope/down-slope (daily) cycle on signal correlation. This cycle is a low-level thermodynamic response to the radiative heating cycle [44]. During the day, solar heating over land surface increases the lower-tropospheric temperature (and moisture) and hence instability, leading to development of convective effects in the mixed boundary layer. As a result the variations in mixing ratios in the frontal mixing zone are governed by turbulence. Information about t 0 can be extracted from the analysis of the scale dependent local kurtosis, by dividing the sample into smaller windows of length T w < W L (n points), and then selecting the optimal block width T w for which the average kurtosis K(T w ) over all windows is the closest to that of the normal probability density function (PDF), i.e., K(T w ) = 3 ( Figure 2, right panel). The biased moment estimator for kurtosis can be defined as [45][46][47]: For shorter time windows, T w < 0.01 days (roughly 15 min), a value close to K(T w ) = 1.2 has been found, since for short time periods the value of u(t) is almost constant, characterized by less pronounced tails, with respect to a normal distribution. On the other hand, K(T w ≡ W L ) > 3 coincides with the kurtosis of the whole data set u(t). Finally, for an average window length T w ≈ 1 day, K(T w ) = 3 was found, the effects of mixing act on timescales shorter than 1 day, and at this scale the system completely loses all correlations and become Gaussian. This temporal scale could represent the upper limit of the range where scaling behavior holds.

Structure Function Analysis and Hurst Exponent Estimation for Mesoscale Wind Speed
The structure function analysis (SF) is the principal method used to extract the scaling properties of a physical phenomenon [48]: The relation between ζ(n) and the Hurst exponent H is well known: in classical Kolmogorov theory, in the absence of intermittent corrections, ζ(n) = nH. The scaling exponent S 1 , therefore, coincides with the Hurst exponent itself, ζ(1) = H, while the scaling exponent, S 2 , is directly linked to the power spectral density (PSD) as ζ(2) = β − 1 (β is the power-law exponent of the power spectrum).
H is a measure of the long-term memory of a time series. This represents a key parameter useful in order to classify the time series in terms of whether it is a random, a persistent, or an anti-persistent process. In other words, it is a description of how the "past" increments have on "future" ones. Values in the range H ∈ (0.5, 1] indicate a persistent (correlated) process, while values H ∈ [0, 0.5) are associated with anti-persistent (anti-correlated) processes. H = 0.5 indicates a completely uncorrelated process (e.g., a random walk). H is an important parameter in the description of fluid turbulence since it is related to the auto-correlation function of the velocity time series: the rate at which C( ) decreases with lag . Moreover, for self-similar (monofractal) processes, H is directly related to fractal dimension, D, by the relation: [49]. When H = 0, the average fluctuations δu exhibit a scale dependence.  This result clearly differs from other experimental estimates [37,50], where ζ(1) ≡ H ≈ 0.37. The last panel of Figure 3 shows the comparison among the exponents ζ(n) for MLO and other experimental estimates taken from literature [10,23,38]. The discrepancy in the exponents becomes more and more marked as the moments n increase, and it can be observed for all intervals (Figure 3, last panel). The averaged scaling exponents, and their standard deviations, are: ζ(1) ≈ 0.24 ± 0.01, ζ(2) ≈ 0.46 ± 0.02, ζ(3) ≈ 0.64 ± 0.04, ζ(4) ≈ 0.80 ± 0.06.
The dashed line in the last panel of Figure 3 represents a linear dependence of the scaling exponents on the order n of the form ζ(n) = nH, with H = 1/3. This relation is relative to the classical Kolmogorov K41 prediction n/3 [48], and is reported only here as a reference.

Hurst Exponent Estimation from the Empirical Mode Decomposition
Alternatively, H can be obtained using Empirical Mode Decomposition (EMD) [30]. EMD was developed to process and analyze the evolution of non-stationary data in order to obtain a set of basis functions derived (empirically) from the signal itself [51][52][53][54][55][56][57][58][59]. Since EMD analysis is adaptive (in contrast to traditional Fourier decomposition methods where the basis functions are fixed) and not restricted to stationary data [30,60], the data set may be analyzed without introducing spurious harmonics or artefacts near sharp data transitions, which could appear when using classical Fourier filtering or high-order moments analysis. Indeed, EMD allows local information to be extracted through the instantaneous frequencies which cannot be captured by fixed-frequency methods (such as Fourier or Wavelets). The main consequence is that the frequency is not widely spread (as for Wavelets), with a much better frequency definition and smaller amplitude variation induced frequency modulation [30,61]. Within EMD framework, the data are decomposed into a finite number k (approximately k ≤ log 2 (N) modes where N is the number of data sample) of elementary oscillating basis functions φ j (t), known as intrinsic mode functions (IMFs), characterized by an increasing time scale τ and a residual r k (t), which describes the mean trend, if one exists, given by The decomposition includes two stages: first, the local extrema of u(t) are identified and subsequently connected through cubic spline interpolation. Once connected, the envelopes of local maxima and minima are obtained. Second, the mean m 1 (t) is calculated between the two envelope functions, then subtracted from the original data, is an IMF only if it satisfies the following criteria: (i) the number of local extrema and zero-crossings does not differ by more than 1, and (ii) at any point t, the mean value of the extrema envelopes (e max (t) and e min (t)) is zero. When h 1 (t) does not meet the criteria above, the sifting procedure is repeated using h 1 (t) as the new raw data series, and h 11 (t) = h 1 (t) − m 11 (t) is generated, where m 11 (t) is the mean of the envelopes. The sifting procedure is repeated m times until h 1m (t) satisfies the criteria above. A general rule to stop the sifting is introduced by using a standard deviation σ, evaluated from two consecutive steps: The iterative process stops when σ is smaller than a threshold value typically of the order of [30,62]. Another widely used criterion is based on 3 thresholds α, θ 1 and θ 2 , in order to guarantee globally small fluctuations in the mean while taking into account locally large variations [63]. The sifting procedure is iterated until the parameter represents the mode amplitude) for a prescribed fraction of the total duration (1 − α), and σ(t) < θ 2 for the remaining fraction. The values for the three parameters used in this work are: α = 0.05, θ 1 = 0.05, and θ 2 = 10θ 1 , in accordance with their typical values [63]. The IMFs obtained from the decomposition of the sample April-June 2011 are reported in Figure 4. When the decomposition is applied on dataset possessing certain feature, such as noise time series, fractional Gaussian noise, turbulent time series, fractal time series, random walks, EMD acts intrinsically as a dyadic filter bank [31,[64][65][66]. Each IMF captures a narrow spectral band in frequency space [15,16,39,67] and their superposition behaves as Figure 4). By comparison with the Fourier PSD, each IMF can be interpreted according to its characteristic time scale. In particular, as visible in Figure 4, IMFs φ j |j ∈ [5, 12] capture a power-law dynamic, which seems to mimic a turbulent cascade in a fashion similar to the inertial sub-range [21,48], which develops over different time scales (roughly from tens of minutes up to 1 day).
Ideally, the ABL above MLO would be stratified, characterized by 4 principal zones: sub-cloud, fog (cloud), transitional (inversion), and free atmosphere [68]. Experimental observations [69] and numerical studies have shown that ABL at MLO can present both stable (during nighttime, with Richardson number Ri g < 0.2) and unstable (during daytime, Ri g > 0.2) conditions [70], and in case of strongly stratified (very stable) conditions (nighttime) ABL can be catheterized by a refractive index structure parameter of the order of C 2 n < 10 −14 [70]. However, on the mesoscale, due to the continuous thermal forcing from the surface (shallow convection), the influence of the descending limb of the Hadley cell, which leads to the formation of the trade wind inversion, and the interplay between synoptic-scale (North Pacific anticyclone) and mesoscale (thermal up-slope/down-slope winds) circulation systems, the heights and strengths of these atmospheric layers, are subject to strong fluctuations, which may leads to continuous mixing [68,70]. A possible explanation for the power-law behavior in the PSD (Figure 4) is that during this mixing phase, the large scale eddies under incessant motions undergo a fragmentation process, responsible for the formation of smaller scales. In the classical, self-similar Kolmogorov theory [21,48], the spectral slope β is linked to the second order SF via the relation (2). However, a strong discrepancy (Figure 3,4) is seen when comparing the two slopes: β − 1 and ζ(2), respectively.

Hurst Exponent Estimation from IMF Scaling Relations
Once the EMD decomposition has been performed, and the IMFs have been identified, an estimate of the Hurst exponent H can be obtained by relating the variance of the j-th IMF V j with its characteristic period T j [71], via the relation: The mean period of each IMF mode is estimated by calculating the local extrema points (N j,max , N j,min ) and zero-crossing (N j,0 ) [65,72]: where L is the temporal length of the data and j represents the IMF index.
Modes j ∈ [4,10], (approximately from T ≈ 28 min to T ≈ 1 day) have the same mean period for each sample (Figure 5, log-linear plot), and follow an exponential law, i.e.: in which α ∈ [10 −3 , 3 × 10 −4 ] and γ ∈ [1.9, 2.1] are obtained via a least squares fitting algorithm. The uncertainty associated with T j has been evaluated as the standard deviation of the inverse instantaneous frequency of the j-th IMF: and φ j (t) is the Hilbert transform of the j-th IMF. The first 3 modes do not follow the predicted scaling.
The value of γ close to 2 indicates the quasi dyadic filter bank property of the EMD algorithm for these series, as found in other works [64,65,67]. The average value of H, extracted from the 24 sections of the original data set, is shown in the left panel of Figure 6. The scaling relation 5 holds from T j ≈ 25 min up to the large eddy turnover time T j ≡ t 0 ≈ 1 day, in agreement with the the definition of the mesoscale range, and far beyond the inertial-range of turbulence. High-frequency IMFs (1 ≤ T j < 25 min, Figure 6 left panel) create a curvature that does not follow the scaling relation between V j and T j . Dashed and dotted lines are the average V j value and the 95% confidence bounds, respectively. The right-hand panel of Figure 6 shows the value of H for all the 90 day sections, from 2010 to 2013, with the associated average value and the 95% confidence bounds. These values are in good agreement with the classical estimation of the Hurst exponent reported in the literature [10,23]. Interestingly, similar values have also been reported in the inertial sub-range [37,[73][74][75].

Interoccurrence Times Statistics in Mesoscale Wind
A further method used to obtain information on the scaling properties of a stochastic processes, is theIOT τ analysis [33][34][35][36]. IOTs are a measure of the time between the occurrence of two or more subsequent extreme events in the data that exceed a fixed threshold Q (Figure 7, upper panel). For every Q, an average τ and standard deviation σ τ are defined. By increasing Q, τ and σ τ become larger. The higher the Q, the rarer or more extreme are the events. Furthermore, there is a one-to-one correspondence between Q and the τ , σ τ values (Figure 7, lower panels). In addition, it is known that if a correlation exists in the data, then the IOTs sequence is also correlated [76][77][78][79][80][81].
The behavior of τ and σ τ can be described by an exponential law of the form τ = τ 0 exp (α 1 Q) σ τ = σ 0 exp (α 2 Q), where τ 0 and σ 0 are the values obtained at Q = 0. Figure 7 (lower panels) shows the two exponential laws for τ and σ τ obtained from the average values of the two parameters α 1,2 ; all curves have a value of α 1 ∈ [1.05, 1.21] and α 2 ∈ [1.1, 1.3]. In terms of τ the values of T j obtained through the EMD lies in the interval Q ∈ [0.1, 3]. However, in certain samples the values of τ for thresholds Q > 2.8 seem to increase faster (Figure 7 lower panels). For this reason, the upper limit for the IOTs analysis was fixed at Q = 2.5, for a total of N Q ≈ 250 IOTs measured.
Empirical studies have shown that the PDFs of the IOTs can be described by a universal Tsallis q-exponential function, if a long-term memory exists in the data (i.e., power-law decay of the autocorrelation function) [34][35][36]82] of the form: where α is a normalization factor and q is a measure of the deviation from an exponential distribution (and q > 1 indicates a long-tailed distribution). The limit of validity for the parameter q lies in the range 1 ≤ q ≤ 2. The upper limit arises from the normalization condition of P(τ) to the unit area and the requirement that the normalization constant, α, is positive [82,83]. The PDFs of the IOT, P(τ), from the normalized wind speed data have been evaluated for a number of bins, N bins (N bins ∈ [10 ÷ 20]), logarithmically spaced between the minimum and maximum values of τ, of the width A bins . As it is the PDFs that are calculated, the number of bins and their width are irrelevant to the analysis. The sampling uncertainty on each bin of P(τ) is calculated from the statistical (Poisson) error on the histogram Err[N(x)] = N(x), which, following the error in the probability density becomes: Err[P(τ)] = P(τ)/(A bins N Q ). For threshold values Q ≥ 2.5, the statistics is too poor (N Q < 100) to permit accurate analysis, since the relative statistical error becomes too large. In order to compare the different distributions of the IOT, the PDFs of the normalized IOTs P(τ/σ τ ) were evaluated. This procedure ensures that for every Q the values of τ/σ τ always lies in the same interval. The first two panels of Figure 8 show the PDF P(τ/σ τ ) (symbols), obtained from the 90 days data sets, for the years 2012 and 2014, respectively, at various thresholds Q. The relative q-exponential fit (dashed line) obtained from Equation (8), at the maximum threshold Q reported in legend, is also shown. All the PDFs collapse on to the same q-exponential distribution characterized by an average q ≈ 1.65 ± 0.03, in a range of thresholds Q ∈ [0.3, 2.3], which, in terms of IMF periods, are T j ≈ 30 min up to T j ≈ 0.5 days, demonstrating the universality of the process as reported in [35,36]. This q-value is significant since a similar value has been reported in a different context, namely velocity fluctuations in fully-developed turbulence experiments for the inertial sub-range [35], where q ≈ 1.62 was obtained from the IOT statistics. The IOTs analysis shows that the mesoscale atmospheric dynamics, could share similar features than those observed in laboratory experiments or inertial sub-range turbulence. The q-exponential fitting was performed by minimizing the χ 2 statistics by varying the model parameters of Equation (8): q within the closed interval q ∈ [1 ÷ 2], while the other parameters in the open interval α, β ∈ (0 ÷ ∞).
To check the dependence of the correlated structure in the data set, all the 90 days dataset were randomly shuffled for a large number of trials, keeping the random seed fixed. The last panel of Figure 8 reports the PDF of the shuffled data set in a log linear plot. The random permutation destroys the correlation in the data and, moreover, the distribution becomes independent of the threshold Q and collapses onto the same exponential shape: P(τ) = σ −1 τ exp [−τ/σ τ ], demonstrating their uncorrelated nature.
The value of q in Equation (8) is, therefore, linked to the correlation properties of the data. Following [35], this value can be linked to the Hurst exponent as: q = (3 − H)/(2 − H), where H = 0.46 ± 0.06 (average and standard deviation obtained from different q-exponential fits), which, within the error bar, is compatible with EMD results. IOTs scaling leads to an estimate of H closer (within the error bars) to the values reported in literature for the inertial sub-range [37,[73][74][75], compared to the scaling of the structure function S 1 . The slightly higher value obtained via IOTs measurement could be due to the finite size of the sample. In other words, in order to minimize the error in the estimate, a very large number of IOTs (and, consequently, a longer dataset) is required.

Scaling of High-Order Moments: Intermittency and Arbitrary Order Hilbert Spectral Analysis
To correctly obtain scaling information for wind speed data sets and reduce the impact of the up-slope/down-slope cycle, arbitrary-order Hilbert Spectral Analysis (HSA) [16,72,84] was applied to the wind speed samples. HSA represents an extension of the classical Hilbert-Huang transform [31] designed to characterize scale-invariant properties of a time series. Starting from the EMD, once the IMFs have been obtained, the subsequent step consists of evaluating the Hilbert transform of each mode where p is the Cauchy principal value and φ j (t) is the j-th IMF. The combination of φ j (t) and φ j (t) defines the so called analytical signal Z = φ j + iφ j = A j (t)e iθ(t) , where A j (t) is the time-dependent amplitude modulation and θ(t) is the phase of the mode oscillation [85].
For each mode, the Hilbert spectrum, defined as H(ω, t) = A 2 (ω, t) (where ω = dθ/dt is the instantaneous frequency), provides energy information in the time-frequency domain [67]. A marginal integration of H(ω, t) provides the Hilbert marginal spectrum h(ω) = T −1 T 0 H(ω, t)dt, defined as the energy density at frequency ω [30,86]. In addition, from the Hilbert spectrum, a joint probability density function P(ω, A) can be extracted, using the instantaneous frequency ω j and the amplitude A j of the j-th IMF. This allows the Hilbert marginal spectrum h(ω) to be written as which corresponds to a second order statistical moment [67]. Equation (10) can be generalized to the arbitrary order q ≥ 0 by defining the ω-dependent qth-order statistical moments In particular h(ω) ≡ L 2 is similar to the Fourier spectral energy density, and can be interpreted as the energy associated to the frequency ω [67]. However, it should be pointed out that the definition of frequency in h(ω) is different from the definition in the Fourier framework. The interpretation of the Hilbert marginal spectrum should be given with more caution [30,86,87]. Figure 9 shows three examples of L n (ω) up to the 5-th order, obtained from Equation (11) using the samples January-March 2010, April-June 2012 and July-September 2013. The resulting L n (ω) demonstrates power-law behavior L n (ν) ∝ ω −β n for all n, in the range of frequencies ω ∈ [1.6 × 10 −1 , 4 × 10 1 ] days −1 (approximately from 30 min to 6 days), and the slope β 2 ≈ 1.67 ± 0.02 is compatible with the slope of the Fourier spectrum ( Figure 4). This range is wider than the previous estimate ( Figure 6), due to the local nature of EMD and HSA, the strong daily modulation, as well as the possible "non-stationarity" due to ramp-cliff structures, can be constrained, isolating the properties of the cascade from the possible effects of the larger scale forcing and residual structures [15,16,39]. The last panel of Figure 9 shows the evolution of the slope β 2 as a function of time, with the associated average value and standard deviation: β 2 = 1.64 ± 0.03, compatible with the β ≈ 5/3 characteristic of the inertial sub-range turbulence [21,48].
At higher frequencies, for ω ≈ 50 − 60 days −1 (T ≥ 30 min) a spectral break is observed, characterized by an abrupt slowing of the spectral slope, which tends to reach a saturation at higher order n.
As discussed in Section 3, the relation between the spectral slope and the scaling exponent of the second order SF is E(ω) ∼ ω −β → β − 1 = ζ (2). By extending this relationship to any arbitrary order q, a family of generalized scaling exponents ξ(q) can be introduced through the generalized Hilbert spectra [15,16] as ξ(q) ≡ β q − 1. The exponents ξ(q) are the Hilbert analogues of the standard scaling exponents ζ(q) obtained through the structure functions or through Extended Self-Similarity (ESS) [10,23,38,88]. Equation (11), therefore, is an alternative to the structure function scaling exponents to quantitatively estimate the level of intermittency in the turbulent cascade [21], with the advantage of constraining the effects of noise and large-scale structures [15]. Examples of scaling exponents ξ(n) obtained via least square fits of the generalized Hilbert spectra are shown in Figure 10. The same figure also shows different values of the exponents taken from literature, for laboratory experiments, and atmospheric turbulence. hydrodynamic turbulence using extended self-similarity (ESS) (squares) [38], atmospheric boundary layer in the mesoscale regime (diamonds and stars) [10,23], and other estimates in hydrodynamic turbulence [35,89]; the dashed line represents the theoretical expectation n/3, as estimated from dimensional analysis in the absence of intermittency [48]. From the results, it is easily observed that the departure from a linear (monofractal) scaling is captured in the mesoscale wind data. Moreover, the values of ξ(n) are similar to the exponents obtained in other experiments in comparable range of scales (ten min up to six h).
Moment n = 4 deviates from the scaling of homogeneous and isotropic flow ( Figure 10, left and central panel), and converges to the values reported in Reference [10], without the necessity of employing the ESS procedure. In the above reference, the generalized scaling exponent ζ 4,3 = ζ 4 ζ −1 3 is evaluated in terms of ESS, and in case of complex orography. It is important to note that, in the presence of strong shear, the applicability of ESS seems to be somewhat limited [91,92].
Higher orders are affected by the finite sample size. In other words, the sample size is not sufficient to ensure the convergence of the exponents; heuristically, the convergence of the results is guaranteed up to order n max ≈ log 10 [W L /∆t] − 1 ≈ 4.
The last panel of Figure 10 shows the evolution of the intermittency parameter µ obtained from a simple log-normal cascade model [93,94]: As shown in the last panel of Figure 10, the intermittency effects are clearly present in each sample, but the values of the parameter µ varies from case to case. The figure also shows (dashed line), the reference for the inertial sub-range µ ≈ 0.02. The intermittency parameter obtained for the MLO wind speed, is slightly lower the isotropic case, and the characteristic average value and standard deviation is µ ≈ 0.0165 ± 0.0031. The value of µ is in good agreement with other estimates obtained from the classical multifractal spectrum [90] and from magnitude covariance analysis [7,8]. Intermittency effects could suggest the existence of a universal cascade mechanism associated with the energy transfer between synoptic motions and turbulent microscales in the atmospheric boundary layer. The values of µ have been obtained via least square fitting of the exponents ξ(n) determined from each sample; the orders n > 4 have been excluded in order to make allowance for the finite sample size effects.

Conclusions
The scaling properties of the mesoscale wind speed fluctuations at MLO have been analyzed with different methods, and all the results show empirical evidence that large scales phenomena possess a certain degree of universality and share a number of features and characteristics with fully developed turbulence in the inertial sub-range. The findings presented in this work, specifically the statistics for intermittency and cascade models, can be applied to atmospheric fields, such as, wind resource assessment, extreme event characterization, and short-term wind predictions. Moreover, most numerical weather prediction models, in their standard operational configuration, are unable to capture, or reproduce, the scaling behavior observed in the mesoscale regime, which is required to improve turbulence closure parameterizations [23], i.e., the power-law scaling of the PSD characterized by an exponent β 2 ≈ 1.64 ± 0.03, characteristic of a cascade dynamic.
The EMD framework, yielded a value of the Hurst exponent, H ≈ 0.32 ± 0.04, close to the classical value of the fully developed turbulence, over a temporal range of T j ∈ [0.02, 1] days (from 10-20 min up to 1 day), in agreement with the values reported in the literature [7,8,10]. For a scale of T j = 1 day the wind speed data are normally distributed and characterized by a local kurtosis K(T w ) = 3, exposing the convective effects within the mixed boundary layer [44].
By using the statistics from IOTs, it was found that above a fixed threshold Q ∈ [0.3, 2.3], the PDFs of the normalized IOTs are described by a Tsallis q-exponential function [82].
The range of timescales associated with the IOTs, at these threshold values, is compatible with those obtained via EMD (T j ∈ [0.02, 1]). The value of the parameter q at the MLO, q ≈ 1.65 ± 0.03, is in good agreement with the value found for the inertial sub-range [35]. The dependence of the correlated structures in the data set has been shown by simply randomly shuffling the data in the time series. Randomization destroys the correlation in the data, and the effect is observed in the PDFs of the normalized IOTs, which become independent of the threshold Q and collapse onto the same exponential shape, demonstrating their complete lack of correlation. The Hurst exponent extracted from the q-exponential fit has a slightly higher value than that found using EMD, or from estimates in the literature. It is the finite sample size in this case that influences the estimate of the Hurst exponent.
The scaling of the higher order moments were studied with the HSA. At the mesoscale, departure from linear (monofractal) scaling was observed. Intermittency effects were clearly present in each sample, however the values of the intermittency parameter µ varied from case to case. The average value, with the associated standard deviation, was found to be of the order of µ ≈ 0.0165 ± 0.0031, in good agreement with estimates obtained from the classical multifractal spectrum [90], and from magnitude covariance analysis [7,8]. The scaling exponents ξ(n) obtained thorough HSA are similar to the exponents obtained in other experiments through ESS, up to the fourth order n = 4.
These results suggest the existence of a universal cascade mechanism associated with the energy transfer between synoptic motions and turbulent microscales in the atmospheric boundary layer. Funding: This research was funded by the European Commission-H2020, the ERA-PLANET programme (www.era-planet.eu) (Contract.no. 689443) within the IGOSP project (www.igosp.eu).

Conflicts of Interest:
The authors declare no conflict of interest.