2. Theory and methodology

As introduced above, near-wall coherent structures manifest as high- and low-drag events in the turbulent boundary layer (Örlü & Schlatter Reference Örlü and Schlatter2011). Their persistent nature implies that the forces acting on a particle are not random but strongly time correlated. Accordingly, we focus on characterising the time scales, amplitudes and memory properties of these events to link particle motion within the viscous sublayer to the underlying flow topology. We analyse DNS data from the John Hopkins Turbulent Database for a fully developed turbulent channel flow at $\textit{Re}_{\tau } = u_{\tau }h/\nu \approx 1{\,}000$ (where h is the channel half-height) (Graham et al. Reference Graham2016). The simulation domain covers $8\pi h \times 2h \times 3\pi h$ in the streamwise, x, wall-normal, y, and spanwise, z directions, respectively, and fully resolves scales from energy-containing eddies down to the dissipative Kolmogorov scales. The extracted time series utilises the 4000 stored database snapshots, which have a non-dimensional database time step of $\delta t = 0.0065$ . This provides a fine temporal resolution of $\Delta t^+ = 0.325$ in viscous units. Given that the stochastic resuspension framework is formulated as a single-particle model, the particle concentration is assumed to be sufficiently dilute so that the particles do not interact with each other or significantly alter the flow. Under these conditions, the particles respond to the turbulent velocity field without feeding back onto it, consistent with the one-way coupling assumption. We model tungsten particles with $d_{\!p}=13 \rm \,{\unicode{x03BC}} m$ diameter, corresponding to a particle Reynolds number or dimensionless diameter of $\textit{Re}_{\!p}=d_{\!p}^+=d_{\!p}u_\tau /\nu = 0.013$ and a Stokes number of $\textit{St}^+=\rho _{p}(d_{\!p}^+)^2/18\rho _{\!f} = 0.18$ (where $\rho _{p}$ and $\rho _{\!f}$ are the particle’s and fluid’s density, respectively) in wall units. Given these small values of $d_{\!p}^+$ and $\textit{St}^+$ , the one-way coupling assumption remains valid because the system is in the dilute, single-particle limit, where momentum transfer to the bulk fluid is negligible. This is supported by the experimental findings of Quibeuf, Charru & Lacaze (Reference Quibeuf, Charru and Lacaze2020), who demonstrated that heavy, inertial particles (e.g. zirconium) with $\textit{Re}_{\!p}\lt 10$ follow fluid velocity fluctuations within the viscous sublayer without significantly altering the local flow structure. Furthermore, the baseline Markovian stochastic framework (Ben Shlomo et al. Reference Ben Shlomo, Almog, Klausner, Fattal and Berkovich2024) was originally validated against measurements of tungsten particles, confirming the robustness of the one-way coupled stochastic approach for dense particles in this regime. In our stochastic-Lagrangian framework, the simulation is initialised with the particles at rest (their initial velocity is zero) directly on the surface. Because their dimensionless diameter is significantly smaller than the sublayer height, the particles are deeply immersed within the viscous sublayer.

We identify events using an amplitude criterion based on the normalised wall shear stress (WSS), defined as the ratio of the instantaneous WSS value to its time-averaged value, $\tau _{w}/\left \langle \tau _{w} \right \rangle$ . An amplitude criterion is adopted, where events are classified as high or low drag when the normalised WSS exceeds a $\pm$ 5 % threshold. This specific value provides a physical balance: it isolates distinct energetic structures from background noise while avoiding the aggressive filtering of larger amplitude thresholds (e.g. $\pm 10\,\%$ ) or strict minimum-time criteria. This allows us to capture the full spectrum of short-time-scale sublayer events relevant to particle resuspension, while improving the statistical robustness (Pereira, Thompson & Mompean Reference Pereira, Thompson and Mompean2019). Notably, while varying these thresholds alters the absolute numerical value of the resulting event decay rates ( $\lambda$ ), the underlying qualitative exponential nature of the event distributions remains consistent. A recent research demonstrated that adopting inner scaling ( $\mathrm{\Delta }t^{+} = \mathrm{\Delta }tu_{\tau }^{2}/\nu$ ) renders the fraction of time spent in high- and low-drag events independent of the friction Reynolds number $\textit{Re}_{\tau }$ (Agrawal et al. Reference Agrawal, Ng, Davis, Park, Graham, Dennis and Poole2020). Based on this finding, we employed inner time scaling in our analysis. To quantify memory effects in the WSS signal, we estimate the Hurst exponent (H). The Hurst exponent characterises the long-range temporal dependence of a signal, where $H=0.5$ indicates an uncorrelated Markovian process and $0.5 \lt H \lt 1$ indicates persistent flow memory. We estimate $H$ using rescaled range (R/S) analysis (Hurst Reference Hurst1951). Because the DNS data represent a statistically stationary flow, R/S analysis provides a robust estimation while avoiding the artefactual curvature and detrending errors known to affect other methods, such as detrended fluctuation analysis, when applied to stationary or weak-nonlinear signals (Bryce & Sprague Reference Bryce and Sprague2012). The R/S analysis is described as follows. The sample length, n, defines a temporal window over which the near-wall turbulent fluctuations are evaluated. For each block of length n, the adjusted range R is computed as the difference between the maximum and minimum of the cumulative deviation from the block mean, capturing the largest persistent excursion of the shear stress relative to the local average. The standard deviation S represents the root-mean-square magnitude of the fluctuations within the same window, reflecting the local turbulent intensity. The ratio R/S thus provides a normalised measure of persistence relative to fluctuation strength, and its scaling with n yields the Hurst exponent H, which characterises the temporal correlation and memory in near-wall turbulence. Hence, H is obtained through a linear regression of the data according to the linear equation $\log (R/S) = \log (c) + H \log (n)$ , where c is a constant. Importantly, because this analysis evaluates the continuous signal, the extracted Hurst exponent is fundamentally independent of the amplitude threshold used for discrete drag event identification.

Our non-Markovian generalisation builds upon our previously introduced Lagrangian stochastic Markovian model, which demonstrated good agreement with experimental data (Ben Shlomo et al. Reference Ben Shlomo, Almog, Klausner, Fattal and Berkovich2024). The model (depicted in figure 1) assumes a particle that rolls in a single direction (Henry et al. Reference Henry, Minier and Brambilla2023), with its motion governed by a moment balance, M, about a pivot point

(2.1) \begin{equation} M \equiv I\frac {\mathrm{d}\omega }{\mathrm{dt}} = \textit{bF}_{\!D} + \frac {a}{2}F_{\!L} - r_{a}F_{a} - \frac {a}{2}\textit{mg}, \end{equation}

where I is the particle moment of inertia, approximated as a sphere and evaluated using the parallel axis theorem as $I = (7/20) m d_{\!p}^{2}$ . The particle diameter and mass are denoted by $d_{\!p}$ and m, respectively, g is the gravitational acceleration, $\omega$ is the particle angular velocity, a is the distance between consecutive surface asperities, r_a is the lever arm of the adhesion force and $b = \sqrt {(d_{\!p}/2)^{2} - (a/2)^{2}}$ . The hydrodynamic drag force is modelled using the Stokes formulation $F_D = 6\pi \rho _{\!f}\nu (d_{\!p}/2)\langle u \rangle f_x$ , where $\langle u \rangle$ is the fluid’s mean velocity in the streamwise direction at the particle’s position and $f_x=1.7$ is O’Neill’s wall-correction factor (Stokes Reference Stokes1851; O’Neill Reference O’Neill1968). The hydrodynamic lift force is modelled using the Mollinger & Nieuwstadt (Reference Mollinger and Nieuwstadt1996) formulation, $F_L = 5.3\pi \rho _{\!f}\nu ^2(d_{\!p}^+/2)^{1.87}f_y$ , where $d_{\!p}^+$ is the dimensionless particle diameter and $f_y=3.39$ is Maude’s wall-correction factor (Maude Reference Maude1963). The mean adhesion force is calculated using the Rabinovich approximation, which describes the interaction potential between a spherical particle and a rough surface $F_{a}=(A_{\!H}d_{\!p}/12z_0^2)[1/(1+29d_{\!p}R_q/a^2)+1/(1+1.82R_q/z_0)^2]$ , where $A_{\!H}$ is the Hamaker constant, $z_0$ is the distance of the closest approach between surfaces and $R_q$ is the root mean square roughness of the surface (Rabinovich et al. Reference Rabinovich, Adler, Ata, Singh and Moudgil2000). The parameters used to calculate the hydrodynamic forces are evaluated from the fluid’s or the particle’s physical properties, while the adhesion force parameters are evaluated from the surface characteristics which can be measured by atomic force microscopy (Peillon et al. Reference Peillon, Gélain, Payet, Gensdarmes, Grisolia and Pluchery2022).

Figure 1.

Schematic diagram of the forces and physical parameters appearing in (2.1).

Within this framework, applying Reynolds decomposition separates the overall angular velocity $\omega$ into mean ( $\langle \omega \rangle$ ) and fluctuating ( $\omega '$ ) components. The mean component is determined by the forces in (2.1), while the fluctuating component is modelled as a stochastic process. Consistent with K41 theory (Kolmogorov Reference Kolmogorov1941), this model uses a standard Ornstein–Uhlenbeck process (OUP) (Risken Reference Risken1989)

(2.2) \begin{equation} \mathrm{d}\omega '(t) = -\omega '(t)\frac {\mathrm{d}t}{T} + \sqrt {\frac {2\sigma ^2}{T}} \, \mathrm{d}W(t), \end{equation}

where T is the integral time scale, $\sigma$ is a parameter that controls the variance and the autocorrelation of the process and $W (t )$ denotes a Wiener process (Risken Reference Risken1989). The random term ${\rm d}W(t)$ is expressed as $\sqrt {{\rm d}t}\xi$ , with $\xi$ being a non-correlated white noise sampled from a normalised Gaussian distribution independently at each time step. The resuspension criterion is fundamentally tied to the particle’s spatial correlation with the surface geometry. During a characteristic time interval $\Delta t$ , a rolling particle undergoes a displacement of $(\omega d_{\!p} / 2) \Delta t$ . If a particle rolls over the pivot, it will leave its original place and the original asperity can no longer affect the rolling action, so the rolling action will lose its correlation with the surface when the distance of the particle displaced is larger than $a$ , marking full detachment and resuspension. This kinematic condition yields the critical angular velocity threshold $\omega _c = 2a / (d_{\!p} \Delta t)$ . While the numerical simulation results naturally depend on $\Delta t$ , this parameter is not an arbitrary integration step. It strictly corresponds to the turbulent viscous time scale, $\Delta t = T_0 = \nu / u_\tau ^2$ . As demonstrated by Fu et al. (Reference Fu, Chao, So and Leung2013), evaluating the threshold over $T_0$ properly accounts for the physical turbulent energy-transfer time scale, rather than numerical artefact. It is important to highlight that, according to Fu et al. (Reference Fu, Chao, So and Leung2013), the parameter $\sigma$ depends on a model constant, $C_{0}$ , via the relationship $\sigma = \sqrt {2 C_0 T_0 \langle M{\prime }^2 \rangle } / I$ . Here, $\langle M'^2 \rangle$ is the variance of the fluctuating hydrodynamic moment, which can be computed from the local fluid Reynolds stresses ( $\langle u^2 \rangle$ , $\langle v^2 \rangle$ , $\langle uv \rangle$ , where $u$ and $v$ are the fluid’s velocity components at the streamwise and wall-normal directions, respectively) evaluated at the particle’s position (the full analytical expansion is provided in the supplementary material, available at https://doi.org/10.1017/jfm.2026.11435). This leaves $C_0$ as the sole free parameter used to calibrate the magnitude of the stochastic forcing. While its precise physical interpretation is unclear, based on our hypothesis, we assume that the physical meaning of this constant reflects the intermittent dynamics in the viscous sublayer, induced by near-wall coherent structures.

The observed penetration of intermittent coherent structures into the viscous sublayer seems to contradict the assumption of uncorrelated white noise, as it may not adequately capture the inhomogeneous dynamics of near-wall turbulence. Faranda et al. (Reference Faranda, Pons, Dubrulle, Daviaud, Saint-Michel, Herbert and Cortet2014) demonstrated that Hurst exponents provide a useful measure of deviations from the K41-based Markovian description in inhomogeneous and anisotropic turbulent flows, linking these deviations to the presence of intermittent coherent structures. Motivated by this insight, we propose a generalised non-Markovian model in which the standard OUP is replaced with a fractional Ornstein–Uhlenbeck process (fOUP) (Cheridito, Kawaguchi & Maejima Reference Cheridito, Kawaguchi and Maejima2003). This formulation is consistent with the refined K62 theory (Kolmogorov Reference Kolmogorov1962) and is obtained by replacing the Wiener process increment $\mathrm{d}W(t)$ in (2.2) with an increment of fractional Brownian motion (fBm), $\mathrm{d}B^{H}(t)$ . The governing stochastic differential equation for the fluctuating component becomes

(2.3) \begin{equation} \mathrm{d}\omega '(t) = -\omega '(t)\frac {\mathrm{d}t}{T} + \sqrt {\frac {2\sigma ^2}{T}} \, \mathrm{d}B^{H}(t), \end{equation}

Here, the Hurst exponent H (estimated from DNS data in the results section) governs the memory of the process. For 0.5 $\lt$ H $\lt$ 1, the process becomes persistent, meaning that fluctuations tend to continue in the same direction over time, thus capturing the influence of coherent structures (Faranda et al. Reference Faranda, Pons, Dubrulle, Daviaud, Saint-Michel, Herbert and Cortet2014; Hung & Tsai Reference Hung and Tsai2024). The fBm can be represented as a cumulative integral of a Gaussian random process

(2.4) \begin{align} B^{H}(t) = \frac {1}{\varGamma \left ( H + \dfrac {1}{2} \right )} \int _{-\infty }^{0} \left [ (t-s)^{H-\frac {1}{2}} - (-s)^{H-\frac {1}{2}} \right ] \mathrm{d}W(s) + \frac {1}{\varGamma \left ( H + \dfrac {1}{2} \right )} \int _{0}^{t} (t-s)^{H-\frac {1}{2}} \mathrm{d}W(s), \end{align}

where $\varGamma ( H )$ is the gamma function (Mandelbrot & Van Ness Reference Mandelbrot and Van Ness1968). The derivation of the conditional probability density function (PDF), $P(\omega ', t \mid \omega '_0, t_0)$ (with $t \gt t_{0}$ ), follows directly from the Gaussian nature of the underlying noise. Because the governing linear stochastic differential equations are driven by standard and fBm, the resulting fluctuating angular velocity $\omega '$ is strictly a Gaussian process. By taking the exact time derivative of the Gaussian PDF and rearranging the terms (Risken Reference Risken1989), one obtains a generalised Fokker–Planck equation. The Fokker–Planck equation is not restricted to Markov processes, as both the Gaussian Markov process and fBm satisfy the same 1-point Fokker–Planck partial differential equation (McCauley Reference McCauley2012)

(2.5) \begin{equation} \frac {\partial }{\partial t} P(\omega ', t \mid \omega '_0, t_0) = \left [ -\frac {\partial }{\partial \omega '} D_1(\omega ', t) + \frac {\partial ^2}{\partial {\omega '}^2} D_2(\omega ', t) \right ] P(\omega ', t \mid \omega '_0, t_0), \end{equation}

Here, the drift coefficient $D_1$ represents the deterministic relaxation of the system, taking the form $D_1(\omega ', t) = -\omega '/T$ . The effective dispersion coefficient $D_2$ is explicitly dictated by the time derivative of the process variance, $V(t) = \langle {\omega '}^2(t) \rangle$ , defined as $D_2(t) = (1/2) {\rm d}V(t)/{\rm d}t + V(t)/T$ . Given the noise amplitude of $\sqrt {2\sigma ^2/T}$ , the standard OUP evaluates to a constant dispersion coefficient $D_2 = \sigma ^2/T$ . In contrast, the fOUP is distinguished by a uniquely defined, time-dependent $D_2(t)$ (scaling as $t^{2H-1}$ ) that encapsulates the anomalous dispersion and temporal memory of the sublayer. For additional details on the effect of the Hurst parameter on the particle dynamics and the numerical algorithm of the non-Markovian model, see the supplementary material.