1. Introduction
The distribution of turbulent kinetic energy (TKE) across wavenumbers remains as a central topic in fluid dynamics. The streamwise velocity spectrum of turbulent flows exhibits distinct canonical regimes: the energy-containing range at low wavenumbers, dominated by large-scale eddies responsible for generating TKE; the inertial subrange at intermediate wavenumbers, governed by Kolmogorov’s
$-5/3$
scaling law; and the dissipation range at sufficiently high wavenumbers, where viscosity fully suppresses TKE (Kolmogorov Reference Kolmogorov1941). In the classical Kolmogorov picture, the energy cascades without loss between the energy-containing and dissipative scales. However, this canonical picture is partial in wall-bounded flows and geophysical turbulence, where spectral scaling at large scales deviates significantly. At low wavenumbers, the dynamics of turbulent boundary layers is shaped by anisotropy and long-range dependencies, making it analytically intractable. In particular, a
$k^{-1}$
regime has been consistently observed, and debated over the past few decades, within the energy-containing range of the velocity spectrum sometimes referred to as the shear-production subrange (Tchen et al. Reference Tchen, Larsen, Pecseli and Mikkelsen1985), where
$k$
is the streamwise wavenumber. Empirical evidence for this scaling has been observed in the inertial sublayer of laboratory turbulent boundary layers in pipes and wind tunnels (Nickels et al. Reference Nickels, Marusic, Hafez and Chong2005) and in the atmospheric surface layer, typically defined as the lowest 10 % of the atmospheric boundary layer (ABL) (Katul & Chu Reference Katul and Chu1998).
Theoretical foundations for the
$k^{-1}$
scaling were first introduced by Tchen (Reference Tchen1953) based on a spectral budget. Subsequent developments connected this behaviour to Townsend’s attached eddy hypothesis (AEH) (Townsend Reference Townsend1976; Marusic & Monty Reference Marusic and Monty2019), which describes the energy-containing motions in the logarithmic region as arising from a hierarchy of self-similar eddies attached to the wall. This concept was formalised by Perry, Henbest & Chong (Reference Perry, Henbest and Chong1986), who used dimensional arguments to predict a
$k^{-1}$
scaling in the intermediate-wavenumber range, arising from the overlap between the inner (wall-normal distance-based) and outer (boundary-layer height-based) scaling ranges at sufficiently high Reynolds numbers.
Alternative explanations have been proposed as well. Nikora (Reference Nikora1999) introduced a conceptual model based on superimposed eddy cascades originating at various wall-normal positions. Katul, Porporato & Nikora (Reference Katul, Porporato and Nikora2012) synthesised this view into a unified framework that integrates Tchen’s spectral analysis, Nikora’s eddy cascade picture and Heisenberg’s eddy viscosity concept. Fundamental arguments by Katul and colleagues suggest that the
$k^{-1}$
scaling can be rationalised in the limit of infinite Reynolds number and under certain boundary conditions (Katul et al. Reference Katul, Porporato and Nikora2012). This scaling is attributed to the balance between energy production by large eddies and energy transfer to smaller scales, with the wall acting as a source of anisotropy and organisation. Experimental data from ABL studies, including measurements over flat terrain and within canopy sublayers, consistently support the presence of the
$k^{-1}$
regime (Katul et al. Reference Katul, Chu, Parlange, Albertson and Ortenburger1995; Cava & Katul Reference Cava and Katul2012).
Perry et al. (Reference Perry, Henbest and Chong1986) also linked the
$k^{-1}$
region to the logarithmic variation of streamwise turbulence intensity, although recent results by Hwang, Hutchins & Marusic (Reference Hwang, Hutchins and Marusic2022) suggest that such a logarithmic trend can emerge even without a fully developed
$k^{-1}$
range, particularly at moderate Reynolds numbers. Furthermore, experimental studies by Rosenberg et al. (Reference Rosenberg, Hultmark, Vallikivi, Bailey and Smits2013) and Deshpande, Monty & Marusic (Reference Deshpande, Monty and Marusic2021) have shown that non-self-similar contributions may mask the appearance of the
$k^{-1}$
region in the energy spectrum. Davidson, Nickels & Krogstad (Reference Davidson, Nickels and Krogstad2006) and Davidson & Krogstad (Reference Davidson and Krogstad2009) argued that a hierarchy of space-filling eddies with kinetic energy scaling as
$u_*^2$
necessarily leads to a
$k^{-1}$
spectral band, where
$u_*$
is the friction velocity. However, their formulation assumes scale-wise independence, a prescribed Gaussian correlation kernel for eddies of a given size and does not account for long-range dependencies. In the physical space, Davidson et al. (Reference Davidson, Nickels and Krogstad2006) and Davidson & Krogstad (Reference Davidson and Krogstad2009) associated the
$k^{-1}$
scaling with a logarithmic form of the second-order structure function, consistent with field observations in ABL turbulence (Chamecki et al. Reference Chamecki, Dias, Salesky and Pan2017; Ghannam et al. Reference Ghannam, Katul, Bou-Zeid, Gerken and Chamecki2018).
We show that a specific correlation structure is not a prerequisite for the emergence of
$k^{-1}$
scaling. Instead, we employ asymptotic theorems from random field theory, which rigorously establish a direct link between covariance functions and their spectral counterparts. Our formulation mixes clusters of correlated eddies, with each cluster sharing a common integral length scale (cluster size), rather than isolated, uncorrelated eddies. The only statistical premise is a log-uniform distribution of cluster sizes, a continuous analogue of Townsend’s equal-population condition; no wall attachment, Gaussianity or inertial-range
$k^{-5/3}$
scaling is imposed. We verify the theory with riverine and coastal ABL data where Townsend’s AEH is not strictly applicable.
By identifying the minimal statistical requirement, a log-uniform distribution of energy-containing scales, for the appearance of the
$k^{-1}$
law, our theory offers a more general approach that can inform sub-grid turbulence models and aid in interpreting spectra from complex geophysical and environmental flows, such as those found in tidal channels (Thomson et al. Reference Thomson, Polagye, Durgesh and Richmond2012). The proposed approach accounts for long-range correlations and does not rely on Kolmogorov scaling in the inertial subrange as a foundational premise. We further demonstrate its practical utility both within and beyond the surface layer through analysis of experimental riverine and wind datasets, thereby relaxing the assumptions of Townsend’s AEH. Section 2 presents the stochastic formulation and analytic derivation, § 3 validates the theory against tidal-channel acoustic Doppler current profiler/acoustic Doppler velocimeter measurements, and a multi-level sonic-anemometer array on the ABL at the Yucatán coast, spanning a large range of Reynolds numbers and § 4 summarises key findings.
2. Proposed formulation
The central premise of our formulation is that there exists a hierarchy of eddy clusters characterised by different integral length scales and that the statistical superposition of these clusters results in the observed streamwise velocity fluctuations. Let the integral length scale corresponding to a cluster be identified by a representative parameter
$c$
and that there is a continuous range of clusters with the parameter
$c$
varying from
$c_1$
to
$c_2$
(see figure 1). The resulting streamwise (longitudinal) velocity fluctuation,
$u$
, at a point
$\boldsymbol{x}_1=(x_1,y_1,z_1)$
can then be expressed as
\begin{align} u(\boldsymbol{x}_1,\omega ) = \int _{c_1}^{c_2} \mathfrak{u}(\boldsymbol{x}_1,c,\omega )\,{\textrm{d}}c, \end{align}
where
$\mathfrak{u}(\boldsymbol{x}_1,c,\omega )$
represents the contribution from eddy clusters of scale
$c$
, and
$\omega$
denotes a specific realisation in the probability space. The streamwise covariance function
$\varPsi (\boldsymbol{x}_1,\boldsymbol{x}_2)=\langle u(\boldsymbol{x}_1,\omega )u(\boldsymbol{x}_2,\omega )\rangle$
, with
$\boldsymbol{x}_2=(x_2,y_1,z_1)$
can be obtained as follows:
\begin{align} \varPsi (\boldsymbol{x}_1,\boldsymbol{x}_2) &= \left \langle \int _{c_1}^{c_2} \mathfrak{u}(\boldsymbol{x}_1,c,\omega )\,{\textrm{d}}c \int _{c_1}^{c_2} \mathfrak{u}(\boldsymbol{x}_2,c,\omega )\,{\textrm{d}}c \right \rangle \nonumber \\ &= \int _{c_1}^{c_2} \int _{c_1}^{c_2} \langle \mathfrak{u}(\boldsymbol{x}_1,c,\omega ) \mathfrak{u}(\boldsymbol{x}_2,c^{\prime },\omega ) \rangle \,{\textrm{d}}c\,{\textrm{d}}c^{\prime }. \end{align}
Here,
$\langle \boldsymbol{\cdot }\rangle$
represents an ensemble average. Now, we assume that eddy clusters characterised by different values of
$c$
are uncorrelated, leading to
\begin{align} \varPsi (\boldsymbol{x}_1,\boldsymbol{x}_2) & = \int _{c_1}^{c_2} \int _{c_1}^{c_2} F(\boldsymbol{x}_1,\boldsymbol{x}_2,c,c^{\prime })\delta (c-c^{\prime })\,{\textrm{d}}c\,{\textrm{d}}c^{\prime } = \int _{c_1}^{c_2} F(\boldsymbol{x}_1,\boldsymbol{x}_2,c)\,{\textrm{d}}c. \end{align}
Here,
$F(\boldsymbol{x}_1,\boldsymbol{x}_2,c)$
is the covariance associated with eddy clusters of scale
$c$
. This assumption differs from models that assume no correlation between eddies of different sizes (e.g. Davidson et al. (Reference Davidson, Nickels and Krogstad2006)).
Figure 1. Concept illustrating a superposition of eddy clusters with integral scale
$c$
, whose population is inversely proportional to
$c$
. Smaller clusters cascade at higher wavenumbers, while larger ones contribute to lower wavenumbers.
Let the streamwise velocity fluctuation be modelled as a wide-sense stationary and isotropic random field. Under statistical isotropy, the streamwise covariance function
$\varPsi (\boldsymbol{x}_1,\boldsymbol{x}_2)$
depends only on the distance
$x=\left \| \boldsymbol{x}_2-\boldsymbol{x}_1 \right \|=|x_2-x_1|\gt 0$
. For notational simplicity, we refer to
$\varPsi (\boldsymbol{x}_1,\boldsymbol{x}_2)$
as
$\varPsi (x)$
, with its spectral counterpart given by
$\widehat {\varPsi }(k)$
. Here,
$\widehat {\varPsi }(k)$
denotes the isotropic one-dimensional energy spectral density of streamwise velocity fluctuations. Under these restrictions,
$F(\boldsymbol{x}_1,\boldsymbol{x}_2,c)$
has to be of the form
$F(x,c)$
. An ideal choice of
$F(x,c)$
would result in
$\varPsi (x)$
such that
$\widehat {\varPsi }(k)$
exhibits appropriate scaling characteristics across the full wavenumber spectrum. This includes any long-memory-dependent scaling in the low-wavenumber regime,
$k^{-1}$
scaling in the shear-production subrange (if present), Kolmogorov’s scaling in the inertial subrange and the appropriate decay in the dissipation range. In this work, we propose a choice of
$F$
such that all these scaling regions are accurately captured, with the exception of the dissipative range.
By extending Kolmogorov scaling beyond the inertial subrange and disregarding dissipative effects, which is common in classical turbulence models (e.g. Von Kármán & Lin Reference Von Kármán and Lin1951), the inertial subrange scaling can be interpreted in terms of the fractality of
$u$
. The fractal behaviour is characterised by the fractal dimension
$D$
. At the same time, any long memory dependence of
$u$
, which governs the scaling in the low-wavenumber regime, is described by the Hurst parameter
$H\in (0,1)$
. For an isotropic random field,
$D$
and
$H$
can be estimated from the asymptotic behaviour of the covariance function. If
\begin{equation} 1 - \varPsi (x) \sim x^{\alpha }, \qquad x\to 0, \end{equation}
for some
$\alpha \in (0, 2]$
, then
$D = d + 1 - \alpha /2$
, where
$d$
is the spatial dimension of the random field. Since we characterise only one component of the velocity field in this work, we take
$d=1$
in this work, yielding
$1 \leq D \lt 2$
. On the other hand, if
\begin{equation} \varPsi (x) \sim x^{-\beta }, \qquad x\to \infty , \end{equation}
for some
$\beta \in (0, 1)$
, then
$H = 1 - \beta /2$
. Here,
$g\sim h$
indicates that the function
$g$
behaves like
$h$
under the specified asymptotic limit. The fractal dimension
$D$
determines the high-wavenumber asymptotic behaviour
$\widehat {\varPsi }(k) \sim k^{2D-5}$
. For isotropic turbulence,
$D = 5/3$
, recovers the well-known Kolmogorov
$k^{-5/3}$
scaling. Conversely, the Hurst parameter
$H$
governs the low-wavenumber asymptotic behaviour, yielding
$\widehat {\varPsi }(k) \sim k^{1-2H}, k \to 0$
. These asymptotic properties are rigorously derived from the Abelian and Tauberian theorems within a measure-theoretic framework, as detailed in § 2.8 of Stein (Reference Stein2012). Note that, in conventional random field models,
$D$
and
$H$
are linearly related through a self-similarity condition. In this work, we relax this assumption based on recent covariance models that treat the fractal dimension and the Hurst exponent as independent parameters (e.g. Jetti, Porcu & Ostoja-Starzewski (Reference Jetti, Porcu and Ostoja-Starzewski2023)).
We now propose
$F$
to be of the form
$\psi (x/c)/c$
, with
$\psi (0) = 1$
, and
$\psi (x/c)$
reproduces the same long- and short-range scaling properties as
$\varPsi (x)$
. However, its spectral counterpart
$\widehat {\psi }(k)$
does not exhibit the
$-1$
spectral scaling that may be present in
$\widehat {\varPsi }(k)$
. We later demonstrate that this choice of
$F$
accurately recovers the scaling behaviour across all wavenumber regions, except in the dissipative range. In this choice of
$F$
, the contribution of eddy clusters is weighted inversely with respect to their characteristic scale
$c$
, ensuring that the larger eddies appear less frequently than the smaller ones. This scaling is in line with Townsend’s AEH, which posits a hierarchy of self-similar eddies whose population is inversely proportional to their size (Townsend Reference Townsend1976). Moreover, this scaling naturally reflects the hierarchical organisation of turbulence and provides a way to capture the statistical structure of velocity fluctuations. Building on this premise, we propose that the covariance function
$\varPsi (x)$
is given by the following integral formulation:
\begin{equation} \varPsi (x) = \frac {\sigma ^2}{\log {(c_2/c_1)}}\int _{c_1}^{c_2}\frac {\psi (x/c)}{c} \,{\textrm{d}}c, \quad x \geq 0. \end{equation}
Note that
$\sigma ^2/\log {(c_2/c_1)}$
serves as a normalisation term to ensure that the variance of
$u$
is
$\sigma ^2$
. This formulation models turbulent flow as an ensemble of eddy clusters that span a hierarchy of scales, capturing the statistical structure of turbulence while allowing for variability in eddy sizes (see figure 1). The weighting function follows
$1/c\,{\textrm{d}}c = d(\log {c})$
, indicating that an alternative interpretation of the formulation involves eddy cluster sizes being distributed logarithmically.
This process establishes a mechanism through which the
$-1$
spectral scaling emerges in the velocity spectrum. Specifically, when the range of integral scales satisfies
$c_1 \neq c_2$
, a finite
$-1$
scaling region appears in the spectrum
$\widehat {\varPsi }(k)$
, as will be demonstrated next. In contrast, when
$c_1 = c_2$
, the formulation reduces to a single-scale covariance function,
$\varPsi (x) = \sigma ^2 \psi (x/c_1)$
, which does not exhibit
$-1$
scaling, as
$\psi (x/c)$
lacks this behaviour inherently. The range of length scales is bounded by
$c_1$
and
$c_2$
, flow-dependent parameters. This scale distribution governs the emergence of the
$-1$
spectral region, since eddy clusters associated with smaller
$c$
values contribute to higher wavenumbers, while those with larger
$c$
values populate the lower-wavenumber range (see figure 1). The bounds
$c_1$
and
$c_2$
may be associated with the characteristic scales of large-scale motions and very-large-scale motions (VLSMs), which have been shown to roughly demarcate the
$k^{-1}$
spectral regime (e.g. Kim & Adrian (Reference Kim and Adrian1999)). While recent studies suggest that VLSMs may not exhibit strict structural self-similarity with smaller-scale motions (e.g. Deshpande et al. (Reference Deshpande, Monty and Marusic2021)), this does not preclude the utility of scale-invariant assumptions as a leading-order approximation. Nonetheless, the extent to which such deviations influence model fidelity warrants further investigation in future work.
We assume that
$\psi (x/c)$
retains the same long-range and short-range scaling characteristics as
$\varPsi (x)$
. Specifically, the spectral counterpart
$\widehat {\psi }(k)$
follows the asymptotic behaviour
\begin{align} \widehat {\psi }(k) &\to C_{00} c (kc)^p, \quad k \to 0, \end{align}
\begin{align} \widehat {\psi }(k) &\to C_{10} c (kc)^q, \quad k \to \infty , \end{align}
where
$p=1-2H$
,
$q=2D-5$
,
$C_{00}$
and
$C_{10}$
are positive constants independent of
$c$
. In isotropic turbulence,
$C_{10} = C_K \epsilon$
, where
$C_K$
is the universal Kolmogorov constant and
$\epsilon$
is the mean energy dissipation rate. Physically, as
$k \to \infty$
,
$\widehat {\psi }(k)$
should transition to the dissipative scaling regime. However, the proposed formulation does not explicitly incorporate dissipation effects and, instead, it assumes that the Kolmogorov scaling extends beyond the inertial subrange, as noted earlier. A recently introduced parametric family of isotropic covariance functions,
$\psi (\boldsymbol{\cdot })$
, has been shown to capture both the long-range and short-range scaling characteristics required for turbulent flow characterisation (Jetti et al. Reference Jetti, Porcu and Ostoja-Starzewski2023, Reference Jetti, Cheng, Porcu, Chamorro and Ostoja-Starzewski2025b
). Interestingly, a recent study by Jetti et al. (Reference Jetti, Malyarenko and Ostoja-Starzewski2025a
) demonstrated that, under statistically isotropic conditions, the fractal and Hurst characteristics of the lateral velocity components are identical to those of the streamwise velocity component. This provides a natural extension of the approach proposed here to the modelling of other velocity components.
For analytical tractability, we approximate
$\widehat {\psi }(k)$
using a piecewise formulation
\begin{equation} \widehat {\psi }^{\textit{app}}(k) = {\begin{cases} C_{00} c (kc)^p, & k \leq k_T, \\ C_{10} c (kc)^q, & k \gt k_T. \end{cases}} \end{equation}
This approximation provides a practical means of capturing the dominant spectral trends while preserving the essential scaling behaviour at both extremes.
The transition wavenumber,
$k_T$
, where the two asymptotic forms of
$\widehat {\psi }(k)$
intersect, is given by
\begin{equation} k_T=\dfrac {1}{c}\left (\dfrac {C_{10}}{C_{00}}\right )^{1/(p-q)}\!. \end{equation}
Physically,
$k_T$
represents the wavenumber at which energy transitions between the two scaling regimes for a given eddy cluster characterised by
$c$
. As discussed earlier, larger values of
$c$
correspond to smaller
$k_T$
, indicating that larger-scale eddy clusters cascade at lower wavenumbers, while smaller
$c$
values lead to higher
$k_T$
, reflecting energy transfer at finer scales.
The approximate form of the spectrum
$\widehat {\varPsi }(k)$
, corresponding to the covariance function in (2.6), can be derived using the piecewise approximation in (2.9). For
$k \leq k_T(c_2)$
, spectra across
$c_1 \leq c \leq c_2$
exhibit low-wavenumber behaviour, leading to
\begin{equation} \widehat {\varPsi }^{\textit{app}}(k) =\frac {\sigma ^{2}}{\log (c_{2}/c_{1})} \int _{c_{1}}^{c_{2}} C_{00}\,c\,(kc)^{p}\,d(\log c) =\frac {\sigma ^{2}C_{00}k^{p}\big (c_{2}^{\,1+p}-c_{1}^{\,1+p}\big )} {(1+p)\,\log (c_{2}/c_{1})} \sim k^{p}. \end{equation}
Conversely, for
$k \geq k_T(c_1)$
, spectra across
$c_1 \leq c \leq c_2$
exhibit high-wavenumber behaviour, yielding
\begin{equation} \widehat {\varPsi }^{\textit{app}}(k) =\frac {\sigma ^{2}}{\log (c_{2}/c_{1})} \int _{c_{1}}^{c_{2}} C_{10}\,c\,(kc)^{q}\,d(\log c) =\frac {\sigma ^{2}C_{10}k^{q}\big (c_{2}^{\,q+1}-c_{1}^{\,q+1}\big )} {(q+1)\,\log (c_{2}/c_{1})} \;\sim \;k^{q}. \end{equation}
Within the transition region,
$k_T(c_2) \lt k \lt k_T(c_1)$
, we obtain
\begin{align} \widehat {\varPsi }^{\textit{app}}(k) &= \frac {\sigma ^{2}}{\log (c_{2}/c_{1})} \left( \int _{c_{1}}^{c_{k}} C_{00}\,c\,(k c)^{p}\,d(\log c) +\int _{c_{k}}^{c_{2}} C_{10}\,c\,(k c)^{q}\,d(\log c) \! \right)\nonumber \\ &= \frac {\sigma ^{2}}{\log (c_{2}/c_{1})} \left [ \frac {C_{00}}{1+p} \left (\frac {C_{10}}{C_{00}}\right )^{\!\frac {1+p}{\,p-q}}\,k^{-1} -\frac {C_{10}}{q+1} \left (\frac {C_{10}}{C_{00}}\right )^{\!\frac {q+1}{\,p-q}}\,k^{-1} \right . \nonumber\\ & \quad \left . -\frac {C_{00}}{1+p}\,c_{1}^{\,1+p}\,k^{\,p} -\frac {C_{10}}{q+1}\,c_{2}^{\,q+1}\,k^{\,q} \right]\!. \end{align}
Within
$k_T(c_2) \lt k \lt k_T(c_1)$
, the dominant term is the coefficient of
$k^{-1}$
. Since
$\widehat {\varPsi }^{\textit{app}}(k) \geq 0$
, it follows that
$\widehat {\varPsi }^{\textit{app}}(k) \sim \mathcal{O}(k^{-1})$
in this range. For isotropic turbulence, Kolmogorov’s scaling fixes
$D=5/3$
. The Hurst parameter
$H$
varies between
$0.5$
and
$1$
, where higher values indicate greater persistence (long-range dependence). As
$H$
increases, the second term becomes significant and, in the limit
$H \to 1$
, it also exhibits the
$k^{-1}$
scaling, reinforcing the dominant behaviour. This result demonstrates the emergence of a distinct
$-1$
spectral scaling regime. Figure 2(a) illustrates
$\widehat {\varPsi }^{\textit{app}}(k)$
for a representative set of parameters (
$D=5/3$
,
$H=0.75$
), highlighting the
$-1$
spectral scaling between
$k_1 := k_T(c_1)$
and
$k_2 := k_T(c_2)$
.
The real-space analogue of the
$-1$
spectral scaling is the logarithmic dependence of the structure function (Davidson et al. Reference Davidson, Nickels and Krogstad2006). The second-order structure function,
$S_2(x)$
, is defined as
\begin{equation} S_2(x) = \big \langle (u(\boldsymbol{x}_1) - u(\boldsymbol{x}_2))^2 \big \rangle = 2(\sigma ^2 - \varPsi (x)), \end{equation}
where
$x = \left \| \boldsymbol{x}_2-\boldsymbol{x}_1 \right \|$
. The proposed formulation in (2.6) inherently produces this logarithmic dependence, as illustrated in figure 2(b). The detailed derivations for the logarithmic dependence follow similar steps to those used in the derivation of
$k^{-1}$
scaling. Therefore, these derivations are provided in Appendix A. This result establishes that the emergence of
$-1$
spectral scaling directly results from the logarithmic variation of
$\varPsi (x)$
, as prescribed by (2.6). The interpretation of the
$-1$
spectral scaling as an emergent property of eddy interactions across multiple scales aligns with the perspective of Nikora (Reference Nikora1999). However, our formulation demonstrates that Townsend’s AEH, a fundamental assumption in Nikora (Reference Nikora1999), can be relaxed. Specifically, the occurrence of
$-1$
spectral scaling does not necessitate the strict presence of attached eddies but instead depends on the statistical distribution of eddy sizes within a cluster, as described by (2.6). Since the proposed model depends only on the distribution of integral-like scales and not on wall attachment or mean shear, it may also apply to flows such as mixing layers and wakes at high Reynolds numbers, where broad and self-similar eddy populations can produce an intermediate
$k^{-1}$
range.
This relaxation of AEH is particularly relevant for modelling tidal and riverine flows, where we show that the
$-1$
velocity spectrum scaling can persist even outside the surface layer, extending the applicability of the proposed framework to broader geophysical settings.
Figure 2. For a set of parameters (
$D=5/3$
and
$H=0.75$
) (a) the function
$\widehat {\varPsi }^{\textit{app}}(k)$
with
$-1$
spectral scaling between
$k_1$
and
$k_2$
. (b) Structure function
$2(\sigma ^2-{\varPsi }^{\textit{app}}(x))$
with logarithmic dependence between
$x_1\simeq 1/k_1$
and
$x_2\simeq 1/k_2$
. Here,
$D=5/3$
corresponds to the Kolmogorov’s
$k^{-5/3}$
scaling at high wavenumbers
$(k\geq k_1)$
, while
$H=0.75$
indicates long-range dependence and yields the
$k^{-0.5}$
scaling at low wavenumbers
$(k\leq k_2)$
.
3. Experimental validation
Next, we analyse a tidal current dataset to demonstrate the utility of the proposed model beyond the overlap region where inner and outer scaling laws converge. The model remains applicable in regions exhibiting a hierarchical structure of eddy clusters. The dataset was collected at Nodule Point, located on the eastern side of Marrowstone Island, where the U.S. Navy deployed a small array of Verdant Power
$^{\textrm {TM}}$
turbines. Velocity profiles were recorded within 10–17 February 2011, using an acoustic Doppler current profiler mounted on one of the legs of the tidal turbulence tripod. Also, velocity measurements were obtained during the spring tide (17–21 February 2011) using an acoustic Doppler velocimeter (ADV) positioned at the apex of the tidal turbulence tripod, 4.7 m above the seabed. The site has a depth of 22 m, with a maximum hub-height current of 1.8 m s−1 (Thomson et al. Reference Thomson, Polagye, Durgesh and Richmond2012). The turbulent velocity fluctuation data undergo a series of processing steps to ensure spectral consistency. First, the streamwise,
$u$
, and transverse,
$v$
, velocity fluctuation components are extracted via coordinate transformation. To enforce stationarity, a robust detrending approach based on empirical mode decomposition is applied. This method effectively mitigates low-frequency trends while preserving the integrity of the turbulence spectra by removing the residual component and the largest-scale intrinsic mode function from the velocity fluctuations (Cheng et al. Reference Cheng, Neary and Chamorro2024a
).
We analysed the mean velocity profile to identify the overlap region, where the
$-1$
spectral scaling is typically observed. The ADV measurements were taken at 4.7 m above the seabed, corresponding to approximately
$0.25\delta$
, where
$\delta$
is the boundary-layer thickness. This places the measurements beyond the typical surface layer, which is generally confined within
$\sim 0.1\delta$
. Consequently, Townsend’s AEH no longer holds in this region. This conclusion is further supported by an analysis of the shear stress profile, which confirms that the measurement location lies outside the overlap region.
The experimental streamwise velocity spectra,
$\varPhi _u(f)$
, and their compensated counterparts,
$f\varPhi _u(f)$
, at this location are presented in figure 3(a,b). We observe the persistence of
$-1$
spectral scaling beyond the overlap region, a finding of particular significance. Since the data were collected in the temporal domain, time and frequency arguments replace spatial and wavenumber counterparts in the formulation by using the Taylor’s frozen field hypothesis. The turbulence intensity, defined as the ratio of the root-mean-square velocity fluctuation to the mean flow, is around 9 %, a reasonable level for using this hypothesis. For illustration, the function
$\psi (\boldsymbol{\cdot })$
in the model (2.6) is considered to be a class of covariance functions proposed in Jetti et al. (Reference Jetti, Porcu and Ostoja-Starzewski2023, Reference Jetti, Cheng, Porcu, Chamorro and Ostoja-Starzewski2025b
). The dashed-line representation in figure 3(a,b) indicates
$\widehat {\varPsi }^{\textit{app}}(k)$
with the piecewise approximation in (2.9). The proposed model, (2.6), captures both the low-frequency and high-frequency spectral behaviours and, most importantly, reproduces the
$-1$
spectral scaling region with high fidelity. For comparison, we also present the structure function data in figure 3(c). At small time lags
$\tau$
, the data follow Kolmogorov scaling, and it plateaus at larger lags. It also shows an intermediate region where a logarithmic behaviour emerges as indicated by the grey line.
Figure 3. (a–c) Streamwise velocity statistics for a tidal current 4.7 m above the seabed. (a) Spectrum
$\varPhi _u(f)$
with
$f_1$
,
$f_2$
indicating the
$f^{-1}$
region and shaded area marking saturation. (b) Compensated spectrum
$f\varPhi _u(f)$
. (c) Structure function
$S_2(\tau )$
with lag normalised by plateau value
$\tau _p$
. (d–f) Corresponding ABL wind measurements. Black: experiment; red: proposed model; dashed: approximation from (2.9).
We also analysed ABL wind turbulence measurements at a coastal interface of the Yucatán Peninsula, where the convergence of land and ocean surfaces creates a diverse range of surface roughness conditions. This site is particularly well suited to investigate the emergence of the
$-1$
spectral scaling. The proposed formulation accurately captures the
$-1$
spectral scaling when present and robustly models spectra that do not exhibit this scaling, demonstrating the model’s versatility. Figure 3(d–f) presents the streamwise velocity spectrum, compensated spectrum and second-order structure function from a representative ABL wind dataset recorded on 11 May 2011 (22:07–23:54). The data exhibit a clear Kolmogorov scaling at high frequencies and a plateau at low frequencies, indicative of short-range dependence. An extended
$-1$
spectral region is also present. The proposed model accurately captures all three regimes, demonstrating its capacity to represent multiscale behaviour in atmospheric turbulence. Details comparing additional experimental data with the proposed formulation, along with a computationally efficient method for determining its parameters, are provided in Appendix B. There, we illustrate several cases that exhibit extended
$-1$
spectral scaling regions, including data with both short- and long-range dependences. We also examine cases where such scaling is absent. Our results demonstrate the robustness of the model in accurately identifying the
$k^{-1}$
region when present and correctly recognising its absence when unsupported.
4. Conclusions
We have presented a theoretical approach that derives the
$k^{-1}$
spectral scaling in turbulent velocity spectra using random field theory without relying on Kolmogorov’s inertial-range assumptions, specific eddy correlation forms or Townsend’s AEH. By applying asymptotic theorems, we establish a direct connection between the streamwise velocity covariance and its spectral form, showing that the statistical mixing of eddy clusters with distributed integral scales leads to
$k^{-1}$
scaling as an intermediate asymptotic. The model also recovers the logarithmic form of the second-order structure function, consistent with its spectral analogue. Validation against wind and riverine turbulence data confirms that this scaling can persist beyond the surface layer and in the absence of geometrically attached eddies. These insights suggest a more general mechanism for the emergence of
$k^{-1}$
scaling and motivate further investigation into the role of surface roughness, stratification and unsteadiness in shaping spectral energy distributions in geophysical and engineering flows.
Acknowledgements
We thank V. Neary and B.F. Espinoza for providing field data.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Logarithmic dependence of the structure function
The asymptotic behaviour of the covariance function for a fractal dimension
$D$
and Hurst parameter
$H$
is
\begin{align} 1-{\psi }(x) &\to S_{00} (x/c)^{-(1+q)}, \quad x \to 0, \end{align}
\begin{align} {\psi }(x) &\to S_{10} (x/c)^{-(1+p)}, \quad x \to \infty . \end{align}
Note that
$p=1-2H$
and
$q=2D-5$
, as defined in § 2. The constants
$S_{00}$
and
$S_{10}$
are specific to the covariance function
$\psi (x)$
and are related to
$C_{00}$
and
$C_{10}$
. To demonstrate the logarithmic dependence of the combined structure function, we adopt an approximate form of
$\psi (x)$
for analytical tractability
\begin{equation} 1-{\psi }^{\textit{app}}(x) = {\begin{cases} S_{00} (x/c)^{-(1+q)}, & x \leq x_T, \\ 1-S_{10} (x/c)^{-(1+p)}, & x \gt x_T. \end{cases}} \end{equation}
Here,
$x_T$
is the transition point where the two asymptotic forms of
${\psi }(x)$
intersect, satisfying
\begin{equation} S_{00} (x_T/c)^{-(1+q)} = 1-S_{10} (x_T/c)^{-(1+p)}. \end{equation}
The approximate form of the combined structure function,
$2(\sigma ^2 - \varPsi (x))$
, is derived using the piecewise approximation in (A3). For
$x \leq x_T(c_1)$
, the behaviour of all the structure functions corresponding to different
$c$
values in the range
$c_1 \leq c \leq c_2$
is dictated by the fractal dimension
$D$
. Thus
\begin{eqnarray} 2(\sigma ^2 - {\varPsi }^{\textit{app}}(x)) &=& \frac {2\sigma ^2}{\log (c_2/c_1)} \int _{c_1}^{c_2} S_{00} \left (\frac {x}{c}\right )^{-(1 + q)} \nonumber\\[6pt] d(\log c) &=& -\frac {2\sigma ^2 S_{00} x^{-(1 + q)} \big (c_1^{1+q} - c_2^{1+q}\big )}{(1+q)\log (c_2/c_1)} \sim x^{-(1+q)}. \end{eqnarray}
For isotropic turbulence, where
$D = 5/3$
, it follows that
$\forall x \leq x_T(c_1)$
,
$2(\sigma ^2 - \varPsi ^{\textit{app}}(x)) \sim x^{2/3}$
. Similarly, for
$x \geq x_T(c_2)$
, the structure function behaviour is independent of the fractal dimension and is instead determined solely by the Hurst parameter
$H$
. Thus
\begin{eqnarray} 2(\sigma ^2 - \varPsi ^{\textit{app}}(x)) &=& \frac {2\sigma ^2}{\log (c_2/c_1)} \int _{c_1}^{c_2} (1 - S_{10} (x/c)^{-(1+p)} ) \nonumber\\[6pt]d(\log c) &=& 2\sigma ^2 (1 - A x^{-(1+p)}) \sim 2\sigma ^2, \end{eqnarray}
where
$A = {S_{10} (c_2^{1+p} - c_1^{1+p} )}/{(1 + p) \log (c_2 / c_1)} $
. For
$x_T(c_1) \lt x \lt x_T(c_2)$
, the combined structure function follows
\begin{eqnarray} 2(\sigma ^2 - \varPsi ^{\textit{app}}(x)) &=& \frac {2\sigma ^2}{\log (c_2 / c_1)} \left ( \int _{c_1}^{c_x} ( 1 - S_{10} (x/c)^{-(1 + p)} ) d(\log c) \right. \nonumber\\[6pt] && \left. +\, \int _{c_x}^{c_2} S_{00} (x/c)^{-(q + 1)} d(\log c) \right)\!, \end{eqnarray}
where
$c_x$
is the value of
$c$
for which
$x = x_T$
from (A4). Since an explicit expression for
$c_x$
cannot be obtained from (A4), we approximate
$x_T \simeq 1/k_T$
(see (2.10)), ensuring consistency between the spectral and real-space transition regions (
$k_2 \leq k \leq k_1$
corresponding to
$x_1 \leq x \leq x_2$
). This leads to
\begin{align} 2(\sigma ^2 - \varPsi ^{\textit{app}}(x)) &= \frac {2\sigma ^2}{\log (c_2 / c_1)} \left ( \log \left (\frac {x}{x_1}\right ) + \frac {S_{10}}{1 + p} \left ( \frac {C_{10}}{C_{00}} \right )^{\frac {1 + p}{p - q}} \left ( \left ( \frac {x}{x_1} \right )^{-(1 + p)} - 1 \right ) \right . \nonumber \\ & \left . - \frac {S_{00}}{q + 1} \left ( \frac {C_{10}}{C_{00}} \right )^{\frac {q + 1}{p - q}} \left ( 1 - \left ( \frac {x}{x_2} \right )^{-(q + 1)} \right ) \right)\!. \end{align}
Here, we define
$x_1 := x_T(c_1) \simeq 1/k_1$
and
$x_2 := x_T(c_2) \simeq 1/k_2$
, we note that for
$x_1 \lt x \lt x_2$
, the first and third terms in (A8) are positive, while the second term is negative. However, the overall expression remains positive, with the first term (logarithmic dependence) dominating over the third term (power-law dependence). Thus, the structure function in the transition region can be approximated as
$2(\sigma ^2 - {\varPsi }^{\textit{app}}(x)) \sim B\log {x} + A$
. Equation (A8) provides a close approximation to (A7). For isotropic turbulence, Kolmogorov’s scaling fixes
$D = 5/3$
. As the Hurst parameter approaches 1, the second term in (A8) also begins to exhibit logarithmic dependence on
$x$
, as seen by examining its asymptotic form
\begin{align} \frac {S_{10}}{1 + p} \left (\frac {C_{10}}{C_{00}}\right )^{\frac {1 + p}{p - q}} \left ( \left (\frac {x}{x_1}\right )^{-(1 + p)} - 1 \right ) \ \to \ - S_{10} \left (\frac {C_{10}}{C_{00}}\right )^{\frac {1 + p}{p - q}} \log \left (\frac {x}{x_1}\right)\!. \end{align}
This result follows from the Taylor series expansion
\begin{align} \log {z} = (z-1) - \frac {(z-1)^2}{2} + \frac {(z-1)^3}{3} + \ldots , \end{align}
which converges for
$|z - 1| \leq 1$
. As
$H \to 1$
, we have
$(x/x_1)^{-(1 + p)}=(x/x_1)^{-2(1-H)} \to 1$
, ensuring that the first-term approximation of the logarithm is highly accurate.
Appendix B. Field validation of
$\boldsymbol{k}^{\boldsymbol{-1}}$
scaling
We introduce a model covariance class that captures fractal and Hurst-type behaviours in
$\psi (\boldsymbol{\cdot })$
, enabling robust validation of the proposed formulation (2.6) using experimental data. To extract
$\psi (\boldsymbol{\cdot })$
from covariance data, we consider the first derivative of
$\varPsi (x)$
\begin{equation} \varPsi ^{\prime }(x) = \frac {\sigma ^2}{x\log {(c_2/c_1)}} \left ( \psi (x/c_1) - \psi (x/c_2) \right ), \quad x \geq 0. \end{equation}
Recovering
$\psi (\boldsymbol{\cdot })$
from
$\varPsi ^{\prime }(x)$
is an ill-posed problem due to degeneracy in extracting individual covariance functions from their superposition. To address this, we adopt a covariance model from Jetti et al. (Reference Jetti, Porcu and Ostoja-Starzewski2023, Reference Jetti, Cheng, Porcu, Chamorro and Ostoja-Starzewski2025b
) that accounts for the long-range dependence and fractal properties. Unlike conventional spectral models (e.g. von Kármán, Kaimal), this formulation simultaneously captures both features, yielding a more comprehensive representation of turbulence. Its validity is demonstrated through comparisons with field data from tidal currents and ABL flows, reproducing the expected asymptotic spectral scaling. We consider the isotropic covariance functions introduced in Jetti et al. (Reference Jetti, Porcu and Ostoja-Starzewski2023):
$\psi _{\alpha ,\beta ,\gamma }(x) = ( 1 - ( 1 + (x/c)^{-\gamma } )^{-\alpha } )^{\beta }$
,
$x \geq 0$
.
A constructive proof establishes that
$\psi _{\alpha ,\beta ,\gamma }(\boldsymbol{\cdot })$
is a valid isotropic covariance function in
${\mathbb{R}}^d$
for all
$d = 1,2,\ldots$
if
$\alpha \in (0,1]$
,
$\beta \gt 0$
and
$\gamma \in (0,2]$
. While our focus is on the one-dimensional case (
$d=1$
), we present generalised formulas for completeness. The function
$\psi _{\alpha ,\beta ,\gamma }$
exhibits a short-range-dependence (SRD) when
$\gamma \beta \gt d$
and long-range-dependence (LRD) when
$\gamma \beta \leq d$
. The corresponding fractal dimension is given by
$D = d + 1 - (({\gamma \alpha })/{2})$
, while the Hurst parameter is
$H = 1 - (({\gamma \beta })/{2})$
for LRD cases (
$\gamma \beta \gt d$
).
The spectral density
$\widehat {\psi }_{\alpha ,\beta ,\gamma }(k)$
exhibits the following asymptotic behaviour:
Low wavenumbers (
$k \to 0$
)
\begin{align} &\widehat {\psi }_{\alpha ,\beta ,\gamma }(k) \to C_{00} c (kc)^{\gamma \beta - d}, \quad \text{if } \gamma \beta \in \left (\frac {d-1}{2}, d\right)\!, \end{align}
\begin{align} &\widehat {\psi }_{\alpha ,\beta ,\gamma }(k) \to K_0 c, \quad 0 \leq K \lt \infty , \quad \text{if } \gamma \beta \gt d. \\[9pt] \nonumber \end{align}
High wavenumbers (
$k \to \infty$
):
$ \widehat {\psi }_{\alpha ,\beta ,\gamma }(k) \to C_{10} c (kc)^{-d-\gamma \alpha }$
. The constants
$C_{00}$
,
$K_0$
and
$C_{10}$
depend on
$\alpha$
,
$\beta$
,
$\gamma$
and
$d$
but remain independent of
$c$
. See details in Jetti et al. (Reference Jetti, Porcu and Ostoja-Starzewski2023, Reference Jetti, Cheng, Porcu, Chamorro and Ostoja-Starzewski2025b
). We adopt
$\psi _{\alpha ,\beta ,\gamma }(x)$
to solve (B1), reducing the optimisation problem to determining the optimal values of
$\alpha$
,
$\beta$
,
$\gamma$
,
$c_1$
and
$c_2$
. The best-fit parameters for
$\varPsi ^{\prime }(x)$
are obtained by
\begin{equation} \varPsi ^{\prime }(x) = \frac {\sigma ^2}{x\log {(c_2/c_1)}} ( \psi _{\alpha ,\beta ,\gamma }(x/c_1)-\psi _{\alpha ,\beta ,\gamma }(x/c_2)), \quad x \geq 0. \end{equation}
Next, we first validate this approach using experimental ABL turbulence data from a coastal environment.
We analyse a dataset from a meteorological tower at the coastal interface of the Yucatán Peninsula (
$90^{\circ }02^{\prime }47^{\prime \prime }$
W,
$21^{\circ }09^{\prime }53^{\prime \prime }$
N). We found that sea winds consistently exhibit the
$-1$
scaling for all time scales and seasons, while land winds lack a clear
$-1$
region. We searched the three-year dataset for near-neutral, stationary wind segments with a stable wind direction for at least 60 min. A coordinate transformation aligned the velocity components with the mean wind direction, resulting in streamwise (
$u$
) and cross-wise (
$v$
) components. Stationarity was enforced using empirical mode decomposition (Cheng et al. Reference Cheng, Neary and Chamorro2024a
), which removes low-frequency trends by eliminating the residual and largest-scale intrinsic mode function (Cheng et al. Reference Cheng, Neary and Chamorro2024b
).
For a stationary streamwise velocity profile, the optimal parameters are determined to satisfy (B1), ensuring the positive definiteness of
$\psi _{\alpha ,\beta ,\gamma }(\boldsymbol{\cdot })$
. Here,
$\tau _1$
and
$\tau _2$
replace
$c_1$
and
$c_2$
from the original formulation. Kolmogorov’s scaling is enforced by setting
$\gamma \alpha = 2/3$
while optimising (B4). The spectra are then analysed for long-range dependence, imposing
$\gamma \beta \gt 1$
for short-range-dependent velocity profiles.
We analyse sea and land winds across seasons, focusing first on stationary sea wind datasets that exhibit a well-defined
$-1$
spectral scaling. Figure 4 presents the streamwise velocity spectra,
$\varPhi _u(f)$
, and the compensated spectra,
$f\varPhi _u(f)$
, for three periods. The velocity profile shows short-range dependence, with a spectral plateau at low frequencies. The spectral bounds,
$f_1 := 1/\tau _1$
and
$f_2 := 1/\tau _2$
, effectively define the
$-1$
scaling region, with
$f_u$
marking the integral time scale. The optimisation relies only on the covariance function’s first derivative, without explicitly imposing the
$-1$
scaling, demonstrating its natural emergence.
Figure 4. Streamwise velocity spectra of winds from the sea (SRD) on (a) 12 January 2011 (09:05–11:05), and (b) 11 May 2011 (22:07–23:54).
Next, we analyse long-range-dependent stationary sea winds exhibiting
$-1$
spectral scaling, shown in figure 5 for three periods. The estimated bounds
$f_1$
and
$f_2$
are slightly less precise but still capture the dominant features of the
$k^{-1}$
region. This reduced accuracy probably stems from the difficulty in precisely defining the long-range-dependence constraint and the smoother transition from the
$k^{-1}$
regime to the inertial subrange compared with previous cases.
Figure 5. Streamwise velocity spectra of winds from the sea (LRD) (a) from 01 January 2012 (23:48) to 02 January 2012 (03:50), and (b) from 07 April 2012 (19:11) to 08 April 2012 (00:14).
Finally, we examine land winds with short-range dependence that lack significant
$-1$
spectral scaling, shown in figure 6 for three periods. The optimisation yields nearly identical
$f_1$
and
$f_2$
values, converging to the characteristic frequency
$f_u$
, which corresponds to the inverse of the integral time scale. This demonstrates the robustness of the model in correctly identifying the
$k^{-1}$
region when present, as in sea winds, while also accurately recognising its absence in land winds.
Figure 6. Streamwise velocity spectra of winds from the land on (a) 09 May 2011 (09:22–11:30), and (b) 25 May 2011 (08:30–09:42).
Open access
