<figure><img src="image_1.png" /></figure>

HAL
open science

# Hybrid stochastic-structural modelling of particle-laden turbulent flows based on wavelet reconstruction

Roxane Letournel, Clément Morhain, Marc Massot, Aymeric Vié

► To cite this version:

Roxane Letournel, Clément Morhain, Marc Massot, Aymeric Vié. Hybrid stochastic-structural mod-
elling of particle-laden turbulent flows based on wavelet reconstruction. 2025. hal-05300188v2

HAL Id: hal-05300188

https://hal.science/hal-05300188v2

Preprint submitted on 15 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Hybrid stochastic-structural modelling of particle-laden turbulent flows based on wavelet reconstruction

Roxane Letournel<sup>a,b,c,*</sup>, Clément Morhain<sup>a,b,c,**</sup>, Marc Massot<sup>c</sup>, Aymeric Vié<sup>a,b</sup>

<sup>a</sup>Laboratoire EM2C UPR 288, CNRS, CentraleSupélec, Université Paris-Saclay, 3, rue Joliot-Curie 91192
Gif-sur-Yvette cedex France

<sup>b</sup>Fédération de Mathématiques de CentraleSupélec, CNRS, 3, rue Joliot-Curie 91192 Gif-sur-Yvette cedex
France

<sup>c</sup>Centre de Mathématiques Appliquées, CNRS, École polytechnique, Institut Polytechnique de Paris, Route de Saclay, 91120 Palaiseau Cedex, France

## Abstract

Reduced-order modelling and simulation of turbulent particle-laden flows is required in numerous configurations, where the whole spectrum of turbulent scales through DNS is out of reach. Whereas structural or stochastic models have been derived in order to provide a synthetic turbulent model for the non-resolved scales of the gaseous flow field, reproducing preferential concentration is challenging because it requires capturing both spatial and temporal correlations. We present a novel reduced-order framework that overcomes this limitation by combining wavelet-based structural modelling with stochastic evolution. Using compactly supported divergence-free wavelets within a multiresolution analysis, the method provides direct control over spatial structures and correlations of synthetic multiscale velocity fields. In particular, a dedicated procedure enables to enforce a prescribed turbulent energy spectrum despite the non-local contribution in Fourier space of the wavelet basis functions. The stochastic evolution of wavelet coefficients further ensures consistent temporal correlations. The proposed framework is evaluated in homogeneous isotropic turbulence under a fully reduced setting, where all turbulent scales must be provided by the model. Results show that it accurately reproduces preferential concentration across a wide range of Stokes numbers, achieving closer agreement with DNS data than classical Fourier-

*Now at SAFRAN Tech, Châteaufort, France

**Corresponding author

Email address: clement.morhain@centralesupelec.fr (Clément Morhain)
based kinematic simulations. This establishes a versatile and physically consistent
turbulence model that combines structural fidelity with stochastic dynamics, offering a
new tool to investigate particle–turbulence interactions.

***Keywords:*** Synthetic turbulence, kinematic simulation, divergence-free wavelets, stochastic modelling, structural modelling, particle-laden flows

**1. Introduction**

Due to the multiscale nature of turbulent particle-laden flows, their modelling
and simulation require model-reduction strategies to limit the number of degrees of
freedom in the system, since DNS of such flows is often out of reach, especially in
applications. Nowadays, Large Eddy Simulation (LES) and Reynolds-Average Navier
Stokes Simulation (RANS) are the most widely used Reduced-Order modelling and
Simulations (ROS) to understand and predict turbulent flows. However, averaging
or filtering operations hinder particle interactions with unresolved fluid structures,
whereas the dynamics of the disperse phase depends on its coupling with small-scale
turbulence at the location of individual particles. This loss of small-scale information
therefore necessitates additional modelling efforts. Challenges for such models include
obtaining accurate statistics of the disperse phase and reproducing particle preferential
concentration, hence correctly accounting for the coupling terms between the two phases
[1, 2]. These coupling terms are crucial in turbulent two-phase flows as it governs the
interaction between the two phases, which is often of interest for the studied physics
(such as liquid spray atomization and mixing, droplet combustion, heat transfer by solid
particles, bubble-induced turbulence, etc.). As this coupling can act at all flow scales
[3], unresolved-scale models must accurately predict the particle spatial distribution at
both resolved and unresolved scales. The most classical modelling approach combines
a single-phase ROS model to describe the gas phase with an additional sub-filter scale
(SFS) model for the particulate phase [4]. These models are suitable for both LES and
RANS approaches [5], and we thus refer to them as ROS. Different modelling strategies,
such as stochastic and structural models, can be employed to account for the effects of
unresolved small-scale interactions.
Stochastic models aim at mimicking the effect of sub-filter scale flow on particle paths by introducing an additional random process acting directly on the particle velocity, acceleration or dissipation along their trajectories. From this arises the general formalism for PDF, which has been developed by many authors [6–11]. In most cases, stochastic models are initially developed for fluid particles to reproduce specific statistics of velocity fluctuations. Subsequently, these models are extended to inertial particles by modelling the fluid velocity seen by particles on their trajectory [12–20]. Such models can retrieve the right amount of particle kinetic energy and any one-point statistics by construction. Indeed, stochastic model parameters can tune the modelled processes’ expectations, variances, and covariances. Letournel [2] has verified the adequacy of existing models to Kolmogorov theories of turbulence and proposed a classification based on the order of the model (velocity, acceleration, etc.) and the type of associated stochastic model (Gaussian, conditioned Gaussian, non-Gaussian).

However, stochastic models rely on uncorrelated Gaussian noises and therefore assume a Markovian nature for the velocity seen by the particles or their acceleration [6, 21, 22]. Since they cannot enforce spatial correlations, very close particles can see different noises. Pozorski and Apte [17] observed that for large-inertia particles, the diffusion introduced by the stochastic models could compensate for the lack of randomization from small scales, roughly reproducing the preferential concentration observed in Direct Numerical Simulation (DNS). Furthermore, these models cannot accurately recreate small-scale patterns using coarse sub-filter information alone for particles with low inertia. Squires and Eaton [23] suggested that two-point statistics can only be found in structural models to ensure particle-fluid and particle-particle velocity correlations. This motivates the need for a structural approach to reproduce turbulence spatial and temporal organisation.

Structural models consist of reconstructing a divergence-free velocity field that mimics the spatial and temporal correlations of the fully resolved turbulent flow, allowing particles to interact with coherent structures across scales. Unlike stochastic models, this implies reconstructing information in the whole space, not only for each particle. This method can preserve the spatial correlation, so neighbouring particles
can experience similar fluid velocities. Structural modelling has been envisaged in the
literature using the following options:

* The Approximate Deconvolution Method (ADM) uses a deconvolution operator to capture fine-scale turbulent structures typically unresolved in coarse-grid simulations. Studies using ADM as a subgrid closure for particle motion equations have reported a significant improvement in specific statistics of fluid velocity components, such as the root-mean-square [24] and partial recovery of subgrid kinetic energy [25]. However, ADM's efficiency diminishes in flows with higher Reynolds numbers, necessitating larger filter widths and cut-off length scales, which increases computational costs. While ADM has shown applicability to incompressible homogeneous isotropic turbulence, its effectiveness in other flow configurations or phenomena like multiphase interactions may vary, requiring further considerations or modifications. Finally, ADM's capability to reconstruct scales is limited, as it can only effectively capture those close to the cut-off scale, resulting in potential challenges in accurately predicting flow characteristics such as preferential concentration [26].

* Fractal interpolation consists in reconstructing the velocity field on a finer mesh by interpolating the filtered field. The interpolation, which depends on the fractal dimension of the signal, should either be known or assumed a priori. More details on this technique can be found in the literature [27, 28]. Marchioli et al. [24] considered models based on fractal interpolation and concluded that this technique cannot systematically recover the amount of sub-filter scale turbulent kinetic energy for the fluid velocity field and hence has a hard time recovering particle segregation.

* Kinematic simulation (KS) is a computationally attractive method to reconstruct the small-scale field. This model has been introduced by Kraichnan [29] and extended by Fung et al. [30]. By imposing a turbulent kinetic energy spectrum, a random velocity field is generated using orthogonal and random Fourier modes. Malik and Vassilicos [31] conducted an exhaustive comparison between DNS
and KS and validated the model regarding fluid particle dispersion in homogeneous isotropic turbulence. Other studies focused on particle one-, two- and multi-particle dispersion with successful results of the KS [32–34]. KS has been extended to LES application and adapted to be a sub-grid model for particles by Flohr and Vassilicos [35], also combined with ADM [36]. They have coupled the LES and KS at the energy level by matching the dissipation rate. In this context, Ray and Collins [37] have shown excellent prediction for clustering of inertial particles with $St \ge 2$. Furthermore, Hausmann et al. [38] have extended this approach to two-way coupling, showing improved results in preferential concentration. However, the lack of sweeping of smaller-scale eddies by the larger-scale flow hinders the recovery of particle segregation for low-inertia particles. Additionally, in the study conducted by Murray et al. [39], the temporal evolution of the synthetic fluid has not been introduced to avoid unwanted oscillations in Lagrangian statistics. They investigated Eulerian and Lagrangian statistics of the flow field given by KS and also observed that such an approach fails to retrieve low-inertia particle segregation.

To summarize, while stochastic models can reproduce specific one-point statistics of velocity fluctuations, they have limitations in capturing spatial correlations and accurately recreating small-scale patterns. On the other hand, structural models offer improved capabilities in reproducing two-point statistics of fluid velocity components and preserving spatial correlations. However, their ability to accurately predict preferential concentration can be limited, especially for low-inertia particles. Therefore, while these models provide valuable insights into turbulent flows and particle behaviour, additional considerations and modifications are necessary to fully reproduce particle statistics in a wide range of flow scenarios.

In the present contribution, we introduce a novel reduced-order framework that overcomes the previously identified limitations by combining wavelet-based structural modelling with stochastic evolution. This new framework relies on compactly supported divergence-free wavelets [40] within a multiresolution analysis. Wavelets, with
compact support, are a powerful tool [41, 42] that, similar to sinusoidal functions, allow decomposing any function in a Hilbert space [43, 44]. The basis of such a theoretical approach and related assumptions are recalled. Relying on a few key mathematical properties, which are investigated and presented, the method provides direct control over spatial structures and correlations of synthetic multiscale velocity fields. Despite the non-local spectral contribution of the wavelet basis, it is possible to accurately reproduce a prescribed energy spectrum by modelling the wavelet coefficients as in-dependent random variables. The stochastic evolution of wavelet coefficients further ensures consistent temporal correlations. Our new approach is evaluated on fully re-duced turbulence cases for which all interactions between particles and turbulence have to be modelled. The impact of the new Wavelet-based Stochastic-Structural Model (Wb-SSM) is first justified in 2D and 3D using turbulence statistics. Then, the disperse phase statistics are investigated, and we focus on particle preferential concentration (PPC). Results are compared with DNS and Fourier-based Kinematic Simulation (Fb-KS) and thus assess the potential of the proposed framework. This approach paves the way for compressing the velocity field with controlled error and reducing the number of degrees of freedom, as well as for introducing intermittency [45], which is left for future investigation and discussed in the conclusion.

Let us underline such an approach relying on wavelets has already been used in the literature. Benzi et al. [46] are probably the first to use wavelets with random coefficients to generate a 1D multiaffine signal in the context of turbulence. In the same spirit, other attempts only focusing on either temporal or spatial 1D reconstruction followed [47–49]. The key initial building block of our work is to be found in Letournel [2], who extended the use of wavelets to the 3D spatio-temporal reconstruction of turbulence. The proposed strategy involves incrementally incorporating key features of turbulence to construct the simplest possible model capable of reproducing low-order statistical properties with fidelity and has led to the present contribution. Demonstrating the versatility of this subgrid-scale modelling framework, a similar rationale was later adopted by Hausmann et al. [50], who developed subgrid-scale models for large eddy simulations (LES) using the same wavelet-based foundation. However, in contrast to the present approach, they opted to a linearised version of the Navier-Stokes equations
themselves, employing a Galerkin-type projection to model the evolution of the wavelet
coefficients. They concluded that their results for particle statistics with wavelet-
enriched LES are somehow improved in comparison with the Fourier-enriched LES.
Compared to their Galerkin-type projection of linearised Navier-Stokes, which requires
several simplifications to be used in practice, we propose a rigorous analysis of the
use of wavelet to reproduce a given energy spectrum as well as a simpler temporal
evolution of the velocity field using a Langevin stochastic equation for the amplitude of
each wavelet, allowing a direct control on time correlations.

The paper is organised as follows: the baseline modelling framework of turbulent
particle-laden flows is presented in Section 2, and we point out the advantages and
drawbacks of Fb-KS to establish a list of specifications for the new framework. Within
Section 3, we present the multiresolution analysis framework, we will rely on in con-
structing the Wb-SSM. Section 4 introduces the Wb-SSM, detailing the construction of
the velocity field, the procedure to enforce a prescribed energy spectrum, the statistical
characterisation of the resulting fields and their connection to turbulence scales, and
the stochastic evolution of wavelet coefficients. The proposed framework is validated
in Section 5, where the Wb-SSM is first assessed as a surrogate for single-phase turbu-
lence. Point-particle behaviour in such flows is then examined, and particle preferential
concentration is compared against DNS and Fb-KS, thus proving the potential of the
approach.

## 2. Turbulent particle-laden flows and the challenges of ROS

2.1. *Baseline modelling framework*

We consider monodisperse solid spherical particles with diameter $d_p$ and density $\rho_p$. The particle diameter is smaller than the Kolmogorov length scale; thus, particle-resolved DNS is unnecessary [51], and a point-particle approximation is adopted. Its position is defined in space by coordinates $\mathbf{x}_p$ of the particle centre of mass, and it has a velocity $\mathbf{v}_p$. The density ratio between the particles and the fluid is large, and the Reynolds number of the particle is smaller than one. The volume fraction is small enough to neglect particle-particle collisions. Therefore, many of the forces in the
original Basset-Boussinesq-Oseen equation are assumed negligible compared to the drag force. Furthermore, if the drag force obeys Stokes' law, the equation of motion is linear in the velocity difference between fluid and particle. Then the motion equations of particles are written as:

$$ \left\{ \begin{aligned} \frac{d\mathbf{x}_p}{dt} &= \mathbf{v}_p, \\ \frac{d\mathbf{v}_p}{dt} &= \frac{\mathbf{u}_{f@p}(t) - \mathbf{v}_p(t)}{\tau_p}, \end{aligned} \right. \qquad (1) $$

where $\mathbf{u}_{f@p}(t)$ is the undisturbed fluid velocity evaluated at particle position. The particle relaxation time $\tau_p$ is defined as $\tau_p = \frac{\rho_p d_p^2}{18\mu}$ where $\mu$ is the dynamic viscosity of the fluid. The Stokes number is defined as $St = \tau_p / \tau_\eta$ where $\tau_\eta$ is the Kolmogorov timescale.

In the point-particle limit and in the context of dilute flows, the undisturbed fluid velocity can be considered equal to the fluid velocity at the particle location, i.e., $\mathbf{u}_{f@p}(t) = \mathbf{u}(\mathbf{x}_p(t), t)$. The fluid velocity follows the Navier-Stokes equation, which we here consider in the incompressible limit ($\partial_{x_i} u_i = 0$) for a Newtonian fluid:

$$ \frac{\partial u_i}{\partial t} + \frac{\partial}{\partial x_j} (u_i u_j) = -\frac{1}{\rho} \frac{\partial P}{\partial x_i} + \frac{\mu}{\rho} \left( \frac{\partial^2 u_i}{\partial x_j \partial x_j} \right), \qquad (2) $$

where $\rho$ is the carrier phase density and $P$ its pressure.

## 2.2. Reduced-Order modelling and Simulation (ROS)

The multiscale nature of single-phase turbulent flows does not allow to simulate all the scales of the problem, and Direct Numerical Simulation (DNS) is therefore too expensive. Reduced-order modelling strategies have therefore been developed to reduce computational costs while reliably retaining the main features of the flow. This reduction is performed using ensemble-averaging or filtering operator leading respectively to the well-known Reynolds Averaged Navier-Stokes (RANS) or Large Eddy Simulation (LES) approaches [52].

The common ground of all reduced-order models is to decompose the flow field into a resolved part $U$ and an unresolved part $u'$: In the present contribution, we will focus on the most complex scenario in which the flow field is fully reduced, and all scales have
to be modelled. As the focus here is on being able to predict particles evolution, we first look at the impact of the reduction on Eq. (1):

$$ \frac{d\mathbf{v}_p}{dt} = -\frac{\mathbf{v}_p(t)}{\tau_p} + \frac{\mathbf{u}(\mathbf{x}_p, t)}{\tau_p}. \quad (3) $$

It is the last term of the RHS that requires modelling. For stochastic models, this term is replaced by a random contribution. By specifying the adapted model, retrieving the one-point statistics of the disperse phase is possible. However, they cannot retrieve two-point statistics, which are inherent to the spatial correlation of the fully resolved field. On the other hand, structural models focus on a model for the Eulerian velocity field $\mathbf{u}(\mathbf{x}, t)$. As presented in the introduction with the related state of the art, one of the most promising approaches to retrieve two-point statistics relies on Kinematic Simulations (KS), which fully reconstruct $\mathbf{u}(\mathbf{x}, t)$. By doing so, spatial correlations can be retrieved, a key issue to correctly reproduce preferential concentration of particles. In the literature, KS essentially represents turbulence as a divergence-free sum of Fourier modes with prescribed energy spectrum and oscillatory time evolution. In the following section, we outline the principles of this model and discuss its main limitations.

## 2.3. Fourier-based kinematic simulation

The Fourier-based KS (Fb-KS) Kraichnan [29], Fung et al. [30] has been widely studied in the literature as a model for stationary Homogeneous Isotropic Turbulence (HIT) [39, 53–56]: it consists of a superposition of sine and cosine functions, such as Fourier series synthesis:

$$ \mathbf{u}(\mathbf{x}, t) = \sum_{n=1}^{N_{\kappa}} \sum_{i=1}^{n_n} a_{n,i} \cos(\kappa_{n,i} \cdot \mathbf{x} + \omega_n t) + b_{n,i} \sin(\kappa_{n,i} \cdot \mathbf{x} + \omega_n t), \quad (4) $$

where $N_\kappa$ is the number of distinct wave numbers, $\kappa_n$ and $\omega_n$ are the n-th wave number vector and frequency, and $|a_{n,i}|$ denotes the amplitude of the i-th mode in the n-th shell. The wave number vectors $\kappa_{n,i} = \kappa_n e_i$ are chosen according to a geometric sequence for $\kappa_n$, ranging in the inertial range and $e_i$ is a unit vector in a random direction (i.e, a vector uniformly distributed over the unit sphere). Coefficients $a_{n,i}$ and $b_{n,i}$ are chosen to have random directions orthogonal to $e_i$ to ensure divergence-free of the velocity field.
In order to reproduce an energy spectrum $E(\kappa)$ prescribed by a model, the distribution of modal amplitudes in the Fourier-based KS can be derived analytically using Parseval's identity. In this framework, the total turbulent kinetic energy contained within a spectral shell centred at $\kappa_n$ is directly related to the squared amplitudes of the corresponding Fourier modes. Specifically, exploiting the statistical equivalence between $|a_{n,i}|$ and $|b_{n,i}|$, Parseval's theorem yields:

$$ \sum_{n=1}^{N_{\kappa}} E(\kappa_n) \Delta\kappa_n = \sum_{n=1}^{N_{\kappa}} \sum_{i=1}^{n_n} \mathbb{E}(|a_{n,i}|^2) = \sum_{n=1}^{N_{\kappa}} \sum_{i=1}^{n_n} \mathbb{E}(|b_{n,i}|^2), \quad (5) $$

where $\mathbb{E}(\cdot)$ denotes the expectancy. Now, considering only the $n_n$ Fourier modes contributing to the energy at a given wave number shell centred on $\kappa_n$, we write:

$$ E(\kappa_n) \Delta\kappa_n = \sum_{i=1}^{n_n} \mathbb{E}(|a_{n,i}|^2). \quad (6) $$

Assuming an equipartition of energy among these $n_n$ modes, we arrive at the following formula for the variance of the normally distributed mode amplitudes:

$$ \mathbb{E}(|a_{n,i}|^2) = \frac{E(\kappa_n) \Delta\kappa_n}{n_n}, \quad (7) $$

ensuring that the synthesised field reproduces the target energy spectrum $E(\kappa)$ in expectation.

Temporal coefficients $\omega_n$ represent the eddy turnover frequency associated with each wave number. The most common expression is the following [30, 57]:

$$ \omega_n = \lambda \sqrt{\kappa_n^3 E(\kappa_n)}. \quad (8) $$

More recently, Osborne et al. [58] suggested another form to improve the scalings of Eulerian one-point two-time velocity structure functions:

$$ \omega_n = \lambda \sigma_u \kappa_n. \quad (9) $$

In both models, the parameter $\lambda$ controls the overall eddy turnover time. For very small values of $\lambda$, all the eddy turnover times are very large and turbulence is almost frozen.

The use of such a KS, especially in the context where it is used to mimic the turbulent carrier phase seen by inertial particles, reaches some limitations:
<figure><img src="image_2.png" /><figcaption>Figure 1: Autocorrelation of the velocity of fluid particles along their trajectory on a KS of fluid for different values of eddy turnover times.</figcaption></figure>

* Figure 1 shows the autocorrelation of fluid particle velocity for different values of $\lambda$. The oscillating behaviour observed for $\lambda \neq 0$ is highly unrealistic since the autocorrelation function of the particle velocity is typically well described by a decreasing exponential function [52]. In particular, fluid particle velocities should be completely uncorrelated after a long time, and the autocorrelation should tend towards 0. The oscillating and divergent behaviour of the Lagrangian statistics could be attributed to the sinusoidal time evolution of the synthetic field. To mitigate this issue, Murray et al. [39] recommended using $\lambda = 0$

* However, a frozen fluid is not an appropriate turbulence model: without sweeping, the nature of particle clustering is fundamentally changed. In the absence of collisions, low inertia particles will therefore accumulate indefinitely around the equilibrium points of the flow. The lack of temporality could also explain the difficulty in mapping Eulerian and Lagrangian scales [39].

* Enforcing a periodic condition on such a synthetic field considerably reduces the realizable phase space for the $\kappa_n$. The simulation method here considers only $\kappa_n$ wavevectors with integer coordinates. The number of degrees of freedom for the random sampling of the small $\kappa_n$ is very limited in this case, and the obtained energy spectra are discrete and very noisy.

* The Eulerian spectrum of the synthetic fluid is sparse by construction, which is
not the case for real turbulence

The previous list of limitations of the existing approaches naturally suggests specifications of a new framework for KS:

(i) The synthetic fluid is spatially correlated with compactly supported modes;

(ii) It is a divergence-free flow, to recover the correct fluid particle behaviour;

(iii) It can reproduce any energy spectrum;

(iv) The temporal correlation can be fitted to map Eulerian scales with Lagrangian ones;

(v) Particles one and two-point DNS statistics can be recovered;

These specifications will be reached through the use of Multi-Resolution Analysis, the quintessence of which is presented in the next section before introducing the Wb-SSM.

### 3. Multiresolution analysis and wavelet basis

The introduction of multiresolution analysis (MRA) and wavelets constitute a versatile framework to build a ROS satisfying all the previous properties and addressing the limitations identified in Sec. 2.3. This section focuses on the theoretical fundamentals required to build the divergence-free wavelet basis that will be used in the Wb-SSM.

#### 3.1. *Context and theoretical fundamentals of divergence-free wavelet bases*

In the literature, wavelets and multiresolution analysis have been originally introduced at the convergence of fundamental issues of harmonic analysis and in order to tackle the analysis of multiscale and potentially non-smooth fields or time signals encountered in several applications [59–62], one of them being the analysis of turbulent flow fields [44, 62, 63] and their multifractal character. It is naturally connected to the theory of frames in signal processing and Riesz bases [62].

The purpose of wavelets and MRA is to provide an original basis decomposition for function spaces with localized support in space and frequency variables, based on
a hierarchy of subspaces invariant by translation and related to one another through dilatation, thus allowing a multiscale decomposition. Some wavelets can be built in order to have a compact support, and the decomposition of a field on such a basis gives access to local regularity indices (and thus to compression with error control - a key issue in image compression), whereas Fourier coefficients only provide global ones.

**Definition 3.1.** A Riesz basis {$b_n, n \in \mathbb{N}$} of a real Hilbert space E is a set of elements of E verifying:

$$ \forall f \in E, \quad \exists (c_n)_{n \in \mathbb{N}} \in E^{\mathbb{N}} \quad \text{such that} \quad \lim_{n \to +\infty} \left\| f - \sum_{k=0}^{n} c_k b_k \right\| = 0 \quad (10) $$

and the frame inequality:

$$ \exists A, B \in \mathbb{R}_+^* \quad \text{such that} \quad A \sum_{n=0}^{+\infty} |c_n|^2 \le \|f\|^2 \le B \sum_{n=0}^{+\infty} |c_n|^2. \quad (11) $$

Riesz basis functions provide a complete and orthogonal set of basis functions that span the function space. For the purpose of readability and exposition, we consider the space of one-dimensional square-integrable functions $L^2(\mathbb{R})$. MRA is defined as a sequence of embedded sub-spaces that goes from the finest (small scale or high resolution) to the coarsest (large scale or low resolution) and such that if a function belongs to one of these spaces $V_j$, the function scaled by a factor of 2 belongs to the next space. More precisely, we have:

**Definition 3.2.** A sequence of closed sub-spaces $(V_j)_{j \in \mathbb{Z}}$ of $L^2(\mathbb{R})$ is a MRA if $\forall j \in \mathbb{Z}$,

$$ \bullet V_j \subset V_{j+1} $$

$$ \bullet \bigcap_{j \in \mathbb{Z}} V_j = \{0\} \quad \text{and} \quad \overline{\bigcup_{j \in \mathbb{Z}} V_j} = L^2(\mathbb{R}) $$

$$ \bullet \text{ Dilatation: } f \in V_j \Leftrightarrow f(2\cdot) \in V_{j+1} $$

$$ \bullet \text{ Translation: } f \in V_j \Rightarrow \forall k \in \mathbb{Z}, f(-k) \in V_j $$

The invariance by dilatation and translation implies that there exists a function $\varphi \in V_0$, called scaling function, such that the family {$\varphi(\cdot -k)}_{k \in \mathbb{Z}}$ is a Riesz basis. Therefore,
each space $V_j$ is generated by the family of functions $\{\varphi_{j,k}; k \in \mathbb{Z}\}$ with $\varphi_{j,k}(x) = 2^{j/2}\varphi(2^j x - k)$, where the parameter $j$ defines the level of detail.

It follows that any function $f_j \in V_j$ can be decomposed into the basis of functions:

$$f_j(x) = \sum_{k \in \mathbb{Z}} c_{j,k} \varphi_{j,k} = 2^{J/2} \sum_{k \in \mathbb{Z}} c_{j,k} \varphi(2^j x - k). \quad (12)$$

Let us define the basis of the supplementary spaces $W_j$:

$$V_{j+1} = V_j \oplus W_j, \quad (13)$$

where the sum $V_j \oplus W_j$ is direct. Let us introduce a function $\psi$, called mother wavelet such that the family of functions $\{\psi(\cdot - k); k \in \mathbb{Z}\}$ generates a space $W_0$. Here again, the successive spaces $W_j$ can be generated by the families of functions $\psi_{j,k} = 2^{j/2}\psi(2^j \cdot -k)$, the so-called daughter wavelets. Iterating this decomposition of $V_j$ provides the MRA $(V_j)$ and the associated wavelet spaces $W_\ell$:

$$V_j = V_0 \bigoplus_{\ell=0}^{j-1} W_\ell, \quad (14)$$

which allows to write when $j$ tends to infinity, the wavelet decomposition of the whole space $L^2(\mathbb{R})$ for any $J \in \mathbb{Z}$:

$$L^2(\mathbb{R}) = V_0 \bigoplus_{\ell=0}^{+\infty} W_\ell. \quad (15)$$

Thus, it follows that any function $f \in L^2(\mathbb{R})$ can be decomposed into the basis of functions $\{\varphi_k, \psi_{j,k}; j \ge 0, k \in \mathbb{Z}\}$, with $\varphi_k = \varphi(\cdot - k)$ and $\psi_{j,k} = 2^{j/2}\psi(2^j \cdot -k)$:

$$f = \sum_{k \in \mathbb{Z}} c_k \varphi_k + \sum_{j \ge 0} \sum_{k \in \mathbb{Z}} d_{j,k} \psi_{j,k}, \quad (16)$$

where $(d_{j,k})_{k \in \mathbb{Z}}$ are traditionally called the details at level $j$.

Such a wavelet basis and related MRA can be extended in multidimension. The details are not necessary at this level and are provided to the interested reader in Appendix A.

Whereas wavelets and MRA have also been used as an efficient tool for the discretisation and simulation of PDEs<fnref n="1" />, our purpose is neither discretization nor analysis

<footnote n="1">¹Wavelets have been employed as function space bases in various numerical applications. For instance,</footnote>
of velocity fields but synthesizing such fields in ROS. However, general wavelets do not provide divergence-free fields for the velocity description in turbulent flows. Recovering the incompressibility condition on the flow field could be done through the resolution of a Poisson equation [68].

## 3.2. Divergence-free wavelet basis

A more natural strategy relies on divergence-free wavelets [40] and has been successfully used to analyze 2D/3D turbulent flows [41, 69]. With compact support, they can naturally take into account boundary conditions, contrary to classical approaches only applicable for periodic flows [70].

These divergence-free wavelets have been originally introduced by Lemarié-Rieusset [40] in order to provide bases suitable to represent the incompressible Navier-Stokes solution, in two and three dimensions. In that case, Lemarié-Rieusset [71] has shown that compactly supported wavelets, used for non-periodic problems, can only be supported by biorthogonal wavelet bases. This means that the sum $V_j \oplus W_j$ is direct, but not orthogonal, and the choice of $W_j$ is not unique. For the definition and properties of biorthogonal MRA and associated wavelets, the reader is referred to Kahane [72].

Thus, there exists a large body of literature and theory giving access to various compactly supported divergence-free wavelet bases in 2D and 3D and related MRA, in order to analyse / synthesise velocity fields with potentially a sparse / compressed representation. The general formulation of divergence-free wavelets we are working with is given in Appendix A.3. Any divergence-free vector field $\mathbf{u} \in (L^2(\mathbb{R}^n))^n$ can be then decomposed into the divergence-free wavelet basis:

$$ \mathbf{u} = \sum_j \sum_k \sum_\epsilon d_{\text{div}j, \mathbf{k}}^\epsilon \Psi_{\text{div}j, \mathbf{k}}^\epsilon \quad (17) $$

<footnote>they have been used to solve the two-dimensional Navier-Stokes equations via a Galerkin method [64]. In the broader context of numerical analysis, wavelets also facilitate the compression of solutions to partial differential equations (PDEs) with controlled error [65–67]. This approach significantly reduces memory usage and computational time, particularly for problems characterized by front propagation or the presence of singularities.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: (a-d) Scaling functions and divergence-free 2D wavelets for the MRA spline functions</figcaption></figure>

The wavelets are indexed by $\varepsilon \in \{0, 1\}^n$ and are located at position $\ell_j \mathbf{k}$, with $\mathbf{k} \in [[1; 2^j]]^n$. In the following, we consider the biorthogonal splines functions, introduced by Kahane [72] and also extensively studied by Deriaz and Perrier [41]. An example of these divergence-free wavelets $\Psi_{\text{div},j,\mathbf{k}}^{\varepsilon}$ is given in Fig. 2. Those two-dimensional wavelets are obtained by combining spline functions, following Appendix A.1. They have a compact support of characteristic size $\ell_j = \mathcal{L}/2^j$.

In summary, multiresolution analysis and divergence-free wavelet bases provide a rigorous and versatile framework for representing incompressible velocity fields across multiple scales. Their compact support ensures localization in both space and scale, while the divergence-free construction guarantees consistency with the incompressibility constraint of turbulence. These properties make them particularly well suited for the synthesis of turbulent-like flows. In the next section, we exploit these tools to construct a new wavelet-based stochastic-structural model (Wb-SSM) for turbulent-like flows.

## 4. The Wavelet-based Stochastic-Structural Model (Wb-SSM)

The contribution of this work is to propose a new paradigm for synthetic turbulence simulation by replacing the Fourier decomposition of KS with an MRA, and the deterministic oscillatory time evolution with a stochastic process defined on the wavelet coefficients. The use of compactly supported wavelets, unlike global sinusoidal functions, preserves spatial locality, avoids spurious periodic oscillations in temporal statistics, and removes the strong constraint of periodic boundary conditions. This makes the approach naturally adaptable to heterogeneous flows and complex geometries. Moreover, it opens the way to compressed representations of turbulent fields with
sparse data, while maintaining essential statistical properties and enabling an efficient and physically consistent temporal evolution.

In this section, we first introduce the general formulation of the model based on the divergence-free wavelet basis. Then, we explain how to choose the coefficients of the wavelet decomposition in order to reproduce a target energy spectrum, hence allowing the generation of a frozen turbulent-like field. Subsequently, we introduce a temporal evolution of the field, allowing the synthetic velocity to reproduce both spatial and temporal turbulent-like statistics. This strategy emulates key features of the Navier-Stokes dynamics while avoiding the explicit resolution of the governing PDEs. Finally, we examine how the scales represented in the Wb-SSM correspond to those found in turbulent flows, and assess its ability to reproduce key features of turbulence phenomenology.

### 4.1. Wavelet representation of multiscale velocity fields

We aim at generating an n-dimensional velocity field in a periodic box of length $L$. Consistently with the arguments presented above, we suggest the following form for the velocity field, replacing the Fourier decomposition with an MRA:

$$ \mathbf{u}(\mathbf{x}, t) = \sum_{j=J_{\min}}^{J_{\max}} \sum_{\mathbf{k}} \sum_{\varepsilon} d_{\mathrm{div}j, \mathbf{k}}^{\varepsilon}(t) \Psi_{\mathrm{div}j, \mathbf{k}}^{\varepsilon}(\mathbf{x}). \quad (18) $$

The velocity field is represented as a superposition of static, divergence-free wavelets $\Psi_{j,\mathbf{k}}^{\varepsilon}$. The time dependence enters through the coefficients $d_{j,\mathbf{k}}^{\varepsilon}(t)$, which will be modeled in the following sections. This formulation implies that the spatial structures of the flow remain fixed, while their amplitude and orientation evolve over time. Although this approach does not capture the sweeping effect of eddies, it preserves the rigorous structure of the multiresolution analysis.

Moreover, turbulence, being inherently multiscale yet bounded in scale range, requires the introduction of minimum and maximum levels of detail, $j_{\min}$ and $j_{\max}$, which define the largest and smallest flow scales, respectively. Within this range, each level of detail $j$ is associated with vortices of size $\ell_j = L/2^j$, where $L$ is the dimension of the box. Each vortex of level $j$ thus has a specific contribution to the global energy spectrum. This is illustrated in Fig. 3(a) and 3(b), where two examples of synthetic
<figure><img src="image_4.png" /><figcaption>Figure 3: Snapshots of a velocity component for 2D flows generated with a single level j = 3 (a) and j = 5 (b).</figcaption></figure>

fluids reconstructed with a single level j are given. Their corresponding energy spectra are plotted in Fig. 4(a), and the maximum of the energy is consistently found for a wave number of the order of magnitude $\kappa_j = 2\pi/\ell_j$. The energy associated with level j thus depends directly on the values of $d_{\text{div}j,\mathbf{k}}^{\varepsilon}$.

Homogeneous Isotropic Turbulence (HIT) is characterized by a statistically zero mean velocity field, while the variance of the velocity fluctuations defines the turbulent kinetic energy. This energy is distributed across scales according to the energy spectrum $E(\kappa)$, which may be represented by canonical models such as the Kolmogorov spectrum. Accordingly, we propose a Gaussian distribution of coefficients $d_{\text{div}j,\mathbf{k}}^{\varepsilon}$ with zero mean and variance depending on the associated level:

$$d_{\text{div}j,\mathbf{k}}^{\varepsilon} \sim \mathcal{N}(0, \sigma_j^2), \qquad (19)$$

where $\sigma_j$ is chosen to reproduce a given energy spectrum. By superposing wavelets of different sizes and positions, one obtains increasingly complex multiscale velocity fields that remain divergence-free, as illustrated in Fig. 5. In the next sections, we will exploit these properties to model wavelet coefficients $d_{\text{div}j,\mathbf{k}}^{\varepsilon}(t)$, so that the resulting velocity field mimics key statistical features of turbulence in both space and time.
<figure><img src="image_5.png" /><figcaption>Figure 4: (a) Energy Spectra cut beyond $3\kappa_j$ and (b) total energy budget for single-level Wb-SSM. Blue line stands for 2D field, red line for 3D field. Dashed and dotted lines are empirical laws given in Eq. (26).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: 2D velocity field generated with random coefficients $d_{\text{div}}^{\epsilon}$ in Eq. (17).</figcaption></figure>
## 4.2. Generating a frozen velocity field with a prescribed energy spectrum

To generate a frozen turbulent-like velocity field, we need to determine the variance $\sigma_j^2$ of the wavelet coefficients $d_{\text{div},j,\mathbf{k}}^{\varepsilon}$ at each level $j$ in order to reproduce a target energy spectrum $E(\kappa)$. Two approaches, one analytical and one empirical, are proposed to define this scale-dependent variance. In the sections that follow, we compare the ability of the Wb-SSM and Fb-KS to reproduce a target energy spectrum, and we highlight the limitations of the Wb-SSM for steep spectra.

### 4.2.1. Reproducing a model energy spectrum with wavelets

In the context of Fourier-based Kinematic Simulation of turbulence, we observed in Eq. (7) that a target energy spectrum can be obtained by distributing uniformly the energy $E(\kappa_n)\Delta\kappa_n$ among the $n_n$ Fourier modes spectrally localized at this wave number. We aim at extending this reasoning to the computation of wavelet coefficients. However, when using non-orthogonal wavelet bases, the total energy is not the sum of the individual contributions due to the absence of mutual orthogonality between basis functions:

$$ \sum_j E(\kappa_j) \Delta\kappa_j \neq \sum_j \sum_k \sum_\varepsilon \left( d_{\text{div},j,k}^{\varepsilon} \right)^2 . $$

This discrepancy arises from spectral interference between basis functions $\Psi_{\text{div},j,k}^{\varepsilon}$. However, if the coefficients $d_{\text{div},j,k}^{\varepsilon}$ are treated as independent, identically distributed (i.i.d.) random variables, one can recover an additive structure for the energy spectrum in a mathematical expectation sense.

Let us define the synthetic velocity field at scale $j$ as

$$ \mathbf{u}_j = \sum_k \sum_\varepsilon d_{\text{div},j,k}^{\varepsilon} \Psi_{\text{div},j,k}^{\varepsilon} . $$

Hence,

$$ \mathbf{u} = \sum_j \mathbf{u}_j \qquad (20) $$

As the wavelet coefficient are assumed to be independent random variables, in particular $(\mathbf{u}_j)_{j \in [0, J]}$ are independent variables. A straightforward calculation (see Appendix B) shows that for two independent random fields $X$ and $Y$, the turbulent energy spectrum
is additive in expectation:

$$ \mathbb{E}[E_{X+Y}(\kappa)] = \mathbb{E}[E_X(\kappa)] + \mathbb{E}[E_Y(\kappa)], \quad (21) $$

Thanks to this property, we obtain

$$ \mathbb{E}[E_{u}(\kappa)] = \sum_{j} \mathbb{E}[E_{u_{j}}(\kappa)]. \quad (22) $$

It is important to note that Eq. (22) is not a trivial result. Indeed, the energy spectrum is not a linear functional of the wavelet-decomposed velocity field, and individual wavelets possess non-localized (i.e., non-punctual) spectral support across wave number space. The validity of this additive property in expectation relies critically on the assumption of statistical independence between the wavelet coefficients across scales and positions.

Even though wavelets spectra are not completely localized, from Fig. 4(a), we observed that the spectral support of $E_{u_j}(\kappa)$ is concentrated around a central wave number $\kappa_j$. To fit a target energy spectrum $E_\odot(\kappa)$, this translates into the following condition:

$$ \mathbb{E}[E_{u_j}(\kappa_j)] \approx E_\odot(\kappa_j), \quad (23) $$

Assuming second-order stationarity and homogeneity, and letting $\mathbb{E}[(d_{\text{div } j, \mathbf{k}}^{\varepsilon})^2] = \sigma_j^2$ (independent of position $\mathbf{k}$), we write the expected spectrum at scale $j$ as

$$ \begin{aligned} \mathbb{E}[E_{u_j}(\kappa)] &= \sum_{\mathbf{k}} \sum_{\varepsilon} \mathbb{E} \left[ (d_{\text{div } j, \mathbf{k}}^{\varepsilon})^2 \right] E_{\Psi_{\text{div } j, \mathbf{k}}^{\varepsilon}}(\kappa) \\ &= \sigma_j^2 \sum_{\mathbf{k}} \sum_{\varepsilon} E_{\Psi_{\text{div } j, \mathbf{k}}^{\varepsilon}}(\kappa). \end{aligned} \quad (24) $$

Due to the wavelet basis normalization and approximate spectral localization, the energy of a single wavelet mode satisfies (cf. Parseval's identity):

$$ \frac{1}{2} \| \Psi_{\text{div } j, \mathbf{k}}^{\varepsilon} \|^2 = \int_{0}^{\infty} E_{\Psi_{\text{div } j, \mathbf{k}}^{\varepsilon}}(\kappa) \, d\kappa \approx E_{\Psi_{\text{div } j, \mathbf{k}}^{\varepsilon}}(\kappa_j) \Delta \kappa_j := C^{\varepsilon}. $$

Since the localization constant $C^\varepsilon$ depends only on the wavelet type $\varepsilon$, we define the total contribution per position as

$$ \sum_{\varepsilon} E_{\Psi_{\text{div } j, \mathbf{k}}^{\varepsilon}}(\kappa_j) \Delta \kappa_j := C. $$
Thus, the expected spectral energy per scale becomes

$$ \mathbb{E}[E_{\mathbf{u}_j}(\kappa_j)]\Delta\kappa_j = \sigma_j^2 \cdot \sum_{\mathbf{k}} C = \sigma_j^2 \cdot 2^{nj} \cdot C, $$

where n denotes the spatial dimension and $2^{nj}$ is the number of wavelets at scale j in a regular multiresolution grid.

Finally, matching this to the target spectrum $E_\odot(\kappa_j)$ as required by Eq. (23), we obtain the scaling law for the coefficient variance:

$$ \sigma_j^2 := \mathbb{E} \left[ \left( d_{\text{div } j, \mathbf{k}}^{\varepsilon} \right)^2 \right] = \frac{1}{C} \cdot \frac{E_\odot(\kappa_j) \Delta\kappa_j}{2^{nj}}. \quad (25) $$

The normalization constant C is ultimately determined to ensure that the total kinetic energy of the synthesised velocity field matches the target turbulent kinetic energy k. Eq. (25) thus provides a practical framework for generating divergence-free velocity fields with a prescribed energy spectrum with wavelets. The next section introduces an alternative way to recover this scaling.

### 4.2.2. An empirical law for the detail coefficients

The scaling law (25) can also be retrieved empirically, allowing the constant C to be determined for the current wavelet basis. To demonstrate this, Figure 4(b) reports the variance $\sigma_{j,\text{ref}}^2$ of single-level wavelet fields (i.e., $j_{\text{min}} = j_{\text{max}}$), for $d_{\text{div } j, \mathbf{k}}^{\varepsilon} = 1$. Based on these data, the following empirical laws are obtained for two- and three-dimensional flows:

$$ \begin{aligned} \sigma_{j,\text{ref}}^2 &= 0.24 \left( \frac{\ell_j}{\ell_{j_{\text{max}}}} \right)^{-2}, & \text{in 2D,} \\ \sigma_{j,\text{ref}}^2 &= 0.54 \left( \frac{\ell_j}{\ell_{j_{\text{max}}}} \right)^{-3.1}, & \text{in 3D.} \end{aligned} \quad (26) $$

where $l_j = L/2^j$. Note that these empirical laws are valid only for the present spline wavelets, and must be recalculated for a different choice of MRA.

Normalizing the target variance $E(\kappa_j)\Delta\kappa_j$ to be recovered at scale j by the reference variance $\sigma_{j,\text{ref}}^2$, we obtain nearly the same scaling as in (25) for $\sigma_j$:

$$ \sigma_j^2 = \frac{E(\kappa_j)\Delta\kappa_j}{\sigma_{j,\text{ref}}^2}. \quad (27) $$
In general this reasoning fails because of the potential overlapping of the wavelet basis function spectra. Yet, as shown previously, independence of the random wavelet coefficients makes it valid in the present model.

In 2D, we recover $n = 2$ as expected, whereas in 3D $n = 3.1$ provides a slightly more accurate fit. Although the difference is hardly visible, we retain the value 3.1 for the subsequent simulations.

In the following, we assess the effectiveness of the proposed spectral fitting strategy. We then introduce an evolution law that governs the temporal dynamics of the synthetic velocity field while preserving the prescribed energy spectrum. Finally, this dynamic extension is applied to the simulation of particle-laden turbulent flows.

### 4.2.3. Spectrum calibration: assessment and comparison with Fb-KS

To evaluate the ability of the Wb-SSM approach to reproduce prescribed energy spectra, we compare its performance against the Fb-KS method. As a reference, we consider the canonical Kolmogorov spectrum:

$$E(\kappa) = \sigma_u^2 \kappa_0^{2/3} \kappa^{-5/3}, \quad \forall \kappa \in [\kappa_{\min}; \kappa_{\max}] \qquad (28)$$

which is commonly used in the literature [39]. Figure 6 compares energy spectra obtained using wavelet and Fourier (Fb-KS) approaches. The energy spectrum obtained from wavelet fields closely matches the prescribed model, except for the extreme values of wave numbers. The wavelet model's degrees of freedom do not allow for independent control of each wave number but rather influence a range of wave numbers associated with a particular wavelet representation. Compared to the Fourier spectrum, the wavelet spectrum appears much smoother. Indeed, in a periodic box, Fourier modes are not continuously distributed in frequency; only discrete wave numbers compatible with the box size are allowed. This discrete sampling leads to oscillations and sparsity in the Fourier spectrum, as observed here. In contrast, the wavelet approach, with its compact support and overlapping frequency content, yields a smoother and more continuous energy spectrum. This difference highlights the inherent dissimilarity between Fourier and wavelet approaches: the Fourier method provides highly resolved frequency information but lacks local spatial details, whereas the wavelet method offers
<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison of energy spectra obtained with the Fb-KS (noisy blue curve) and Wb-SSM (smooth red curve), targeting a -5/3 Kolmogorov spectrum (black line).</figcaption></figure>

less detailed spectral resolution but provides better spatial information. Each wave-
length contributes to the spectrum at its corresponding discrete energy in the Fourier
representation. On the other hand, each wavelet already covers a continuous range of
wavelengths (Fig. 4(a)), and when superimposed, these wavelet spectra populate the
spectrum continuously.

4.2.4. Limitations of the approach for steep spectra

However, the present model does not allow for the reproduction of arbitrary target
spectra. This limitation arises from the assumption of statistical independence among
wavelet coefficients, combined with the limited regularity of the underlying wavelet
basis functions. Figure 7 shows that, in addition to the main energy contribution
centred around the wave number $\kappa_j$, the energy spectrum of a single wavelet exhibits
pronounced high-wave number bumps. The slow decay of these bumps reflects the
limited regularity of the wavelet basis. When a physical velocity field is decomposed
onto a wavelet basis, this phenomenon is not visible. In fact, the hierarchical structure
of the wavelet transform allows small-scale wavelets to compensate for the lack of
smoothness in large-scale components through cross-scale interactions.

In contrast, within the proposed synthesis framework, small-scale wavelet coeffi-
cients are generated independently of the large-scale structure and are prescribed solely
based on the target energy spectrum. As a result, cross-scale compensations of ir-
<figure><img src="image_8.png" /><figcaption>Figure 7: Energy spectrum of a 3D divergence-free wavelet ($j = 3, ε = (1, 0, 0)$).</figcaption></figure>

regularities (such as non-differentiability) are absent from the synthetic field. Since the wavelet coefficients are uncorrelated, these bumps may persist in the total energy spectrum of the synthesised velocity field.

Consequently, the ability to match a target energy spectrum depends on the relative spectral decay rates of the target and the wavelet basis. To highlight this, we can analyse the independent contributions of each scale of a synthesised velocity field $u_i$ to the energy spectrum. This is made possible thanks to Eq. (22). Specifically, the model can approximate a given spectrum only if its decay is less steep than that of the individual wavelet spectra. This explains why the model successfully captures a Kolmogorov $-5/3$ energy spectrum, as illustrated in Fig. 8. In this case, the analysis of scale-by-scale contributions of the wavelet levels reveals that the bumps are hidden beneath the overall spectral decay. However, the model fails to reproduce more stringent spectra shapes, such as the Passot-Pouquet spectrum depicted in Fig. 9. In the latter case, due to the narrow range of scales involved and the steep spectral decay, only a single level of wavelets significantly contributes to the total energy. Consequently, the inability of the model to couple higher level wavelet coefficients to the bumps—due to the independence assumption—prevents an accurate spectral reconstruction beyond a certain cut-off wave number.

To conclude this section, we emphasize that the ability of the wavelet-based model to recover Fourier-based statistical properties is not self-evident. Nonetheless, this raises a broader question regarding the appropriateness of striving to reproduce Fourier
<figure><img src="image_9.png" /><figcaption>Figure 8: Scale-by-scale energy spectra for a Kolmogorov-type inertial range spectrum: $E(\kappa) = \sigma_u^2 \kappa_0^{2/3} \kappa^{-5/3}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Scale-by-scale energy spectra for the Passot-Pouquet model spectrum: $E(\kappa) = \frac{16\sqrt{2}\sigma_u^2}{\sqrt{\pi}\kappa_0} \left(\frac{\kappa}{\kappa_0}\right)^4 \exp\left(-2\left(\frac{\kappa}{\kappa_0}\right)^2\right)$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Snapshot of a 2D velocity field obtained for different values of $j_{max} \in \{3, 5, 7\}$ and fixed $j_{min} = 0$ in a cubic box of size $L = 1$.</figcaption></figure>

statistics using a wavelet framework. Given that compactly supported wavelet modes more closely resemble localized turbulent structures than the globally supported Fourier modes, one may question whether spectral fidelity in the Fourier sense is the most relevant criterion for evaluating the realism of such synthetic turbulence fields.

### 4.3. Characteristic length and timescales

Now that we are able to generate a frozen turbulent-like velocity field, we can investigate how the parameters $j_{min}$ and $j_{max}$ relate to the characteristic length and timescales of turbulence. For the spectrum model given in Eq. (28), these wave numbers are linked respectively to the integral scales $L = 2\pi/\kappa_{min}$ and Kolmogorov scales $\eta \sim 2\pi/\kappa_{max}$. This is illustrated in Figure 10. It represents one velocity component of the synthetic flow described by Eq. (18), with different values of $\ell_{max}$. As expected, the lower $\ell_{max}$, the more details are present in the simulation, revealing the presence of small and very localized vortices.

We can then verify the scaling of the timescales with the model's parameters. In Fig. 11, we study the dependency of the Kolmogorov timescale with the two parameters, $\ell_{min}$ and $\ell_{max}$, which are respectively the smallest and the largest length scales in the flow. We obtain scaling in power laws:

$$ \tau_{\eta} = \sqrt{\nu/(\langle \epsilon \rangle)} \sim \sigma_{u}^{-1} \ell_{min}^{0.4} \ell_{max}^{0.6}, \quad (29) $$

which are close to the theoretical ones, predicted by the Kolmogorov phenomenology [73]:

$$ \tau_{\eta} \sim \sigma_{u}^{-1} L^{1/3} \eta^{2/3}. \quad (30) $$
<figure><img src="image_12.png" /><figcaption>Figure 11: Kolmogorov timescale $\tau_{\eta}$ measured for different values of $L \sim \ell_{min}$ and compared with a scaling law of power 0.4 (a) and different values of $\eta \sim \ell_{max}$, compared with a scaling law of power 0.6.</figcaption></figure>

We recall that the largest timescale in the flow is the autocorrelation time of the fluid particle velocities, i.e., the Lagrangian integral timescale defined by:

$$T_L = \frac{1}{\sigma_u^2} \int_0^{+\infty} R_L^u(\tau) \, d\tau \qquad (31)$$

Since $T_L \sim \sigma_u^{-1} L$, it is also expected to scale with the largest scales in the flow. Figure 12 shows the autocorrelation function of the velocity for different values of $\ell_{min}$ and $\ell_{max}$. The overlap of dashed and full lines confirms the independence of the autocorrelation, hence $T_L$, with the small scales. On the other hand, the parameter $j_{min}$ has a clear impact on the characteristic time of autocorrelation $T_L$.

In summary, the parameters $j_{min}$ and $j_{max}$ effectively control the range of length and timescales represented in the synthetic turbulent field. While $j_{max}$ governs the smallest resolved scales and the associated Kolmogorov timescale $\tau_{\eta}$, $j_{min}$ sets the largest flow structures and the Lagrangian integral timescale $T_L$. These results confirm that the wavelet-based construction captures the expected multiscale characteristics of turbulence, establishing a direct link between the wavelet parameters and the physical scales of the flow.
<figure><img src="image_13.png" /><figcaption>Figure 12: Lagrangian velocity autocorrelation for fluid particles, shown for different values of $j_{min}$ and $j_{max}$ (solid lines: $j_{max}$ = 5; dashed lines: $j_{max}$ = 6).</figcaption></figure>

### 4.4. Evolution of the wavelet field using stochastic processes

To mimic the Navier-Stokes equations, it is required to make the flow field evolve in time, and for that purpose, we propose to write a time evolution law directly on the coefficients. Hausmann et al. [50] proposed to use a Galerkin projection of the linearise Navier-Stokes equation. To make this approach tractable, they relied on several assumptions. First, the solution is a locally optimal solution instead of a globally optimal solution. Second, there is a spectrum enforcement based on the modulation of the transport and diffusive part. Unfortunately, the resulting approach is not applicable to a case with no large-scale resolution, as it is not enabling positive energy forcing in this situation. Therefore, we propose a simpler approach based on stochastic processes on the details of the wavelet reconstruction. The constraints on the evolution model are such that the stationarity of the flow must be preserved, as well as the energy spectrum. Thus, at any time, the distribution of the coefficients must remain a Gaussian of variance $\sigma_j$. We propose a stochastic Langevin model for each coefficient, which guarantees the two properties mentioned above, and introduces an additional control parameter: the eddy turnover time of that level $T_{E,j}$:

$$d(d_{\text{div},j,k_1,k_2}^e) = -\frac{1}{T_{E,j}} d_{\text{div},j,k_1,k_2}^e dt + \sqrt{\frac{2\sigma_j^2}{T_{E,j}}} dW_t. \quad (32)$$
<figure><img src="image_14.png" /><figcaption>Figure 13: Lagrangian velocity autocorrelation, for $j_{\text{min}} = 1$ (red to yellow colours) and $j_{\text{min}} = 3$ (blue to green colours), with fixed $j_{\text{max}} = 6$. From dark to light shades, increasing values of $\lambda \in \{0.01, 0.05, 0.1, 0.5, 1, 2, 5\}$.</figcaption></figure>

The random variables are Gaussian and driven by Ornstein-Uhlenbeck processes. They are stationary by construction with variance $\sigma_j$ and timescales $T_{E,j}$. We introduce the eddy turnover frequency $\omega_j = T_{E,j}^{-1}$ for which some models are already proposed in the literature in the framework of Fb-KS, and in the following, we use the model of Eq. (8).

Figure 13 shows the velocity autocorrelation for different values of $\lambda$. We remark that for very small values of $\lambda$ (up to frozen turbulence), the velocity can be negatively correlated, and a brief oscillating behaviour can occur before reaching the total decorrelation. For larger values of $\lambda$, the decrease of the autocorrelation function is monotonic and more rapid, and the Lagrangian integral timescale $T_L$ is hence reduced. This is expected because the dynamics introduced in the temporal equation decorrelates fluid particles velocities faster.

In conclusion, the Wb-SSM described by Eq. (18) ensures the spatial correlation of the flow. Unlike the Fb-KS, the velocity field is constructed from compactly supported modes. Additionally, the use of divergence-free wavelets guarantees that the resulting field is divergence-free. The above procedure enables the temporal evolution of a velocity field consistent with a prescribed energy spectrum. These features respectively satisfy points (i), (ii), and (iii) of the specifications outlined in Sec. 2.3. The remaining
points, (iv) and (v), will be assessed in the results section.

## 5. Assessment of the wavelet-based model on particle-laden HIT turbulence

In this section, we validate the Wb-SSM by verifying its ability to reproduce statistics of a three-dimensional mildly turbulent flow obtained by DNS with and without particles. The reference fields are simulated on a 128³ mesh, using the DNS Fieldz code [74]. This code solves the Navier-Stokes equations in a three-dimensional box with periodic boundary conditions by applying the divergence-free constraint. A pseudo-spectral approach is used, the time integration is done by the second-order Adams-Bashford method, and the forcing of the flow is ensured with the model of Chen et al. [75] by forcing only the smallest wave numbers to maintain the stationarity of the flow. Snapshots and energy spectra are observed once a stationary state is reached, the turbulence characteristic being summarized in Table 1.

<table><caption>Table 1: Properties for the 3D stationary turbulence obtained with FieldZ using a 128³ mesh.</caption><thead><tr><th>Re<sub>λ</sub></th><th>η</th><th>τ<sub>η</sub></th><th>ε</th><th>Δx</th><th>ν</th></tr></thead><tbody><tr><td>54</td><td>0.047</td><td>0.171</td><td>0.5</td><td>0.049</td><td>0.005</td></tr></tbody></table>

### 5.1. Single-phase turbulence

First, we adapt the variance $\sigma_j$, following Eq. (27) in order to reproduce the reference energy spectrum as closely as possible, see Fig. 5.1. We manage to reproduce a very large part of the spectral scales, except the smallest scales, at the limit of the mesh resolution, for which the wavelet field does not reproduce the decrease of the energy content of the DNS field. As already highlighted in paragraph 4.2.4, these bumps are most likely attributable to the low regularity of the B-splines used in constructing the wavelet basis.

Snapshots of velocity fields obtained in DNS (Fig. 14(a)), with the Wb-SSM (Fig. 14(b)) and with a Fb-KS (Fig. 14(c)) can be compared, and we observe structures of similar size and intensity between the DNS and Wb-SSM. The patterns appearing in
<figure><img src="image_15.png" /><figcaption>Figure 14: Instantaneous cut planes of the velocity magnitude for DNS (a), Wb-SSM (b) and Fb-KS (c).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: Comparison of energy spectra between DNS (black full line), Wb-SSM (red dashed line), and Fb-KS (blue dot-dashed line).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: Autocorrelation function of fluid particle velocity for DNS (black line), Fb-KS (blue lines), and Wb-SSM (red lines) with $\lambda = 0$ (solid lines) and with $\lambda = 0.5$ (dotted lines).</figcaption></figure>

the Fb-KS are very stretched, and this anisotropy probably comes from discrete forcing at certain wave numbers to impose periodicity. Isotropic wavelets [41] allow, on the contrary, to avoid the existence of elongated structures, which do not appear visibly in the DNS reference for HIT.

We now compare the temporal evolution of stationary HIT obtained with the DNS. The temporal evolution is controlled by the parameter $\lambda$ (Eq. (8)), which drives the eddy turn-over time. We therefore investigate its impact with comparison to DNS. Figure 16 compares the autocorrelation functions of DNS, wavelet, and Fb-KS for a frozen and a dynamic flow. Both methods retrieve the exponential behaviour of the DNS autocorrelation for $\lambda = 0$, i.e., with a frozen field. We thus conclude that ensuring the faithful reproduction of the energy spectrum, i.e., the spatial scales, is sufficient to reproduce the small and large temporal character scales (namely $\tau_{\eta}$ and $T_L$). However, $\lambda = 0$ yields a field with non-physical steadiness. For $\lambda = 0.5$, the Fb-KS produces unphysical oscillations of the autocorrelation, as already reported by Murray et al. [39]. The Wb-SSM still presents a physical monotonically-decreasing behaviour, but with a faster decorrelation than for the DNS, thus resulting in a reduced Lagrangian integral timescale as observed in Sec. 4.4. This shows that a global parameter for eddy turnover time is not sufficient to fit all turbulence statistics, thus suggesting a scale-
dependent parameter. This adaptation is straightforward in the Wb-SSM, for which we
can access each scale. However, in the scope of the present paper, we limit ourselves to
a global parameter, which will be shown to be sufficient to recover particle preferential
concentration—the main objective of the present work, as developed in Sec. 5.2.

5.2. Inertial particles in a turbulent flow

This section first aims at presenting and qualitatively characterizing the PPC phe-
nomena that occur for different particle Stokes numbers in a Wb-SSM (Sec. 5.2.1). We
then show in Sec. 5.2.2 the sensitivity of the model to the parameter $\lambda$, which controls
the eddy turn-over time. Finally, we seek to faithfully reproduce the particle segregation
obtained by DNS in Sec. 5.2.3 and compare the results with Fb-KS.

To quantify the PPC, we use the ratio $\langle n(x,t)^2 \rangle / \langle n(x,t) \rangle^2$ [76], which measures the deviation from homogeneity. Here, $\langle \cdot \rangle$ denotes a spatial average over the entire domain. This indicator approaches unity when particles are uniformly distributed and numerous.

5.2.1. *A multiscale model able to reproduce qualitatively PPC*

Wb-SSM enables to recover qualitatively the behaviour of inertial particles in a HIT flow. Figure 17 shows particle number density fields obtained with Wb-SSM populated by increasing inertial particles in frozen turbulence. The heterogeneity of the fields is typical of PPC. This phenomenon has already been well described in the literature [77], and our simulations are consistent with such observations. Inertial particles segregate when their characteristic time $\tau_p$ is of the same order of magnitude as the proper characteristic time of vortices, and more specifically with Kolmogorov vortices of characteristic time $\tau_\eta$ (i.e. $\text{St}_\eta = 1$).

5.2.2. *Influence of the eddy turnover velocity*

In dynamic turbulence, vortices are created, transformed, deformed, and dissipated,
and their characteristic time is much more challenging to estimate. However, we know
that it is necessarily reduced compared to that in frozen turbulence. The particles in a
vortex do not necessarily have the time to reach their maximum segregation, the vortex
having changed in the meantime. Figure 18(a) compares the time evolution of PPC in a
<figure><img src="image_18.png" /><figcaption>Figure 17: Particle number density fields obtained at $t_f = 6$ with particles of different Stokes numbers transported by Wb-SSM. Parameters of the Wb-SSM are $\lambda = 0$, $j_{min} = 1$, $j_{max} = 5$. Eulerian fields are obtained with box projection on a 2$^7$ grid.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18: (a) Temporal evolution of PPC for $\lambda = 0$ (full lines) and $\lambda = 0.5$ (dotted lines), with logarithmic temporal scale. (b) Influence of $\lambda$ on the PPC as a function of the Stokes number $St_\eta$. PPC is measured at time $30T_L$, with a grid resolution of $2^6$.</figcaption></figure>

frozen ($\lambda = 0$) and in a dynamic ($\lambda = 0.5$) Wb-SSM. We observe that the disperse phase reaches a stationary regime for PPC when $\lambda \neq 0$. We expect the stationary regime to be reached for $t > \max(T_{E,j})$, the maximum eddy turnover time of Wb-SSM. On the other hand, with $\lambda = 0$, PPC keeps increasing for particles such that $St_L < 1^2$. Figure 18(b) shows the segregation for different values of $\lambda$, measured at time $30T_L$. We can see that $\lambda$ controls the intensity of PPC. As expected, increasing the eddy turnover time reduces the level of particle segregation.

The model parameters thus impact the two-point statistics of inertial particles consistently with the scale and time changes they represent in the fluid. We can therefore verify that with realistic parameters obtained by reproducing the statistics of a turbulent fluid flow, the behaviour of inertial particles in synthetic flow is correctly predicted.

### 5.2.3. Particles in stationary HIT

Finally, we now seek to validate the three-dimensional Wb-SSM with time evolution, which ensures stationarity of the disperse phase statistics. Energy spectra, snapshots

<footnote>$^2$Particles at St = 70 experiences a small bump because of the time required for the particle to accelerate in the flow, before reaching the statistically-stationary quasi-Brownian limit.</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 19: Evolution of particle preferential concentration with Stokes number obtained with 3D DNS (black full line), Wb-SSM (red dashed line) and Fb-KS (blue dot-dashed line).</figcaption></figure>

of the velocity field, and autocorrelation of fluid particles have already been compared with Fb-KS and Wb-SSM, and the latter shows better agreement with the DNS of HIT.

We plot in Fig. 19 the PPC according to the particles’ inertia and compare the results obtained in DNS with stationary wavelet and Fb-KS<fnref n="3" />. Results are shown for $\lambda = 0.5$, representing the best DNS agreement. On the one side, the segregation obtained with Fb-KS is shifted towards larger values of Stokes numbers. This shift has already been observed in [39]. This shows that small scales probably need to be better reproduced. On the other side, Wb-SSM not only predicts the correct peak of segregation but also correctly reproduces the preferential concentration for all $St_\eta > 1$. The PPC is under-predicted for small Stokes numbers, which might be related to the rapid decorrelation of small scales observed in Sec. 4.4. The eddy turnover of small scales is probably too high, and this excessively reduces the segregation at those specific scales. Further studies should consider other models of $\omega_j$, with a weaker increase in $\kappa_j$ to avoid the rapid forcing of the small scales eddy turnover time. We can see that in this case, there

<footnote n="3">³A limited number of particles was used to perform this comparison, hence the observed offset in the PPC metric. Nevertheless, this limitation is not critical, as the results still demonstrate that the Wb-SSM approach can accurately reproduce the DNS statistics.</footnote>
is no scale separation because the peak is reached at $St_η = 1$, as a consequence of the small Reynolds number of the turbulence.

The collective behaviour of the particles, measured through particle segregation and trajectory crossing, is better predicted qualitatively and quantitatively by Wb-SSM than with Fb-KS. In addition, recall that in the dynamic case, the field obtained with a Fb-KS produces non-physical oscillating Lagrangian statistics.

## 6. Conclusion

In this work, we have designed a new modelling strategy for disperse phase flows which significantly improves the prediction of the particle preferential concentration. The proposed method uses a novel synthetic turbulence model based on divergence-free wavelets and a stochastic Langevin equation on the wavelet coefficients for temporal evolution. This formalism has the same advantages as a Fourier-based KS, namely its analytic formulation and numerical simplicity, the total control of the associated probability space, and its similarity with the Kolmogorov picture of turbulence. However, instead of sinusoidal functions used in the Fourier decomposition, wavelets have compact support and are more adapted to represent other characteristics of real turbulent flows. In particular, we have shown that the model shows improved performance over the Fourier-based approach in capturing important DNS features, such as time autocorrelation. Application to turbulent disperse phase flows shows an excellent agreement regarding particle preferential concentration compared to the DNS, resulting in a great improvement compared to the Fourier-based KS.

Whereas very good results are obtained for the PPC, there is room for improvements in the prediction of single-phase statistics. As mentioned in this work, the main limitation is in the parameter $\lambda$ which controls the PPC as well as the eddy turn over time. In fact, we observed that the parameter value required to match the DNS data for the PPC metric is not compatible with reproducing the DNS time autocorrelation. This parameter could also be set as scale-depend to add another degree of freedom to improve the two-time statistics.
The question of the validity of the wavelet model itself as a single-phase LES model
can be raised. Many properties on the synthetic velocity field remain to be observed,
and analytical quantities such as the dissipation or the turbulent kinetic energy could
be derived. Furthermore, conducting a comparative study with existing functional and
structural LES models (e.g., Smagorinsky) regarding fidelity and numerical cost would
be interesting. Ultimately, applying the model in a full LES simulation constitutes the
next step toward rigorously demonstrating its versatility. To make Wb-SSM compu-
tationally competitive with Fb-KS as a subgrid-scale model, it is essential to exploit
wavelet compression techniques to reduce the dimensionality of the synthetic velocity
field [78].

Using anisotropic wavelets is promising for flows with more singular geometry or
physics—such as shear layers or shocks—and a preferential direction can be better
represented using these wavelets. In particular, the compact support of our basis
functions offers a significant advantage over the Fb-KS approach in handling boundary
conditions. Extending this framework with anisotropic wavelets presents a compelling
direction for addressing the challenges of anisotropic boundary layer turbulence.

While the use of Fb-KS is limited to non-spatially-intermittent flows due to the
infinite support of Fourier modes, the compact support of Wb-SSM modes offers the
possibility to introduce spatial non-gaussianity in the model. This could be achieved by
modulating the wavelet coefficients with a multiplicative cascade process, as already
suggested by different authors [46–48].

Finally, regarding two-phase flows, Wb-SSM could be extended to take into account
two-way coupling, as the forcing term introduced by the presence of particles is easily
reproducible locally and spectrally with wavelets.

Acknowledgements

This work has been supported by a PhD grant from Region Ile-de-France DIM MATHINNOV (R. Letournel). The authors would like to thank Erwan Deriaz for fruitful discussions. This work has also been supported by the Agence Innovation Défense (AID), the CIEDS project OPEN-NUM-DEF (PI. L. Gouarin, M. Massot, T.
Pichard), a PhD Grant from the Ministry of Research at CentraleSupélec (C. Morhain), as well as the HPC@Maths Initiative (PI. L. Gouarin and M. Massot) of the Fondation de l'École polytechnique, which we would like to acknowledge. This work was performed using computational resources from the "Mésocentre" computing centre of Université Paris-Saclay, CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France (https://mesocentre.universite-paris-saclay.fr/). The authors would also like to thank Remi Zamansky for sharing the DNS code FieldZ, and Ludovic Goudenège for fruitfull discussions on stochastic modelling.

# Appendix A. Multi-Resolution Analysis

In this appendix, we briefly review the construction of divergence-free wavelets, using the example of biorthogonal splines.

## Appendix A.1. Biorthogonal splines

Let us consider the example of biorthogonal splines functions [42, 62]. Their scaling functions are B-splines, which are defined piecewise by polynomials, and the associated wavelets can be obtained by linear and translated combinations of B-splines. For example, B-splines of degree 1 and 2, $\varphi_0$ and $\varphi_1$, associated to wavelets $\psi_0$ and $\psi_1$ are represented in Fig. 20(a) and 20(b) while Fig. 20(c) illustrates the dilatation and translation of $\psi_1$. Four possible translations of the wavelet at level $j = 0$ are plotted in dotted lines for a periodic function in $[0, 1]$. At level $j = 1$, the support of the wavelet (in dashed lines) is half the size of the mother wavelet, and its amplitude is twice as large to preserve its norm. The same goes for the higher levels, and we notice that the information is more and more localized in space, and the generated field is less regular.

## Appendix A.2. Multidimensional wavelets

One can extend the notion of multi-resolution analysis to spaces of any dimension. The simplest way is to perform tensor products between wavelets and scaling functions of the same level (isotropic wavelets). In the following, we consider the case of two-dimensional spaces, which can be directly generalized to higher dimensions.
<figure><img src="image_21.png" /><figcaption>Figure A.20: Scaling functions (blue) and wavelets (red) for splines of (a) degree 1 and (b) degree 2. (c) gives an example of translated and dilated wavelets of level 0 (dotted lines), 1 (dashed lines), 2 (full line).</figcaption></figure>
We consider two MRA of $L^2(\mathbb{R})$: $(V_J^0)$ and $(V_J^1)$ respectively associated to the scale and wavelets couple-functions $(\varphi_0, \psi_0)$ and $(\varphi_1, \psi_1)$. The MRA space tensor product $V_J^0 \otimes V_J^1$ is generated by the scaling functions $\{\varphi_{0j,k_1}(x)\varphi_{1j,k_2}(y); (k_1, k_2) \in \mathbb{Z}^2\}$, where $\varphi_{0j,k_1}(x) = 2^{j/2}\varphi_0(2^j x - k_1)$ and the same goes for $\varphi_{1j,k_2}$. The first coordinate is decomposed in the MRA $(V_J^0)$, and the second one into $(V_J^1)$. Then any function $f_J$ of $V_J^0 \otimes V_J^1$ can be written as:

$$f_J(x, y) = \sum_{k_1, k_2 \in \mathbb{Z}} c_{J, k_1, k_2} 2^J \varphi_0(2^J x - k_1) \varphi_1(2^J y - k_2). \quad (A.1)$$

Iterating the decomposition of $V_{J+1}$ using the supplementary spaces, we have:

$$ \begin{aligned} V_{J+1}^0 \otimes V_{J+1}^1 &= (V_J^0 \oplus W_J^0) \otimes (V_J^1 \oplus W_J^1) \\ &= (V_J^0 \otimes V_J^1) \oplus (W_J^0 \otimes V_J^1) \oplus (V_J^0 \otimes W_J^1) \oplus (W_J^0 \otimes W_J^1) \\ &= (V_0^0 \otimes V_0^1) \oplus (W_0^0 \otimes V_0^1) \oplus (V_0^0 \otimes W_0^1) \oplus (W_0^0 \otimes W_0^1) \\ &\qquad \oplus \dots \\ &\qquad \oplus (W_J^0 \otimes V_J^1) \oplus (V_J^0 \otimes W_J^1) \oplus (W_J^0 \otimes W_J^1), \end{aligned} \quad (A.2) $$

which leads to the following decomposition of $f_J$, introducing the daughter wavelets $\psi_{1j,k}$ and $\psi_{0j,k}$:

$$ f_J(x, y) = \sum_{(k_1, k_2) \in \mathbb{Z}^2} c_{k_1, k_2} \varphi_0(x - k_1) \varphi_1(y - k_2) \\ + \sum_{J=0}^{J-1} \left( \sum_{k_1, k_2} d_{j, k_1, k_2}^{(1,0)} \psi_{0j, k_1}(x) \varphi_{1j, k_2}(y) + \sum_{k_1, k_2} d_{j, k_1, k_2}^{(0,1)} \varphi_{0j, k_1}(x) \psi_{1j, k_2}(y) \right. \\ \left. + \sum_{k_1, k_2} d_{j, k_1, k_2}^{(1,1)} \psi_{0j, k_1}(x) \psi_{1j, k_2}(y) \right) \quad (A.3) $$

This decomposition involves three types of wavelets:

* one along the x direction: $\psi^{(1,0)}(x, y) = \psi_0(x)\varphi_1(y)$. Coefficients in the decomposition associated with such terms are called horizontal detail coefficients because the wavelet function is applied in the horizontal direction;

* one along the y direction: $\psi^{(0,1)}(x, y) = \varphi_0(x)\psi_1(y)$, respectively with vertical detail coefficients;
• one along both directions at once: $\psi^{(1,1)}(x, y) = \psi_0(x)\psi_1(y)$, associated with diagonal coefficients;

and one scaling function: $\varphi^{(0,0)}(x, y) = \varphi_0(x)\varphi_1(y)$.

The interest of this wavelet basis lies in the fact that their support is proportional to $2^{-j}$ in each direction and thus preserves a spatial and frequency localization, level by level, in all directions at once. This is why this decomposition is called isotropic.

## Appendix A.3. Divergence-free wavelets

Introduced by Lemarié-Rieusset [40], divergence-free wavelets aim at providing a multi-resolution analysis for the following space:

$$H_{\text{div},0}(R^n) = \{f \in (L^2(\mathbb{R}^n))^n / \text{div} f \in L^2(\mathbb{R}^n), \text{div} f = 0\}. \quad (A.4)$$

**Proposition Appendix A.1.** Let $(V_j^1)$ be a multi-resolution analysis associated to the scaling function $\varphi_1$, such that $\varphi_1 \in C^{1+\epsilon}$ for any $\epsilon > 0$ and a wavelet $\psi_1$. Then there exists an MRA $(V_j^0)$ associated to a scaling function $\varphi_0$ and a wavelet function $\psi_0$ verifying:

$$\varphi_1'(x) = \varphi_0(x) - \varphi_0(x-1), \quad \psi_1'(x) = 4\psi_0(x). \quad (A.5)$$

The isotropic divergence-free wavelets in $\mathbb{R}^n$ are obtained by tensor products of functions $\phi_0, \psi_0, \phi_1$ and $\psi_1$ in the following manner:

**Proposition Appendix A.2.** There exist $(n-1)(2^n-1)$ vector functions $\psi_{\text{divi}}^{\epsilon} \in H_{\text{div},0}(R^n)$, with $i \in \mathbb{I}_{\epsilon} = [[1; n-1]]$ and $\epsilon \in \Omega_n^* = \{0, 1\}^n, \setminus \{0, \dots, 0\}$, compactly supported such that every vector function $u \in H_{\text{div},0}(R^n)$ can be uniquely expanded:

$$u = \sum_{i \in \mathbb{I}_{\epsilon}} \sum_{j \in \mathbb{Z}} \sum_{\mathbf{k} \in \mathbb{Z}^n} \sum_{\epsilon \in \Omega_n^*} d_{\text{divi},j,\mathbf{k}}^{\epsilon} \Psi_{\text{divi},j,\mathbf{k}}^{\epsilon} \quad (A.6)$$

with $\Psi_{\text{divi},j,\mathbf{k}}^{\epsilon}(x) = 2^{nj/2} \Psi_{\text{divi}}^{\epsilon}(2^j\mathbf{x} - \mathbf{k})$.

The generating wavelets $\Psi_{\text{divi}}^{\epsilon}$ take the general following form in $\mathbb{R}^n$:
**Definition Appendix A.3.** Let $\epsilon \in \Omega_n^* = \{0, 1\}^n, \setminus \{0, ..., 0\}$ be given. Let $i_0 \in \{1, ..., n\}$ such that $\epsilon_{i_0} = 1$. Then for every $i \in \{1, ..., n\} \setminus \{i_0\}$, the vector function $\Psi_{\text{divi}}^{\epsilon} = [\Psi_{\text{divi}}^{\epsilon}]_l$ with $l = 1, n$ writes:

$$ [\Psi_{\text{divi}}^{\epsilon}]_l = \begin{cases} 0 & \text{if } l \notin \{i, i_0\} \\ \theta_{0, \epsilon_1} \otimes \dots \otimes \theta_{0, \epsilon_{i-1}} \otimes \theta_{1, \epsilon_i} \otimes \theta_{0, \epsilon_{i+1}} \otimes \dots \otimes \theta_{0, \epsilon_n} & \text{if } l = i \\ \eta_1 \otimes \dots \otimes \eta_n & \text{if } l = i_0 \end{cases} \quad (A.7) $$

where:

$$ \left\{ \begin{array}{ll} \eta_j = \theta_{0, \epsilon_j} & \text{for } j \neq i, i_0 \\ \eta_{i_0} = \psi_1/4 & \\ \eta_i = -\theta'_{1, \epsilon_i} & \end{array} \right. $$

with the notation for $r \in \{0, 1\}$:

$$ \theta_{r, \epsilon_j} = \begin{cases} \phi_r & \text{if } \epsilon_j = 0 \\ \psi_r & \text{if } \epsilon_j = 1 \end{cases} $$

For example, the following two-dimensional scaling function can be defined:

$$ \Phi_{\text{div}}(x, y) = \begin{pmatrix} \varphi_1(x)\varphi'_1(y) \\ -\varphi'_1(x)\varphi_1(y) \end{pmatrix} = \begin{pmatrix} \varphi_1(x)[\varphi_0(y) - \varphi_0(y-1)] \\ -[\varphi_0(x) - \varphi_0(x-1)]\varphi_1(y) \end{pmatrix}, \quad (A.8) $$

along with the divergence-free vector wavelets, given by:

$$ \begin{align} \Psi_{\text{div}}^{(1,0)}(x,y) &= -\frac{1}{4} \begin{pmatrix} \psi_1(x)\varphi'_1(y) \\ -\psi'_1(x)\varphi_1(y) \end{pmatrix} = -\frac{1}{4} \begin{pmatrix} \psi_1(x)[\varphi_0(y) - \varphi_0(y-1)] \\ -4\psi_0(x)\varphi_1(y) \end{pmatrix}, \nonumber \\ \Psi_{\text{div}}^{(0,1)}(x,y) &= \frac{1}{4} \begin{pmatrix} \varphi_1(x)\psi'_1(y) \\ -\varphi'_1(x)\psi_1(y) \end{pmatrix} = \frac{1}{4} \begin{pmatrix} 4\varphi_1(x)\psi_0(y) \\ -[\varphi_0(x) - \varphi_0(x-1)]\psi_1(y) \end{pmatrix}, \nonumber \\ \Psi_{\text{div}}^{(1,1)}(x,y) &= \frac{1}{4} \begin{pmatrix} \psi_1(x)\psi'_1(y) \\ -\psi'_1(x)\psi_1(y) \end{pmatrix} = \frac{1}{4} \begin{pmatrix} 4\psi_1(x)\psi_0(y) \\ -4\psi_0(x)\psi_1(y) \end{pmatrix}. \tag{A.9} \end{align} $$

In three dimension, the decomposition of $(L^2(\mathbb{R}^3))^3$ is given by:

$$ (V_j^1 \otimes V_j^0 \otimes V_j^0) \times (V_j^0 \otimes V_j^1 \otimes V_j^0) \times (V_j^0 \otimes V_j^0 \otimes V_j^1). \quad (A.10) $$

There are $2 \times (2^3 - 1) = 14$ generating divergence-free wavelets $\Psi_{\text{div}}^{\epsilon}$, with the following type:

$$ \epsilon \in \{(0, 0, 1), (0, 1, 0), (1, 0, 0), (0, 1, 1), (1, 1, 0), (1, 0, 1), (1, 1, 1)\}. $$
The construction of these wavelets is not detailed here and the reader is referred to the appendices of the thesis of [42] for their implementation.

## Appendix B. Additivity of the energy spectrum for independent random fields

### Appendix B.1. Energy spectrum of homogeneous isotropic turbulence

For a homogeneous isotropic velocity field $\mathbf{u}$ with Fourier transform $\hat{\mathbf{u}}(\boldsymbol{\kappa})$, the energy spectrum is the angle-integrated power spectral density:

$$ \mathbb{E}[E(\boldsymbol{\kappa})] = \int_{\|\boldsymbol{\kappa}\|=\kappa} \mathbb{E}\left[\frac{1}{2}\|\hat{\mathbf{u}}(\boldsymbol{\kappa})\|^2\right] d\Omega_{\boldsymbol{\kappa}} \quad (B.1) $$

$$ = S_d(\kappa) \mathbb{E}\left[\frac{1}{2}\|\hat{\mathbf{u}}(\kappa)\|^2\right], \quad (B.2) $$

where $\kappa = \|\boldsymbol{\kappa}\|$. Here, $S_d(\kappa)$ is the surface area of a ($d-1$)-sphere of radius $\kappa$:

$$ S_d(\kappa) = \begin{cases} 2\pi\kappa, & \text{in 2D,} \\ 4\pi\kappa^2, & \text{in 3D.} \end{cases} \quad (B.3) $$

Hence, the energy spectrum is proportional to the squared modulus of the Fourier transform of the velocity field.

### Appendix B.2. Additivity property

Let $X(\mathbf{x})$ and $Y(\mathbf{x})$ be two independent random fields. The Fourier transform of their sum is given by

$$ \mathcal{F}[X+Y] = \mathcal{F}[X] + \mathcal{F}[Y]. \quad (B.4) $$

Taking the squared modulus of both sides, we obtain

$$ |\mathcal{F}[X+Y]|^2 = |\mathcal{F}[X]|^2 + |\mathcal{F}[Y]|^2 + 2\Re(\mathcal{F}[X] \overline{\mathcal{F}[Y]}). \quad (B.5) $$

Since $X$ and $Y$ are independent, their Fourier coefficients are uncorrelated. Assuming zero-mean fields, this implies

$$ \mathbb{E}[\mathcal{F}[X] \overline{\mathcal{F}[Y]}] = \mathbb{E}[\mathcal{F}[X]] \mathbb{E}[\overline{\mathcal{F}[Y]}] = 0. \quad (B.6) $$

Taking the expectation of the previous expression, the cross term vanishes, leading to

$$ \mathbb{E}[|\mathcal{F}[X+Y]|^2] = \mathbb{E}[|\mathcal{F}[X]|^2] + \mathbb{E}[|\mathcal{F}[Y]|^2]. \quad (B.7) $$
By definition, the turbulent energy spectrum $E(\kappa)$ is proportional to the expectation of the squared modulus of the Fourier transform (see Appendix B.1). Hence

$$ \mathbb{E}[E_{X+Y}(\kappa)] = \mathbb{E}[E_X(\kappa)] + \mathbb{E}[E_Y(\kappa)]. \qquad (B.8) $$

Thus, the energy spectrum of the sum of two independent random fields is equal to the sum of their individual energy spectra in expectation.

## References

[1] R. Letournel, F. Laurent, M. Massot, A. Vié, Reproducing segregation and particle dynamics in Large Eddy Simulation of particle-laden flows, in: International Conference on Liquid Atomization and Spray Systems (ICLASS), volume 1, 2021.

[2] R. Letournel, Reduced-order modeling and simulation of turbulent disperse two-phase flows: new theoretical and modeling approaches for reproducing intermittency, segregation and two-way coupling, Ph.D. thesis, Université Paris-Saclay, 2022. https://theses.hal.science/tel-03665513.

[3] R. Letournel, F. Laurent, M. Massot, A. Vié, Modulation of homogeneous and isotropic turbulence by sub-Kolmogorov particles: Impact of particle field heterogeneity, International Journal of Multiphase Flow 125 (2020) 103233. URL: https://doi.org/10.1016/j.ijmultiphaseflow.2020.103233. doi:10.1016/j.ijmultiphaseflow.2020.103233.

[4] J. G. Kuerten, Point-Particle DNS and LES of Particle-Laden Turbulent flow - a state-of-the-art review, Flow, Turbulence and Combustion 97 (2016) 689-713. URL: http://dx.doi.org/10.1007/s10494-016-9765-y. doi:10.1007/s10494-016-9765-y.

[5] C. Marchioli, Large-eddy simulation of turbulent dispersed flows: a review of modelling approaches, Acta Mechanica 228 (2017) 741-771. doi:10.1007/s00707-017-1803-x.
[6] M. W. Reeks, On the continuum equations for dispersed particles in nonuniform flows, Physics of Fluids A 4 (1992) 1290-1303. doi:10.1063/1.858247.

[7] J. Pozorski, J.-P. Minier, Probability density function modeling of dispersed two-phase turbulent flows, Physical Review E 59 (1999) 855. doi:10.1103/PhysRevE.59.855.

[8] J.-P. Minier, E. Peirano, The pdf approach to turbulent polydispersed two-phase flows, Physics reports 352 (2001) 1-214. doi:10.1016/S0370-1573(01)00011-4.

[9] A. Innocenti, C. Marchioli, S. Chibbaro, Lagrangian filtered density function for LES-based stochastic modelling of turbulent particle-laden flows, Physics of Fluids 28 (2016). doi:10.1063/1.4967800.

[10] A. Innocenti, R. O. Fox, M. V. Salvetti, S. Chibbaro, A lagrangian probability-density-function for collisional turbulent fluid-particle flows, Journal of Fluid Mechanics 862 (2019) 449-489.

[11] A. Innocenti, R. O. Fox, S. Chibbaro, A lagrangian probability-density-function model for turbulent particle-laden channel flow on the dense regime, Physics of Fluids 33 (2021) 053308.

[12] D. C. Haworth, S. B. Pope, A generalized Langevin model for turbulent flows, The Physics of fluids 29 (1986) 387-405. doi:10.1063/1.865723.

[13] S. B. Pope, Y. L. Chen, The velocity-dissipation probability density function model for turbulent flows, Physics of Fluids A 2 (1990) 1437-1449. doi:10.1063/1.857592.

[14] B. Shotorban, F. Mashayek, A stochastic model for particle motion in large-eddy simulation, Journal of Turbulence 7 (2006) N18. doi:10.1080/14685240600595685.

[15] P. Fede, O. Simonin, Numerical study of the subgrid fluid turbulence effects on the statistics of heavy colliding particles, Physics of Fluids 18 (2006) 45103. doi:10.1063/1.2189288.
[16] M. Bini, W. P. Jones, Large-eddy simulation of particle-laden turbulent flows, Journal of Fluid Mechanics 614 (2008) 207–252. doi:10.1017/S0022112008003443.

[17] J. Pozorski, S. V. Apte, Filtered particle tracking in isotropic turbulence and stochastic modeling of subgrid-scale dispersion, International Journal of Multiphase Flow 35 (2009) 118–128. URL: http://dx.doi.org/10.1016/j.ijmultiphaseflow.2008.10.005. doi:10.1016/j.ijmultiphaseflow.2008.10.005.

[18] R. Zamansky, I. Vinkovic, M. Gorokhovski, Acceleration in turbulent channel flow: Universalities in statistics, subgrid stochastic models and an application, Journal of Fluid Mechanics 721 (2013) 627–668. doi:10.1017/jfm.2013.48.

[19] J. P. Minier, S. Chibbaro, S. B. Pope, Guidelines for the formulation of Lagrangian stochastic models for particle simulations of single-phase and dispersed two-phase turbulent flows, Physics of Fluids 26 (2014) 113303. URL: http://dx.doi.org/10.1063/1.4901315. doi:10.1063/1.4901315. arXiv:1501.02498.

[20] M. Knorps, J. Pozorski, An Inhomogeneous Stochastic Model for Subgrid-Scale Particle Dispersion in LES, in: Direct and large-eddy simulation IX, Springer, 2015, pp. 671–678.

[21] M. W. Reeks, On a kinetic equation for the transport of particles in turbulent flows, Physics of Fluids A 3 (1991) 446–456. doi:10.1063/1.858101.

[22] J. P. Minier, On Lagrangian stochastic methods for turbulent polydisperse two-phase reactive flows, Progress in Energy and Combustion Science 50 (2015) 1–62. URL: http://dx.doi.org/10.1016/j.pecs.2015.02.003. doi:10.1016/j.pecs.2015.02.003.

[23] K. D. Squires, J. K. Eaton, Measurements of particle dispersion obtained from direct numerical simulations of isotropic turbulence, Journal of Fluid Mechanics 226 (1991) 1. URL: http://www.journals.cambridge.org/abstract{\_}S0022112091002276. doi:10.1017/S0022112091002276.
[24] C. Marchioli, M. V. Salvetti, A. Soldati, Appraisal of energy recovering sub-grid scale models for large-eddy simulation of turbulent dispersed flows, Acta Mechanica 201 (2008) 277-296. doi:10.1007/s00707-008-0066-y.

[25] B. Shotorban, F. Mashayek, Modeling subgrid-scale effects on particles by approximate deconvolution, Physics of Fluids 17 (2005) 81701. doi:10.1063/1.2001691.

[26] M. J. Cernick, S. W. Tullis, M. F. Lightstone, Particle subgrid scale modelling in large-eddy simulations of particle-laden turbulence, Journal of Turbulence 16 (2015) 101-135. URL: http://dx.doi.org/10.1080/14685248.2014.969888. doi:10.1080/14685248.2014.969888.

[27] A. Scotti, C. Meneveau, A fractal model for large eddy simulation of turbulent flow, Physica D: Nonlinear Phenomena 127 (1999) 198-232.

[28] E. Akinlabi, M. Waclawczyk, S. P. Malinowski, J. P. Mellado, Fractal reconstruction of sub-grid scales for large eddy simulation, Flow, Turbulence and Combustion 103 (2019) 293-322.

[29] R. H. Kraichnan, Diffusion by a random velocity field, Physics of Fluids 13 (1970) 22-31. doi:10.1063/1.1692799. arXiv:1011.1669v3.

[30] J. C. H. Fung, J. C. R. Hunt, N. A. Malik, R. J. Perkins, Kinematic simulation of homogeneous turbulence by unsteady random Fourier modes, Journal of Fluid Mechanics 236 (1992) 281-318. doi:10.1017/S0022112092001423.

[31] N. A. Malik, J. C. Vassilicos, A Lagrangian model for turbulent dispersion with turbulent-like flow structure: Comparison with direct numerical simulation for two-particle statistics, Physics of Fluids 11 (1999) 1572-1580. doi:10.1063/1.870019.

[32] J. C. H. Fung, J. C. Vassilicos, Two-particle dispersion in turbulent flows, Physical Review E 57 (1998) 1677. doi:10.1103/PhysRevE.57.1677.
[33] M. A. Khan, A. Pumir, J. C. Vassilicos, Kinematic simulation of turbulent dispersion of triangles, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 68 (2003) 13. doi:10.1103/PhysRevE.68.026313.

[34] D. J. Thomson, B. J. Devenish, Particle pair separation in kinematic simulations, Journal of Fluid Mechanics 526 (2005) 277-302. doi:10.1017/S0022112004002915.

[35] P. Flohr, J. C. Vassilicos, A scalar subgrid model with flow structure for large-eddy simulations of scalar variances, Journal of Fluid Mechanics 407 (2000) 315-349. doi:10.1007/978-94-011-5118-4_140.

[36] Z. Z., S. Wang, G. Jin, A structural subgrid-scale model for relative dispersion in large-eddy simulation of isotropic turbulent flows by couplin kinematic simulation with approximate deconvolution method, Physics of Fluids 30 (2018) 105110.

[37] B. Ray, L. R. Collins, A subgrid model for clustering of high-inertia particles in large-eddy simulations of turbulence, Journal of Turbulence 15 (2014) 366-385. URL: http://dx.doi.org/10.1080/14685248.2014.909600. doi:10.1080/14685248.2014.909600.

[38] M. Hausmann, F. Evrard, B. van Wachem, Large eddy simulation model for two-way coupled particle-laden turbulent flows, Physical Review Fluids 8 (2023) 084301.

[39] S. Murray, M. F. Lightstone, S. Tullis, Single-particle Lagrangian and structure statistics in kinematically simulated particle-laden turbulent flows, Physics of Fluids 28 (2016). URL: http://dx.doi.org/10.1063/1.4942815. doi:10.1063/1.4942815.

[40] P. G. Lemarié-Rieusset, Analyses multi-résolutions non orthogonales, commutation entre projecteurs et derivation et ondelettes vecteurs à divergence nulle, Revista Matemática Iberoamericana 8 (1992) 221-237. URL: http://hal.archives-ouvertes.fr/docs/00/38/16/49/PDF/these-deriaz.pdf.
[41] E. Deriaz, V. Perrier, Divergence-free and curl-free wavelets in 2D and 3D, application to turbulent flows, J. of Turbulence 7 (2006).

[42] E. Deriaz, Ondelettes pour la simulation des écoulements fluides incompressibles en turbulence. Ph. D. thesis, Institut National Polytechnique de Grenoble, 2006. URL: http://hal.archives-ouvertes.fr/docs/00/38/16/49/PDF/these-deriaz.pdf.

[43] M. Farge, Wavelet transforms and their applications to turbulence, Annual review of fluid mechanics 24 (1992) 395–458.

[44] C. Meneveau, Analysis of turbulence in the orthonormal wavelet representation, Journal of Fluid Mechanics 232 (1991) 469–520.

[45] R. Letournel, L. Goudenège, R. Zamansky, A. Vié, M. Massot, Reexamining the framework for intermittency in Lagrangian stochastic models for turbulent flows: A way to an original and versatile numerical approach, Physical Review E 104 (2021) 15104. URL: https://link.aps.org/doi/10.1103/PhysRevE.104.015104. doi:10.1103/PhysRevE.104.015104.

[46] R. Benzi, L. Biferale, A. Crisanti, G. Paladin, M. Vergassola, A. Vulpiani, A random process for the construction of multiaffine fields, Physica D: Nonlinear Phenomena 65 (1993) 352–358. doi:10.1016/0167-2789(93)90060-E.

[47] B. A. Zeldin, Representation and Synthesis of Random Fields: ARMA, Galerkin and Wavelet Procedures, Ph.D. thesis, Rice University, 1996.

[48] T. Kitagawa, T. Nomura, A wavelet-based method to generate artificial wind fluctuation data, Journal of Wind Engineering and Industrial Aerodynamics 91 (2003) 943–964. doi:10.1016/S0167-6105(03)00037-0.

[49] C. Yin, T. Wu, A. Kareem, Synthetic turbulence: A wavelet based simulation, Probabilistic Engineering Mechanics 45 (2016) 177–187. doi:10.1016/j.probengmech.2016.05.001.
[50] M. Hausmann, F. Evrard, B. van Wachem, Wavelet based modeling of subgrid-scales in LES of particle-laden turbulent flows, Physical Review Fluids 8 (2023) 104604. arXiv:2305.09521.

[51] K. Fröhlich, L. Schneiders, M. Meinke, W. Schröder, Validation of Lagrangian two-way coupled point-particle models in large-eddy simulations, Flow, Turbulence and Combustion 101 (2018) 317–341. doi:10.1007/s10494-018-9933-3.

[52] S. B. Pope, Turbulent Flows, Cambridge university press, 2000. doi:10.1088/0957-0233/12/11/705.

[53] R. H. A. Ijzermans, E. Meneguz, M. W. Reeks, Segregation of particles in incompressible random flows: singularities, intermittency and random uncorrelated motion, Journal of fluid mechanics 653 (2010) 99–136. doi:10.1017/S0022112010000170.

[54] K. Gustavsson, E. Meneguz, M. Reeks, B. Mehlig, Inertial-particle dynamics in turbulent flows: caustics, concentration fluctuations and random uncorrelated motion, New Journal of Physics 14 (2012) 115017. arXiv:1206.2018v2.

[55] S. Murray, M. F. Lightstone, S. Tullis, Target Lagrangian kinematic simulation for particle-laden flows, Physical Review E 94 (2016) 1–9. doi:10.1103/PhysRevE.94.033303.

[56] D. Mercier, Large eddy simulation of coupled dispersed phase flows: a statistically-consistent formalism. Ph. D. thesis, Université Paris-Saclay, 2020.

[57] N. A. Malik, Structural diffusion in 2D and 3D random flows, in: S. Gavrilakis, L. Machiels, P. A. Monkewitz (Eds.), Advances In Turbulence VI, Springer Netherlands, Dordrecht, 1996, pp. 619–620.

[58] D. R. Osborne, J. C. Vassilicos, J. D. Haigh, One-particle two-time diffusion in three-dimensional homogeneous isotropic turbulence, Physics of Fluids 17 (2005). doi:10.1063/1.1852578.
[59] I. Daubechies, The wavelet transform, time-frequency localization and signal analysis, IEEE Transactions on Information Theory 36 (1990) 961–1005. doi:10.1109/18.57199.

[60] I. Daubechies, Ten Lectures on Wavelets, Society for Industrial and Applied Mathematics, 1992. URL: https://epubs.siam.org/doi/abs/10.1137/1.9781611970104. doi:10.1137/1.9781611970104.arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9781611970104.

[61] Y. Meyer, Wavelets and Operators, Cambridge Studies in Advanced Mathematics, Cambridge University Press, 1995.

[62] J.-P. Kahane, P.-G. Lemarié-Rieusset, Fourier series and wavelets, volume 3, Gordon & Breach Science Pub, 1995.

[63] M. Farge, Wavelet transforms and their applications to turbulence, Annual review of fluid mechanics 24 (1992) 395–458. doi:10.5860/choice.39-6472.

[64] M. Farge, N. Kevlahan, V. Perrier, E. Goirand, Wavelets and turbulence, Proceedings of the IEEE 84 (1996) 639–669.

[65] A. Cohen, Wavelet methods in numerical analysis, in: Solution of Equation in $\mathbb{R}^n$ (Part 3), Techniques of Scientific Computing (Part 3), volume 7 of Handbook of Numerical Analysis, Elsevier, 2000, pp. 417–711. URL: https://www.sciencedirect.com/science/article/pii/S1570865900070046. doi:https://doi.org/10.1016/S1570-8659(00)07004-6.

[66] M. Duarte, M. Massot, S. Descombes, C. Tenaud, T. Dumont, V. Louvet, F. Lau rent, New resolution strategy for multiscale reaction waves using time operator splitting, space adaptive multiresolution, and dedicated high order implicit/explicit time integrators, SIAM Journal on Scientific Computing 34 (2012) A76–A104.

[67] M. Duarte, S. Descombes, C. Tenaud, S. Candel, M. Massot, Time-space adaptive numerical methods for the simulation of combustion fronts, Combustion and Flame 160 (2013) 1083–1101.
[68] N. K. R. Kevlahan, O. V. Vasilyev, An adaptive wavelet collocation method for fluid-structure interaction at high Reynolds numbers, SIAM Journal on Scientific Computing 26 (2005) 1894–1915.

[69] C. M. Albukrek, K. Urban, D. Rempfer, J. L. Lumley, Divergence-Free Wavelet Analysis of Turbulent Flows, Journal of Scientific Computing 17 (2002) 49–66.
doi:10.1023/a:1015184110888.

[70] K. Urban, Wavelet bases for H (div) and H (curl), in: Wavelets in Numerical Simulation, Springer, 2002, pp. 83–107.

[71] P.-G. Lemarié-Rieusset, Un théorème d'inexistence pour les ondelettes vecteurs à divergence nulle, Comptes rendus de l'Académie des sciences. Série 1, Mathématique 319 (1994) 811–813.

[72] J.-P. Kahane, Sur le chaos multiplicatif, Ann. Sci. Math. Québec 9 (1985) 105–150.

[73] A. N. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Cr Acad. Sci. URSS 30 (1941) 301–305.
doi:10.1098/rspa.1991.0075.

[74] R. Zamansky, F. Coletti, M. Massot, A. Mani, Turbulent thermal convection driven by heated inertial particles, Journal of Fluid Mechanics 809 (2016) 390–437.
doi:10.1017/jfm.2016.630.

[75] S. Chen, G. D. Doolen, R. H. Kraichnan, Z. S. She, On statistical correlations between velocity increments and locally averaged dissipation in homogeneous turbulence, Physics of Fluids A 5 (1992) 458–463.
doi:10.1063/1.858897.

[76] A. Kaufmann, O. Simonin, T. Poinsot, Direct numerical simulation of particle-laden homogeneous isotropic turbulent flows using a two-fluid model formulation, in: 5th Int. Conf. on Multiphase Flow, Yokohama, Japan, 443, 2004, p. 101.

[77] J. Eaton, J. Fessler, Preferential concentration of particles by turbulence, International Journal of Multiphase Flow 20 (1994) 169–209.
doi:10.1016/0301-9322(94)90072-8.
[78] M. Farge, K. Schneider, Wavelets: application to turbulence, University Warnick, lectures (2005).