This equation differs from the exact Navier-Stokes equation because the filtered product $\langle \mathbf{U}_f \mathbf{U}_f \rangle_{ls}$ is not the same as the product of the filtered velocities $\langle \mathbf{U}_f \rangle_{ls} \langle \mathbf{U}_f \rangle_{ls}$. In fact, the difference is the residual stress tensor $\tau_f$ (also called the subgrid-stress tensor). This tensor has to be modeled and a range of approaches exist in the literature (such as the famous Smagorinsky model which is an eddy-viscosity model [205]). LES methods are consistent with point-particle approximations but have to be adapted since the exact instantaneous velocity at the position of each particle is not known. In fact, the unresolved part of the velocity field has to be modeled for each particle, using an approach that has to be consistent with the subgrid-scale model used in the fluid simulation.

C. At a more macroscopic level, Reynolds-Averaged Navier-Stokes (RANS) simulations consist in writing mean-flow equations. These methods are based on the Reynolds decomposition of the fluid velocity into its mean and fluctuating parts:

$$ \mathbf{U}_f(t, \mathbf{x}) = \langle \mathbf{U}_f \rangle(t, \mathbf{x}) + \mathbf{u}'_f(t, \mathbf{x}). \quad (42) $$

Applying this true probabilistic averaging operator to the Navier-Stokes equations gives the Reynolds (or mean Navier-Stokes) equation:

$$ \frac{\partial \langle \mathbf{U}_f \rangle}{\partial t} + \langle \mathbf{U}_f \rangle \cdot \nabla \langle \mathbf{U}_f \rangle = -\nabla \cdot (\langle \mathbf{u}'_f \otimes \mathbf{u}'_f \rangle) - \frac{1}{\rho_f} \nabla \langle P \rangle + \nu_f \nabla^2 \langle \mathbf{U}_f \rangle + \langle \mathbf{F} \rangle. \quad (43) $$

where the first term on the right hand side is the Reynolds stress tensor which has to be modeled. In practice, two classes of models are traditionally used: (i) local models, where the Reynolds stress tensor is approximated using an eddy-viscosity concept in which this turbulent viscosity is based on known local values of the fluid (such as the turbulent kinetic energy $k$ and the turbulent dissipation rate $\epsilon$ in the famous $k-\epsilon$ models); (ii) non-local models, which solve additional equations for the transport of each component of the Reynolds stress tensor (like Reynolds-stress models, RSM).

RANS approaches are consistent with point-particle approximations but require an additional model to simulate the instantaneous fluid velocity at each particle position. Indeed, RANS simulations provide only one-point statistical information on the flow fields. This means that the instantaneous fluid velocity at any point in space and time is not directly available but can be reconstructed - in a statistical sense - using information on the average velocity and its fluctuating part (more details in [8, 9, 217]).

### 4.2.2. Surface roughness

**Physical origin:** As displayed in Fig. 28, the surfaces of particles and substrates are rarely smooth and usually present irregularities. In fact, the surfaces of raw materials (e.g., on a stainless steel substrate or on a micrometer-size plastic particle) appear as quite chaotic with roughness features covering a range of sizes (from a few nanometers up to a few micrometers here), whereas more regular roughness, or ‘smooth surfaces’, can only be obtained by using special surface engineering techniques (e.g., micro-patterned surfaces obtained by photolithography) or when handling biological particles (like pollen grains). This explains that surface topologies are regarded as random and are, at best, characterized by a few averaged dimensions such as [218]: the mean roughness (i.e., the arithmetic average height above a reference plane), the rms roughness (i.e., root-mean-square roughness) and the peak-to-peak distance. However, more detailed statistical information is needed to properly account for the effect of surface roughness and this point is taken up in Section 7.1.2.

**Role in resuspension:** Surface roughness plays a central role in resuspension through the following phenomena:

* First, as mentioned previously, surface irregularities lead to the existence of friction forces.

* Second, roughness features with sizes smaller or comparable to the particle size can have a profound impact on adhesion forces (especially VDW forces). Actually, small asperities on the surface lead often to a significant reduction of the actual adhesion force compared to the one that would exist if we were dealing with perfectly smooth surfaces. This effect of surface roughness on adhesion forces is now well documented (see previous reviews, e.g., [222, 11, 10] and specific papers on the topic [223, 219, 187],
<figure><img src="image_36.png" /><figcaption>(a) AFM image showing the surface topography of stainless steel. Reprinted with permission from [219]. Copyright 2010, American Chemical Society.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(c) SEM image of a polybutylene terephthalate particle. Reprinted with permission from [219]. Copyright 2010, American Chemical Society.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(b) SEM image showing arrays of pillars made using photolithographic techniques (hexagonal packing of pillars with a radius of 2.5 µm and a height of 20 µm). Reprinted with permission from [220]. Copyright 2007, American Chemical Society.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(d) SEM image of a ragweed pollen grain glued to the end of a tipless AFM. Reprinted with permission from [221]. Copyright 2009, American Chemical Society.</figcaption></figure>

Figure 28: Images showing the roughness of two substrates (top) and two particles (bottom). Surface roughness can either appear as very chaotic (left, with raw substrate or particles) or highly regular (right, with engineered surfaces or pollen grain).

as well as reference therein) and can be summarized as follows (see also [224]): protruding roughness features keep particles further away from the rest of the surface, often at separation distances over which van der Waals forces are significantly reduced. Hence, the particle-substrate interaction is mostly governed by roughness features that are in close contact, giving rise to adhesion forces much weaker than the ones obtained when assuming perfectly smooth surfaces (which allow to have large surface areas in close proximity). The overall effect is sketched in Fig. 29, which illustrates that surface roughness is responsible for a distribution of adhesion forces between surfaces. Given that surface roughness is described as random, this explains that adhesion forces are also naturally treated as random variables, characterized by their probabilistic distribution or PDF. In practice, various types of distributions have been measured for the adhesion forces between rough surfaces, such as Gaussian or Weibull distributions (e.g., [187]) or even multimodal distributions (as in [225, 224, 125]).

* Third, roughness features with sizes comparable to or larger than the particle size can modify the hydrodynamic forces acting on particles. More precisely, two separate effects can take place. First, large-scale roughness features can shelter particles located in their wake from the flow, effectively reducing resuspension. Second, when roughness features display sizes comparable to the smallest fluid scale (typically the Kolmogorov scale in a turbulent flow), the fluid flow around such geometrical features is modified (e.g., the wind flow around sand dunes [1]). Depending on the relative size between surface roughness and particle diameter, this can either facilitate resuspension or hinder it [227].

**Modeling surface roughness:** When developing resuspension models, it follows from the previous arguments that the effects of surface roughness on adhesion forces have to be carefully included. This means that information at the nano- and micro- scales have to be accounted for, even if these scales are
<figure><img src="image_40.png" /><figcaption>Figure 29: Illustration of the effect of surface roughness on adhesion forces (which can be significantly reduced when roughness features are present between the bulk of each materials). Experimental data taken from [226].</figcaption></figure>

much smaller than the typical ones of a fluid flow (which range from a few millimeters to a few meters). As a result, surface roughness is often simplified in resuspension models and various descriptions can be found. In the following, we briefly review the existing approaches to represent rough surfaces and how these descriptions are involved in the evaluation of adhesion forces (see Fig. 30 and [10, 187, 219] for more details):

* A first nanoscopic approach consists in using the exact information on the surface topography to compute the resulting forces explicitly. This is only possible if the surface topography has been fully characterized using detailed measurement techniques (such as SEM of AFM methods). The exact topography can then be used to compute the adhesion force $F_{adh}$ at any contact point between a spherical particle and this rough substrate [230]. This is usually achieved by discretizing the volume of each body (e.g., using Finite Element Method) and then performing a direct numerical integration of van der Waals forces similar to the Surface Element Integration technique given in Eq. (25) [179].

* Other descriptions of surface roughness rely on more limited information on the surface topology. One way to represent a rough surface, still with a relatively high level of information, is to resort to a Fourier transform of the surface [231]. This allows to generate numerical rough substrates that closely approach the real ones. The adhesion force between such a rough substrate and a particle can be computed using similar numerical integration techniques (such as SEI methods).

* A coarser description of rough surfaces can be made by resorting to fractal surfaces (especially with Weierstrass-Mandelbrot's functions [228, 229]). Once again, this allows to randomly generate various numerically rough surfaces [229] and, using numerical integration techniques, to estimate the distribution of adhesion forces.

* In practice, it is very frequent to have access only to average information on surface roughness (such as the rms roughness $r_{a,rms}$, the peak-to-peak distance, the extremum sizes of roughness features). In that case, the issue is how to estimate the adhesion force. Depending on the level of information required, two types of model have been suggested:

- When the objective is to approximate the complete distribution of adhesion forces, various models have proposed to represent rough surfaces by a smooth plate covered by hemispherical asperities randomly distributed on the surface [223]. The size of these asperities (resp. their spacing) is determined by the available data on the rms roughness (resp. the peak-to-peak distance). Then, using the Hamaker approach, the total adhesion force is simply given by summing the interaction between the particle and all the elements composing the surface (i.e., all the hemispherical asperities plus the smooth substrate beneath):

$$F_{adh}(\text{sphere} - \text{roughS}) = F_{adh}(\text{sphere} - \text{smoothS}) + \sum_{\text{all asperities}} F_{adh}(\text{sphere} - \text{asperity}). \quad (44)$$
<figure><img src="image_41.png" /><figcaption>Figure 30: Sketch of existing descriptions of surface roughness, showing both exact topologies and reduced descriptions (sorted according to the level of information contained). Top 2 figures reprinted with permission from [228]. Copyright 2004, Elsevier. Middle figure reprinted with permission from [229]. Copyright 2007, IOP Publishing.</figcaption></figure>

By repeating the procedure to generate different realizations of the rough surface, the distribution of adhesion forces can be estimated through Monte-Carlo methods [83, 218].

- When the aim is to predict only the average value of the adhesion force, one can rely directly on analytic formulas. These analytic formulas are obtained by summing the interaction between a spherical particle and all the elements forming the rough surface (i.e., the hemispherical asperities and the smooth surface underneath). Yet, instead of placing asperities randomly, the analytic formulas are obtained either by placing one asperity at a fixed position (e.g., right beneath the particle as in [232]) or by regularly placing asperities on the surface (as in [233, 234]). For instance, in the case of multiple asperities present on the surface, the average value of the adhesion force can be estimated with:

$$F_{\text{adh}} = \frac{A_{\text{Ham}} r_p}{6z_0^2} \left( \frac{1}{1 + \frac{58.144 r_p r_{a,rms}}{\lambda_a^2}} + \frac{1}{\left(1 + \frac{1.817 r_{a,rms}}{z_0}\right)^2} \right) \quad (45)$$
where $A_{\text{Ham}}$ is the Hamaker constant (which reflects the strength of the adhesion between two objects), $z_0$ the contact distance, $r_p$ the particle radius, $r_{\text{a,rms}}$ the rms roughness, and $\lambda_a$ the peak-to-peak distance between consecutive asperities.

More recently, a similar model has been developed by describing the rough surface as a sinusoidal wave surface [227]. The amplitude of the sinusoidal function is given by the average roughness size while the frequency is proportional to the peak-to-peak distance. The adhesion force is then simply obtained using contact mechanics theories (like the JKR model in [227]) and considering that the particle is resting on a stable position on the sinusoidal wave surface (i.e., two points of contact in 2D cases or 3 in 3D cases). Similarly, other formulations based on contact mechanics theories have been proposed to describe single or multiple asperity contacts (more details in [187] and references therein).

A few caveats are in order. First, these adhesion models only approximate the real force. For example, in the case of non-spherical particles having irregular shapes and their own roughness (which are more common in realistic situations, as for sand beads), the particle and the surface can come in contact at multiple points [235], making the resulting adhesive force dependent on the orientation of the particle [236]. Second, contact mechanics theories assume that deformation is much smaller than any characteristic dimension of the body. Yet, this assumption is not always true for small colloidal particles and small roughness elements [237]. Third, the simplified representations based on hemispherical asperities have been derived from information extracted from 2D profiles of surfaces. Yet, a recent study [237] has shown that it is not always straightforward to use 3D information, as protuberances on the surface do not necessarily resemble parabolas. Fourth, more general theories that account for plastic deformations are still not available [238].

### 4.3. A wide range of approaches

Although various attempts have been made to suggest classifications of these models (see recent reviews [11, 10, 122, 239]), no consensus has emerged since different points of views can be adopted to categorize modeling approaches. In the following, we draw on the classification proposed in [10] to propose a hierarchy of existing models. These models are sorted according to the level of information contained, starting with the most detailed models up to the more macroscopic ones (interested readers are directed to this previous review for a comprehensive discussion of different model classes). Note however that this classification is extended to include other models that were not accounted for in the previous review.

#### 4.3.1. *Dynamic models based on spatio-temporal integration*

**General description** The leading idea of exact dynamic models is to track a large number of individual particles within a domain by solving explicitly the particle equations of motion.

**Model formulation** These models are designed in a step-by-step manner summarized as follows:

* The first step consists in defining the particle state vector $Z_p$ which gathers the variables attached to each particle and which are selected to describe each particle as a mechanical object. For instance, one can choose to retain the position $X_p$ and the velocity $U_p$ of each individual particle. This corresponds to having $Z_p = (X_p, U_p)$, which is the classical state vector in Statistical Physics when dealing with kinetic approaches. While this choice appears relevant to capture the translational motion of rigid spherical particles (e.g., sliding or lifting motion), it does not allow to track rolling motion of solid spheres on a surface since the rotational velocity $\Omega_p$ is missing. For that reason, most resuspension models are based on the state vector $Z_p = (X_p, U_p, \Omega_p)$. Additional variables can be introduced in the state vector, such as the particle size $r_p$ when dealing with aggregates or the particle orientation $\mathbf{p}_p$ when dealing with non-spherical particles (see Section 7).

* Once the particle state vector $Z_p$ is chosen, the second step consists in writing an equation for the evolution of these quantities in time. For the sake of clarity, we consider only here the simplest case of rigid spherical particles where the state vector is $Z_p = (X_p, U_p, \Omega_p)$ (the particle diameter or radius
is a constant for each particle and is left out of the state vector). In that case, the equations of translational and rotational motion are given by Newton's second law and write:

$$ \frac{d\mathbf{X}_p}{dt} = \mathbf{U}_p, \qquad (46a) $$

$$ m_p \frac{d\mathbf{U}_p}{dt} = \mathbf{F}_{f \to p} + \mathbf{F}_{s \to p} + \mathbf{F}_{p \to p} + \mathbf{F}_{\text{ext}}, \qquad (46b) $$

$$ I_p \frac{d\Omega_p}{dt} = \mathbf{M}_{f \to p} + \mathbf{M}_{s \to p} + \mathbf{M}_{p \to p} + \mathbf{M}_{\text{ext}}, \qquad (46c) $$

with $m_p$ the particle mass and $I_p$ its moment of inertia. The terms on the right hand side of Eqs. (46) correspond to the forces $\mathbf{F}$ and torques $\mathbf{M}$ acting on the particles. This includes terms related to fluid-particle interactions (.)$_{f \to p}$, surface-particle interactions (.)$_{s \to p}$, particle-particle interactions (.)$_{p \to p}$, and external interactions (.)$_{\text{ext}}$ (see also Section 4.1 for a physical interpretation of the forces at play).

In the present hierarchy of models, dynamics models embodied by Eqs. (46) represent the most detailed level of description or, in other words, fine-grained formulations. In order to obtain a closed set of equations, they require however that all the forces and torques appearing on the rhs of Eqs. (46) be explicitly known. Given the expressions of these forces and torques given in Section 4.1, it is useful to emphasize that they rely on detailed information on the instantaneous fluid velocity around each particle for fluid-particle interactions while accurate expressions for the adhesive forces between rough substrate can only be obtained if information on the nanoscale surface topology is available.

Eqs. (46) are expressed for the variables entering the one-particle state vector $\mathbf{Z}_p$, but, since we are tracking N particles simultaneously, the actual state vector is $\mathbf{Z} = (\mathbf{Z}_p^{(1)}, \mathbf{Z}_p^{(2)}, \dots, \mathbf{Z}_p^{(N)})$, where $\mathbf{Z}_p^{(i)}$ stands for the state vector of the particle labeled (i). The governing equations can then be written in a compact form for the N-particle state vector $\mathbf{Z}$ as a general ordinary differential equation (ODE) in a $N \times N_p$ dimensional space ($N_p$ being the number of variables in the one-particle state vector $\mathbf{Z}_p$)

$$ \frac{d\mathbf{Z}}{dt} = \mathcal{A}(t, \mathbf{Z}(t)), \qquad (47) $$

where the function $\mathcal{A}$ gathers all the terms on the rhs of Eqs. (46). In a statistical formulation, we are therefore dealing with a N-particle PDF approach, in which the PDF $p(t; \mathbf{z})$ (where $\mathbf{z}$ represents the possible values taken by the state vector $\mathbf{Z}$) is the solution of a Liouville equation in its corresponding sample space:

$$ \frac{\partial p(t; \mathbf{z})}{\partial t} = - \frac{\partial [\mathcal{A}_k(t; \mathbf{z}) p(t; \mathbf{z})]}{\partial z_k}, \quad k = 1, \dots, N \times N_p. \qquad (48) $$

**Typical output and characteristics** As it transpires from the previous discussion on the choice of a state vector, exact dynamic models can provide very detailed information at the particle scale (e.g., position, translational and rotational velocity of each individual particle). They allow naturally to predict the history of individual particles (see Section 3.1.1) and, by calculating ensemble averages or using other adequate statistical analysis, to obtain various statistical quantities or factors of interest in each situation (see the list in Section 3.1.1). In addition, when specific models are used to compute the adhesion force between rough objects, these approaches yield also microscopic predictions, such as the distribution of adhesion forces. For example, the 10 % quantile threshold velocity for incipient motion can be estimated by running numerical simulations with various fluid velocities and detecting when more than 10 % of the particles have started moving.

More generally, such dynamic models are naturally consistent with Lagrangian approaches for particle transport by fluid flows, which consist in solving the same equations for the dynamics of particles (but not necessarily with particle-surface interactions which act only close to wall boundaries). Another motivation is that they are very flexible. For example, it is relatively straightforward to extend Eqs. (46) to treat the case of non-spherical particles: this requires to add an extra equation for the particle orientation $\mathbf{p}_p$ and to
modify the expressions for the hydrodynamic forces in order to account for the changes in the drag and lift forces due to the particle non-sphericity and its orientation (see developments in Section 7.1.1).

The main drawback of dynamic approaches is the high computational cost, as well as the near-impossibility to have access to the exact topography of wall surfaces in most practical cases. In fact, simulations can hardly be performed in realistic environmental applications (such as resuspension from sand dunes), since they would require to follow the trajectories of an extremely large number of particles (typically $10^{12}$ or even larger), while still suffering from uncertainties due to choices made in the description of surface characteristic geometrical features. Besides, the time step used has to be small enough to capture fluctuations in adhesive forces and the models for the adhesion force between rough objects often require extensive computations.

**Current use** Given that these models are naturally compatible with N-particle tracking within a flow, most of them are coupled to the Distinct Element Method (DEM), which solves equations of motion similar to Eqs. (46) anywhere in the physical domain (and not only on the surface). Due to the high level of information contained in these formulations, they are usually coupled to a Direct Numerical Simulations (DNS) of the turbulent flow carrying these particles. In practice, the two methods described in Section 4.2.1, namely the PR-DNS and PP-DNS, have been used:

* Using the PR-DNS formulation, an Immersed Boundary Method has been used in a recent paper [215] to discretize the surface of each particle and apply a zero-flux condition for the fluid velocity on each particle surface. As displayed in Fig 31, these simulations allow to obtain precise information on the flow around 13500 particles forming a single monolayer. Numerical results showed that a particle in a dense monolayer bed is resuspended by the combination of two processes: first, it can be temporarily dislocated from the surface when collided by another particle transported by the flow; second, as it interacts with a sweep event, this particle can be entirely dislodged from its initial position.

<figure><img src="image_42.png" /><figcaption>Figure 31: Analysis of the fluid velocity around a dense monolayer deposit obtained using finite-size DNS revealing the role of coherent structures in resuspension: (a) Quadrant analysis showing all events reported for the fluctuating velocities (streamwise $u'^c$ and wall-normal $v'^c$ components); (b) iso-surfaces of the vorticity $\omega_x$ for the side view (top) and the top view (bottom); blue: $\omega_x = -0.4$, red: $\omega_x = 0.4$. Reprinted with permission from [215]. Copyright 2016, Taylor & Francis.</figcaption></figure>

Similar results were obtained using a Lattice-Boltzmann approach to compute the fluid flow around finite-size particles [240]. Considering that resuspension is triggered by a simple rupture of balance between the forces acting on each particle, the authors highlighted the role played by surface roughness in the modifications of the near-wall velocity and its potential consequences on particle resuspension.

* In the frame of the PP-DNS formulation, the hydrodynamic drag and lift forces are evaluated using one of the existing formulas presented in Section 4.1. At this stage, it is worth mentioning that the
precision of the value for the velocity of the fluid seen $\mathbf{U}_s(t)$ (with $\mathbf{U}_s(t) = \mathbf{U}_f(t, \mathbf{X}_p(t))$) depends on how it is evaluated numerically from the DNS simulations. In fact, particles are distributed within the computational domain and are not necessarily located at one of the grid points used in the DNS simulation. Hence, an exact value can be estimated by resorting to the Fourier transform (the same as the one used to solve the Navier-Stokes equation in the DNS simulation). However, only an approximate value is obtained when resorting to an interpolation method: in that case, the accuracy is driven by the precision of the interpolation algorithm (e.g., tri-linear interpolation or tri-cubic interpolation [241]).

Recently, similar N-body simulations have been carried out to investigate the role of impact entrainment in complex multilayer deposits. For that purpose, 15000 particles were initially deposited to form a deposit with more than 10 layers (see Fig 32). Particle motion is then computed neglecting buoyancy, lift, and shielding effects while inter-particle collisions are taken into account using a simple model based on the friction coefficient and the restitution coefficient (i.e., the amount of energy of the impacting particle remaining after the collision). This has allowed to confirm that repletion and saltation motion mainly differ due to differences in the kinetic energy upon collision: saltating particles display much more pronounced hoping motion, with higher energies upon collision with the surface, that can result in much more dislocations of deposited particles. Statistical analysis tools can then be used to extract more macroscopic information such as the threshold for particle incipient motion or the Shields parameter (see for instance [242]).

<figure><img src="image_43.png" /><figcaption>Figure 32: Visualizations of impact entrainment events for turbulent bedload (left) and saltation transport (right) from direct sediment transport simulations near the threshold conditions. The purple and green colors indicate the trajectory of particles before and after impact, respectively. Particles move from the left to the right. Reprinted with permission from [243]. Copyright 2017, American Physical Society.</figcaption></figure>

It is essential to remember that these models require very detailed information on all the forces acting on particles. In terms of level of information, it is thus more consistent to couple these approaches with DNS simulations of the fluid flow, but statistical information can also be derived from them provided that data is available on velocity distributions (this point is taken up in Section 6.4.1).

### 4.3.2. *Dynamic PDF approaches*

**General description** Broadly speaking, dynamic PDF models consist in coarse-grained descriptions of exact dynamic models in which particle variables are modeled as stochastic processes. In line with the general ideas of Synergetics, they correspond to reduced statistical descriptions where the fast degrees of freedom are eliminated to leave white-noise type of stochastic processes as their trace in the evolution equations of the remaining slow degrees of freedom [8, 244]. This amounts to adopting a probabilistic point of view [245, 8, 9], where the objective is to approximate relevant particle-related statistics rather than compute the exact trajectories of individual particles. In a weak formulation, this is equivalent to approximating the probability density function (PDF) of the variables selected to describe particle dynamics, hence the name of the approach. These methods retain a Lagrangian formulation but we are now dealing with stochastic particles which represent samples of this underlying PDF and are to be regarded as random versions of real
ones. In principle, it is possible to track a large number of particles, or pairs of particles, or triplets of particles, etc., which corresponds to handling a one-particle PDF, or a two-particle PDF, or a three-particle PDF, etc. (respectively). Yet, for practical applications in general non-homogeneous and wall-bounded flows, only one-particle PDF are well-developed at the moment. This means that they produce only one-point statistical outcomes (note that this point will resurface in Section 6.4.1).

**Model formulation** Like the exact dynamical model described above, the formulation of PDF approaches is achieved through a two-step construction. However, both steps are less straightforward and involve specific issues that need to be addressed to clarify what these PDF models can and cannot capture.

* The first step concerns the selection of the variables entering the particle state vector $Z_p$ for particles transported by a flow described by a turbulent model, which means that only limited statistical information, such as the one-point first- and second-order velocity moments, is available. To bring out the issues at stake, it is sufficient to consider the simple, but reference, case of small discrete particles heavier than the carrier fluid and described as point-wise particles. Then, the particle equations of motion have the form (retaining only drag and gravity forces)

$$ \frac{d\mathbf{X}_p}{dt} = \mathbf{U}_p, \qquad (49a) $$

$$ \frac{d\mathbf{U}_p}{dt} = \frac{\mathbf{U}_s - \mathbf{U}_p}{\tau_p} + \mathbf{g}, \qquad (49b) $$

where $\tau_p$ is the particle relaxation timescale. The drag force on the rhs of Eq. (49b) is written as a return-to-equilibrium term towards the fluid velocity seen $\mathbf{U}_s$, which is the instantaneous fluid velocity ‘sampled’ by discrete particles as they move across the fluid flow (i.e., $\mathbf{U}_s(t) = \mathbf{U}_f(t, \mathbf{X}_p(t))$ where $\mathbf{U}_f(t, \mathbf{x})$ is the instantaneous fluid velocity field). If we retain the classical kinetic particle state vector, $\mathbf{Z}_p = (\mathbf{X}_p, \mathbf{U}_p)$, then $\mathbf{U}_s$ appears as an external driving force. In the context of DNS and the exact dynamical models considered before, $\mathbf{U}_s$ is a known deterministic signal, leading to the well-posed Liouville formulation. In the frame of modeled turbulent flows (like LES or RANS), $\mathbf{U}_s$ appears, however, as a stochastic process with non-zero correlation timescales (due to the non-zero space and time correlations of turbulent fluctuations). This deviation from ‘heat-bath properties’ leads to severe difficulties and results even in ill-based probabilistic formulations (this key issue is addressed in detail in [9, 244]). In line with the notion of slow/fast variables mentioned above, it is much better to include the fluid velocity seen in the particle state vector, so that $\mathbf{Z}_p = (\mathbf{X}_p, \mathbf{U}_p, \mathbf{U}_s)$. This means that the external driving force is now the fluid acceleration which, according to Kolmogorov turbulence theory, can be regarded as a fast process and has thus better chances to be safely replaced by a stochastic diffusion process, thereby ensuring well-posed probabilistic formulations.

* With the selection of the particle state vector, the second step consists in devising a stochastic model for the velocity of the fluid seen. This has been discussed in several reviews [8, 9, 246] and present state-of-the-art formulations consist in modeling $\mathbf{U}_s$ as the solution of a generalized Langevin equation. Basically, these models are extensions of the ones developed to simulate the dynamics of fluid-like particles for turbulent reactive flows [205, 247], providing a unified framework [217] which allows also new methodologies to be proposed [248]. The specific details of such Langevin equations can be found in previous works and, for our present concern, it is sufficient to write the evolution equations for the retained particle state vector as

$$ d\mathbf{X}_p = \mathbf{U}_p dt, \qquad (50a) $$

$$ d\mathbf{U}_p = \mathbf{D}_p(t, \mathbf{Z}_p) dt, \qquad (50b) $$

$$ d\mathbf{U}_s = \mathbf{D}_s(t, \mathbf{Z}_p, \mathcal{F}[\langle\mathbf{Z}_p\rangle], \langle\Phi\rangle) dt + \mathbf{B}_s(t, \mathbf{Z}_p, \mathcal{F}[\langle\mathbf{Z}_p\rangle], \langle\Phi\rangle) d\mathbf{W}. \qquad (50c) $$

In these equations, $\mathbf{D}_p$ typically represents the drag and gravity forces, as in Eq. (49b), but other forces can also be included. The vector $\mathbf{D}_s$ and the matrix $\mathbf{B}_s$ in Eq. (50c) are the drift and diffusion
coefficients of the stochastic diffusion process retained for the velocity of the fluid seen, while **W** is a vector of independent Wiener processes. In the drift and diffusion coefficients in Eq. (50c), a general notation has been used to indicate that these coefficients can depend on the value of the state-vector **Z**<sub>p</sub> but also on functionals of the mean fields which are calculated from the simulation of that state-vector, written as $\mathcal{F}[\langle\mathbf{Z}_p\rangle]$, as well as on some external fields represented by $\langle\Phi\rangle$. A typical example of $\mathcal{F}[\langle\mathbf{Z}_p\rangle]$ is the particle mean-velocity field while the fluid mean-pressure is another example of what $\langle\Phi\rangle$ can stand for. From Eqs. (50), it follows that **Z**<sub>p</sub> is a stochastic diffusion process whose trajectories are nowhere differentiable but still continuous. This is in line with the idea of modeling discrete particle agitated by the actions of random but continuous fluid motions. Note that for variables undergoing abrupt changes and having therefore discontinuous trajectories in sample space, such as particle velocities upon collisions, Poisson stochastic processes are more likely model candidates but this aspect is outside the scope of this presentation (interested readers are referred to [249, 250] among other references). Using then a compact notation, the general evolution equations in dynamical PDF models can therefore be formulated as:

$$d\mathbf{Z}_p = \mathcal{A}(t, \mathbf{Z}_p, \mathcal{F}[\langle\mathbf{Z}_p\rangle], \langle\Phi\rangle) dt + \mathcal{B}(t, \mathbf{Z}_p, \mathcal{F}[\langle\mathbf{Z}_p\rangle], \langle\Phi\rangle) d\mathbf{W}_t, \quad (51)$$

where $\mathcal{A}$ is the drift vector and $\mathcal{B}$ the diffusion matrix modeling the time-rate-of-change of $\mathbf{Z}_p$.

Contrary to the exact dynamical model where we are considering a N-particle PDF evolving through the Liouville equation, cf. Eq. (48), we are handling here a one-particle PDF and, in that sense, the actual particle state vector **Z** is identical to the particle one, i.e. $\mathbf{Z} = \mathbf{Z}_p$. In the present PDF formalism, the true quantity being simulated is the (one-particle) Lagrangian mass density function (MDF) which, for the chosen state vector $\mathbf{Z}_p = (\mathbf{X}_p, \mathbf{U}_p, \mathbf{U}_s)$, is defined as:

$$F_p^L(t; \mathbf{Y}_p, \mathbf{V}_p, \mathbf{V}_s) = \sum_{i=1}^{N} m_p^{(i)} \delta(\mathbf{Y}_p - \mathbf{X}_p) \delta(\mathbf{V}_p - \mathbf{U}_p) \delta(\mathbf{V}_s - \mathbf{U}_s), \quad (52)$$

where $\mathbf{Y}_p, \mathbf{V}_p$ and $\mathbf{V}_s$ stand for the sample space variables corresponding to $\mathbf{X}_p, \mathbf{U}_p$ and $\mathbf{U}_s$, respectively. The Eulerian MDF is simply defined as the Lagrangian one at a given location, which writes

$$F_p^E(t, \mathbf{x}; \mathbf{V}_p, \mathbf{V}_s) = F_p^L(t; \mathbf{Y}_p = \mathbf{x}, \mathbf{V}_p, \mathbf{V}_s) = \int F_p^L(t; \mathbf{Y}_p, \mathbf{V}_p, \mathbf{V}_s) \delta(\mathbf{Y}_p - \mathbf{x}) d\mathbf{Y}_p. \quad (53)$$

It is important to note that these are one-point MDFs and, therefore, only one-point mean fields can be extracted. Using the same compact notation as in Eq. (51) but leaving out the functional dependencies for the sake of keeping simple notations, the evolution equation in sample space followed by both $F_p^L$ and $F_p^E$ is now a Fokker-Planck equation

$$\frac{\partial F_p^E}{\partial t} = -\frac{\partial [\mathcal{A}_k F_p^E]}{\partial z_k} + \frac{\partial^2 [(\mathcal{B}\mathcal{B}^T)_{kl} F_p^E]}{\partial z_k \partial z_l}, \quad k, l = 1, \dots, N_p, \quad (54)$$

where $\mathcal{B}^T$ is the transpose of the matrix $\mathcal{B}$.

To exemplify how this general formalism is put in practice and to illustrate PDF dynamic models, it is instructive to consider a recent dynamic PDF approach for particle resuspension whose main aspects are briefly recalled here (more details can be found in [251, 83, 218]). This model has been developed to describe the motion of colloidal particles on a rough substrate in the case of a dispersed monolayer. These particles are small enough to be embedded in the viscous sublayer and, consequently, only rolling motion is considered. This rolling motion is assumed to be triggered by hydrodynamic drag forces only, while adhesive forces are acting to prevent it. The small amount of particles on the surface allows to treat each particle independently from one another (i.e., without accounting for particle-particle interactions). The model is based on a three-step scenario for particle resuspension (see also Fig. 33): first, deposited particles are set in motion when the balance of torques is ruptured; second, particles are migrating on the surface (possibly
<figure><img src="image_44.png" /><figcaption>Figure 33: Sketch showing the three-stage scenario used in a recent dynamic PDF approach, where resuspension is assumed to be driven by an incipient motion (when the balance of torques is ruptured) followed by particle migration on the surface and rocking events (where detachment occurs if the kinetic energy is higher than the adhesion energy). Reprinted with permission from [251]. Copyright 2012, American Chemical Society.</figcaption></figure>

accelerating/decelerating); third, particles can detach from the surface during rocking events (where particles impact large-scale roughness features that act as a ramp for detachment).

Assuming a linear relation between the particle angular velocity $\Omega_p$ and the streamwise translational velocity $U_{p,||} \simeq d_p\Omega/2$, there is no need to extend to explicitly include $\Omega_p$ in the state vector, which is kept equal to $Z_p = (X_p, U_p, U_s)$. The translational motion of particles is then given by:

$$I_p \frac{dU_{p,||}}{dt} \simeq R_p M_{drag,||} - \frac{d_p}{2} M_{adh}. \qquad (55)$$

Solving Eq. (55) requires to know the torques due to hydrodynamic drag and to adhesion forces but, contrary to the exact dynamic models, only statistical information are needed in dynamic PDF formulations. In this model, $M_{drag}$ is evaluated by resorting to a stochastic Langevin model to simulate the instantaneous fluid velocity encountered by particles near the boundary [252, 251] while the adhesive forces are approximated with Hamaker approach coupled to a simple model for substrate roughness (represented by a smooth surface covered by hemispherical asperities) [83, 218]. Hence, only statistical information on the fluid velocity (average velocity) and on the surface roughness (mean roughness size and surface coverage) are necessary to obtain a closed formulation and run numerical simulations. The third step in the scenario corresponds to the detachment criteria: when the particle kinetic energy is higher than the adhesive one during a rocking event, a particle is considered as being detached from the surface and, hence, resuspended. This third step requires, therefore, to determine the probability of occurrence of these rocking events. This probability is simply calculated using the surface 'seen' by the particle as it migrates on the surface together with experimental information on the percentage of large-scale roughness features covering the substrate.

**Typical output and characteristics** Although the exact velocity of each particle cannot be evaluated in dynamic PDF methods, the fact that particle velocity is explicitly retained in the particle state vector means that statistical information on particle velocities is readily available. As a result, these formulations provide access to a wide range of factors, like velocity statistics, any quantile threshold velocity for incipient motion, any quantile critical velocity, and, of course, resuspension rates, or remaining/left fractions.

Similarly to the dynamic models, a further interest in resorting to dynamic PDF approaches is that they are flexible enough to include additional information. For instance, following recent observations of collision propagation events in dense monolayer deposits, dynamic PDF approaches are now extended to account for the role of inter-particle collision in triggering the motion of individual particles [110]. For that purpose, the state vector is left unchanged $Z_p = (X_p, U_p, U_s)$ but an additional treatment is introduced to evaluate the probability that a migrating particle collides with another deposited particle and to handle the outcome of the collision. Interestingly, this leads us to evolve beyond the classical Boltzmann-like treatment of collisions since the velocities of such colliding particles are naturally correlated by the same underlying fluid driving force (this point is discussed in more details in [9, Section 10.6]). In addition, it is believed that these approaches can be extended to address multilayer deposits (where inter-particle collisions and collective motion are much more frequent), though such extensions are still at the early stages of development.
<figure><img src="image_45.png" /><figcaption>Figure 34: Validation of a dynamic PDF approach showing the comparison between experimental data and numerical predictions for the adhesion force distribution (left) and remaining fraction on the surface (right).</figcaption></figure>

**Current use** By construction, dynamic PDF approaches are consistent with PDF approaches used to simulate particle transport. They are thus usually coupled with RANS simulation of the fluid flow which provide one-point information on the two first moments of the velocity field, i.e., the mean velocity field and the covariance tensor of the fluctuating velocity components [252, 251]. Special numerical schemes are used to integrate the governing SDEs [245, 83] and statistics of interest are obtained by Monte Carlo methods.

When coupled to proper models for the adhesion between rough surfaces (usually based on hemispherical asperities placed randomly on surface), these methods yield accurate predictions of both the adhesion force distribution (see Fig. 34a) and statistics on particle resuspension (like the remaining fraction displayed in Fig. 34b) to be obtained.

More recently, dynamic PDF approaches have also been coupled to DNS simulations of an impinging jet flow [253]. This allowed to investigate where particles are statistically more likely to be set in motion and to detach. As expected, the large majority of particle incipient motion and detachment events were found to occur in a region close to jet axis but less so as the distance from the jet axis increased (except just underneath the jet axis where the lateral velocity is too small to trigger such events).

### 4.3.3. Integrated kinetic approaches

**General description** Integrated kinetic approaches consist in estimating each resuspension event without computing particle dynamics on surfaces. Hence, these approaches simplify the resuspension process to a kinetic one in which resuspension is assumed to be equivalent to particle detachment, thereby disregarding particle motion along walls. These models provide directly information on the probability of resuspension for each individual particle as well as on statistical quantities (like resuspension rates, threshold velocities and/or critical velocities), but without any account of particle motion prior to detachment.

**Model formulation** Two main integrated approaches have been proposed in the literature.

* Kinetic PDF approaches evaluate directly the probability for a particle to be resuspended from a monolayer deposit by considering only one specific scenario for particle resuspension. Drawing on the simple picture where resuspension results from a balance between motion-inducing and motion-
preventing forces/torques, one can express the instantaneous resuspension rate $T_r(t)$ by:

$$ T_r(t) = \lim_{\Delta t \to 0} \left[ \iint_{\text{resusp. cond.}} p_{\text{ind.}}(s; \mathbf{F}_{\text{ind.}}) p_{\text{prev.}}(s; \mathbf{F}_{\text{prev.}}) d\mathbf{F}_{\text{ind.}} d\mathbf{F}_{\text{prev.}} \right. \\ \left. \qquad |\mathbf{X}_p(t) \cdot \mathbf{e}_{\perp} = 0; t \le s \le t + \Delta t \right] \quad (56) $$

where $p_{\text{ind.}}$ (resp. $p_{\text{prev.}}$) corresponds to the probability to have a certain value of the forces inducing motion $\mathbf{F}_{\text{ind.}}$ (resp. of the forces preventing motion $\mathbf{F}_{\text{prev.}}$). Analytic formulas for $T_r(t)$ can be worked out provided that the distribution of forces acting on the particle is known and provided that a resuspension condition is expressed.

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

(a) Sketch of the resuspension scenario used in the Reeks, Reed and Hall model: energy accumulates (until detachment) as a particle oscillates vertically due to turbulent fluctuations and to surface deformations. Reprinted with permission from [10]. Copyright 2014, Elsevier.

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

(b) Sketch of the resuspension scenario used in the Rock'n'Roll model: a deposited particle oscillates around a pivot point due to turbulent fluctuations until the contact with the pivot point is broken. Reprinted with permission from [10]. Copyright 2014, Elsevier.

<caption>Figure 35: Illustrations of the scenarios used in kinetic PDF approaches: Reeks Reed and Hall model (on the left) and Rock'n'Roll model (on the right).</caption>

To exemplify how this general formalism is put in practice and to illustrate such kinetic PDF approaches, it is instructive to briefly recall the main aspects of the RnR model (more details can be found in [24, 10]). The RnR model assumes that particles are set in motion by their oscillations around a pivot point due to rolling motion, which are induced by the drag and lift forces due to the interaction with the near-wall turbulent flow. In that context, the angle of deflection $\theta$ for small oscillations around the pivot point P (see Fig. 35b is expressed by a damped oscillator model:

$$ \frac{d^2\theta}{dt^2} + \beta_{RnR} \frac{d\theta}{dt} + \omega_{RnR}\theta = I_p^{-1}\Gamma'(t) \quad (57) $$

with $I_p$ the particle moment of inertia, $\omega_{RnR}$ the oscillation frequency, $\beta_{RnR}$ the damping and $\Gamma'$ the fluctuating component of the torque $\Gamma$ exerted on particles. Assuming that the particle is exposed to drag, lift, gravity and adhesion forces, one can write the torque as [10]:

$$ \Gamma = \frac{a}{2} F_{\text{lift}} + r_p F_{\text{drag}} - a F_{\text{adh}} - \frac{a}{2} m \mathbf{g} \cdot \mathbf{n}_{\perp} + r_p m \mathbf{g} \cdot \mathbf{n}_{\|} \quad (58) $$

where $a$ is the lever-arm distance between the pivot point P and the point at which adhesion forces act. The resuspension rate can then be obtained by evaluating the number of particles (per unit time) which reach an extremum (maximum) of the potential energy when the wall-normal fluid velocity is positive (i.e., oriented out of the surface). Assuming a Gaussian distribution, this gives the following expression for the resuspension rate $\tau_r$ (per unit time) of a particle exposed to a given adhesion force:

$$ k_r|_{F_{\text{adh}}} = \frac{\omega_0}{2\pi} \exp\left(-\frac{k(\langle \Gamma |_{F_{\text{adh}}} \rangle)^2}{\langle f^2 \rangle(1 + \eta_{RRH})}\right) \quad (59) $$
with $\omega_0$ a frequency (related to the oscillation frequency $\omega_{RnR}$) and $\eta_{RnR}$ a contribution to the potential energy from the resonant energy transfer (linked with the oscillation frequency $\omega_{RnR}$ and the damping term $\eta_{RnR}$). The term $f$ appearing in Eq. (59) corresponds to the fluctuating component of the hydrodynamic contributions $F(t) = F_{\text{lift}}/2 + r_p F_{\text{drag}}/a$. To obtain the unconditional resuspension rate $\tau_r$ (still expressed per unit time), one has to integrate this conditional rate over all possible adhesion forces (whose distribution is usually prescribed).

More generally, other scenarios have been considered: the model by Wen and Kasper [254] relies on a simple force-balance approach where a particle is detached from a surface as soon as the balance between hydrodynamic lift forces and adhesive forces is ruptured (i.e., $F_{\text{lift}} > F_{\text{adh}}$); the RRH model (proposed later by Reeks et al. [124]) considers an energy accumulation as particles oscillate vertically due to the surface deformations induced by turbulent fluctuations, which can ultimately induces resuspension when enough vibrational energy has been accumulated to overcome the potential well due to adhesive force (see Fig. 35a). Regardless of the resuspension condition considered, the resuspension rate obtained with kinetic PDF approaches is usually expressed in the form of a Boltzmann distribution, i.e., with an exponential factor containing the ratio between the potential well preventing motion and the energy of forces/torques favoring particle incipient motion.

* Impulse criterion models have been developed to evaluate the probability that a particle escapes from a multilayer deposit. As illustrated in Fig. 36, a particle deposited on top of a bed can start moving when the balance between motion-inducing and motion-preventing forces is ruptured. However, this condition is not sufficient for this particle to move freely from the deposit (i.e., either entrained by the fluid or migrating on the surface): this occurs only if the particle is able to escape from the pocket formed by the local deposit configuration. This means that resuspended particles have to acquire enough kinetic energy to overcome the potential barrier of the bed pocket. Hence, the following resuspension condition has been suggested: a particle is resuspended from a pocket bed when the forces acting on it are strong enough to induce motion and when these forces are acting over a long enough duration. This translates into an impulse condition, where the impulse $I_f$ has to exceed a given threshold value $I_{fc}$ [255]:

$$I_f = \int_{t_0}^{t_0+T} \|\mathbf{F}_{s \to p}(s)\| ds \geq I_{fc}. \qquad (60)$$

The hydrodynamic forces entering this equation have to exceed a certain threshold $F_c$, so that particles are set in motion. Expressions for this threshold value $I_{fc}$ have been derived later considering the competition between drag forces, lift forces or gravity [256, 257] and taking into account the slope of the river bed [258].

More recently, models based on an impulse criterion have been replaced by models based on an energy criterion [259]. The idea behind this change is that resuspension is more likely to happen when particles interact with high-intensity turbulent structures rather than with moderate-intensity ones. In fact, the energy transferred from the fluid to the particle is expected to be much higher when the driving hydrodynamic forces significantly exceeds the threshold force for incipient motion (i.e., $\|\mathbf{F}_{s \to p}\| \gg F_c$) compared to the one obtained with near-critical forces (i.e., $\|\mathbf{F}_{s \to p}\| \gtrsim F_c$). This has led to the introduction of a criterion based on the actual energy transferred to the particle [259]:

$$C_{\text{eff}} E_f = C_{\text{eff}} \int_{t_0}^{t_0+T} P_f(s) ds \geq W_c. \qquad (61)$$

where $W_c$ is the minimal amount of work required for complete particle entrainment, $P_f(t)$ is the instantaneous flow power (proportional to the cube of the local flow velocity), and $C_{\text{eff}}$ is the coefficient representing the efficiency of the energy transfer from the fluid to the particle. Again, expressions for the threshold values $W_c$ are available depending on the mode of motion (e.g., rolling or hopping [259]).
<figure><img src="image_48.png" /><figcaption>Figure 36: Sketch showing a pocket geometry: the red particle will be able to migrate on the surface only if the forces acting on it are strong enough to get the particle out of the pocket formed by other surrounding particles. Reprinted with permission from [3]. Copyright 2020, John Wiley and Sons.</figcaption></figure>

**Typical output and characteristics** As it follows from the above description, kinetic approaches yield estimates of the probability of a resuspension event. Yet, given that these models do not account for the migration dynamics, they can overpredict resuspension since some particles set in motion can nevertheless remain on the surface (and come to a stop at a downstream location [260]). In addition, kinetic PDF models have usually been designed considering only a single scenario (either sliding, rolling, lifting motion or even particle collisions). When resuspension occurs through a variety of processes (e.g., with particles rolling and sliding at the same time [261] or with mixed rolling/sliding and hopping motion as displayed in Fig. 9b), these formulations can become more complicated as they require additional information on the relative importance between each event in the overall resuspension process.

Another limitation is related to the input parameters required. In fact, these models require information on the forces/torques acting on the particle, which are often assumed to be described by normal or log-normal distributions (an extension to other distributions was given in [262]). Nevertheless, additional complexity arise when complex multimodal distributions occur (as has been measured recently for the adhesive force between a particle and a surface exposed to outdoor contamination [224]). Another critical parameter of the RnR model is the lever-arm of the adhesion and lift forces. The only experimental estimation of such a parameter was based on centrifuge measurements in [24] for a specific surface/particle pair and for only one particle diameter. The method allowed only for the determination of the mean value of the lever-arm and recent work has highlighted the RnR model sensitivity to the lever-arm parameter [237]. In addition, a basic hypothesis of the RnR model is that the two points of contact between the particle and the surface are symmetric with respect to the particle center. For an actual surface (e.g., like glass with nanometer roughness), a recent study also showed that the assumption that the points of contact are symmetric about the particle center is not always true, which implies that the lever-arm for the adhesion force may not be just twice that of the lift force [237].

**Current use** Due to their easy use and relatively low computational costs, kinetic PDF approaches have been widely used.

* Kinetic PDF approaches: These approaches have been extensively used in the context of atmospheric resuspension and especially in nuclear safety [23]. They provide quick access to the statistics on particle resuspension (including the remaining fraction or the time-averaged resuspension rate) while requiring information only on the distributions of adhesive and hydrodynamic forces. As a result, predictions of useful quantities for atmospheric resuspension can be obtained by measuring experimentally or
computing numerically the mean and standard deviation of these forces. For instance, as displayed in Fig. 37, the Rock'n'Roll model was shown to reproduce the fraction of 20 µm alumina particles remaining on a stainless steel substrate when the adhesion forces are assumed to follow a log-normal distribution (the best fit was obtained with a spread $\sigma_a = 10.4$) [24].

<figure><img src="image_49.png" /><figcaption>Figure 37: Comparison of experimental measurements and predictions with the Rock'n'Roll model for the remaining fraction of 20 µm alumina spheres on a stainless steel substrate exposed to a turbulent flow for a duration of 1 s. Reprinted with permission from [24]. Copyright 2001, Elsevier.</figcaption></figure>

Nevertheless, it should be remembered that characterizing distributions with information on their average value and their standard deviation is not always straightforward since these distributions can have very different shapes and do not necessarily belong to a unique class of PDF [187, 225, 224]: recent studies have indeed highlighted that adhesion forces with outdoor surfaces (contaminated by tiny objects such as dust or sand) display multimodal distributions that are very different from those obtained with clean laboratory surfaces [224].

In addition, since models are designed to capture either rolling, sliding or lifting motion, their range of validity is relatively well established. For instance, for a friction velocity of 0.1 m/s, the Rock'n'Roll model is applicable for particles smaller than 700 µm (resp. 50 µm) when they are suspended in the air (resp. in water). In ventilation ducts, the friction velocity is typically one order of magnitude larger, meaning that the limit of applicability of the RnR model drops by an order of magnitude to particles smaller than 70 µm. In the context of nuclear safety, where LOCA accidents might occur in very high temperature gas-cooled reactors at friction velocities up to 5 m/s, the RnR model will be applicable only to particles smaller than roughly 15 µm.

* Impulse criterion approaches:

These models have been developed in the context of fluvial transport. As a result, they are widely used to predict the Shields parameter for the resuspension of sediments from rivers [3]. Indeed, applying an impulse- or energy-criterion is especially adapted to the turbulent flows encountered in marine systems (with high values of the Reynolds number). This is not the case for static force- or torque-balance approaches, which are more convenient to treat low Reynolds-number flows (such as Stokes flows [263]). For instance, such a method has been used to predict the Shields parameter using a work-based criterion where a particle has to escape a pocket geometry (as the one displayed in Fig. 36) [264].
### 4.3.4. Force-balance/torque-balance approaches

**General description** Force-balance, or torque-balance, approaches consist in evaluating directly the resuspension of each particle by assuming that resuspension corresponds to the rupture of balance between forces/torques acting on the particle. Hence, they are similar to integrated kinetic approaches since they do not consider the dynamics of particles on the surface and also assume that resuspension is synonym of particle detachment. In a sense, they can be regarded as formulations developed in physical space while integrated kinetic PDF models are formulated in sample space.

**Model formulation** Since the resuspension criterion simply amounts to a rupture of the static equilibrium between forces inducing particle motion and forces hindering it, this gives the following deterministic conditions for each of the three resuspension modes (considering only drag, lift, gravity and adhesion forces):

$$ \text{Lift-off} \quad \mathbf{F}_{\text{lift}} > \mathbf{F}_{\text{grav}} + \mathbf{F}_{\text{adh}}, \qquad (62a) $$

$$ \text{Sliding} \quad \mathbf{F}_{\text{drag}} > \mu_{\text{fr}} (\mathbf{F}_{\text{adh}} + \mathbf{F}_{\text{grav}} - \mathbf{F}_{\text{lift}}), \qquad (62b) $$

$$ \text{Rolling} \quad \mathbf{M}_{\text{drag}} + \mathbf{M}_{\text{lift}} > \mathbf{M}_{\text{adh}} + \mathbf{M}_{\text{grav}} \qquad (62c) $$

with $\mu_{fr}$ the static friction coefficient. A closed set of equations is then obtained by assuming that the forces entering Eqs. (62) are given. For that purpose, one can either rely on measurement or on theoretical/empirical expressions for these forces (these expressions are described in Section 4.1).

**Typical output and characteristics** As for integrated kinetic approaches, force/torque balance models provide information on the probability of resuspension for each individual particle as well as to resuspension statistics (like resuspension rates, threshold velocities and/or critical velocities). Yet, no information on particle velocities upon resuspension or migration on the surface is available with such models.

One of their drawbacks lies in the fact that a specific mode of resuspension is often chosen right from the onset. This means that the expression obtained are only applicable to a certain size of particles. In fact, experimental observations have shown that small particles embedded within the viscous sublayer are more prone to rolling motion while larger particles, which are protruding from the viscous sublayer, are more susceptible to interact with near-wall turbulent structure and undergo direct lift-off [10].

**Current use** These models have been widely used, especially in the early stages of resuspension studies. Their popularity is probably related to the observations of threshold velocities for incipient motion, which inspired the idea that resuspension corresponds to a rupture of balance between motion-inducing forces and motion-preventing ones. The exact expressions of these conditions have been applied to a range of configurations, most of them related to the resuspension from monolayer deposits involving spherical particles interacting with smooth substrates [103] or rough substrates/beds [239]. These models were recently extended to compute the resuspension probability by accounting for the distributions in turbulent fluid forces or adhesive forces between rough surfaces. This has led to the use of such force/torque balance models in Monte Carlo approaches [265, 266, 104].

For instance, one model aims at capturing the effect of particle roughness and complex shape by assuming that uniform hemispherical bumps are present on the surface of spherical particles that are removed from the surface by rolling. The particle diameter, the bumps number, and their diameters are treated as random variables sampled from predefined distributions [265]. In particular, an experimentally-measured Gaussian distribution was used for the particle diameter, the number of bumps was assumed to follow a Poisson distribution, and the inverse of a roughness parameter was assumed to follow a log-normal distribution. The adhesion force was then computed using a JKR model (introduced in Section 4.1.4). Under these assumptions, both the adhesion and the aerodynamic forces become random variables. Hence, force-balance/moment-balance approaches can be used to describe single resuspension events for specific combinations of particle, surface, and fluid conditions but they can also capture statistical information on resuspension rates. By resorting to a Monte Carlo method based on torque-balance (for rolling motion),
this formulation was then used to estimate the resuspension rate as function of the fluid shear velocity, as shown in Fig. 38.

<figure><img src="image_50.png" /><figcaption>Figure 38: Comparison of experimental data from [103] and predictions obtained with a Monte-Carlo model based on a torque-balance approach (for rolling motion). The fraction of glass particles (with three sizes) is plotted as a function of the shear velocity. Reprinted with permission from [265]. Copyright 2012, Elsevier.</figcaption></figure>

### 4.3.5. System-scale models based on force/torque balance

**General description** Another category of models consists in predicting directly information at the system scale, i.e., at a macroscopic scale. This means that, instead of solving the resuspension and motion of individual particles, an approximation of the quantity of interest is derived directly from a combination of theoretical arguments and empirical formulas. In that sense, system-scale models are similar to the force-balance approaches presented before. Yet, the main difference concerns the predicted quantity: instead of estimating the average resuspension rate, these models provide information on more macroscopic quantities (such as the threshold velocity for incipient motion or the Shields parameter).

**Model formulation** As indicated, system-scale models consist in proposing an analytical expression for the selected quantity of interest. For instance, a recent review [122] has shown that various approximations of the threshold Shields parameter were obtained in the literature. These expressions vary depending on the theoretical considerations made (e.g., propositions based on force-balance for single rolling, sliding, and/or lifting motion, type of forces involved) as well as on the expressions retained for each of the forces (e.g., hydrodynamic drag force with/without shielding effects). To estimate the Shields parameter, these methods frequently rely on the use of empirical laws which relate a shear stress (or a friction velocity) to the average velocity in a wind-tunnel. The resulting expressions are easy to use and depend on a number of variables (such as the bed slope, the particle density, the particle-based Reynolds number). For instance, considering the resuspension of particles from a bed through sliding motion due to drag & lift & gravity forces and using the statistical theory of turbulence, the threshold Shields parameter $\Theta_c$ has been expressed as [267]:

$$ \Theta_c = \frac{4}{3} \cdot \frac{\mu_{fr}}{C_D + C_L \mu_{fr}} \times \left\{ \frac{10.08}{R_{\star,c}^{10/3}} + \left[ \frac{1}{\kappa} \ln \left( 1 + \frac{4.5 R_{\star,c}}{1 + 0.3 R_{\star,c}} \right) \right]^{-10/3} \right\}^{0.6} \quad (63) $$

with $\mu_{fr}$ the friction coefficient, $C_D$ the drag coefficient, $C_L$ the lift coefficient (more details in Section 4.1), $\kappa$ the von Kármán constant (usually set to $\kappa = 0.41$) and $R_{\star,c} = u_{\star}k_s/\nu_f$ the Reynolds number based on the friction velocity $u_{\star}$ and $k_s$ the bed roughness height (proportional to the particle diameter).
**Typical output and characteristics** The main advantage of such expressions is that they provide quick and relatively reliable information on an observable of interest. For that reason, they are often used as a source term for the amount of resuspended material in numerical simulations of particle transport. However, these expressions should be handled carefully since the range of validity of each model depends on the assumptions made to derive the expression (and possibly on the empirical formulas used for the forces entering the force/torque balance model).

<figure><img src="image_51.png" /><figcaption>Figure 39: Critical Shields parameter Θ_c as a function of the Reynolds number based on particle size and friction velocity R<sub>fc</sub>: a) Comparison between experimental data (yellow), the Shields curve (blue) and the model from Ali and Dey model (red); b) Comparison between various models. Reprinted with permission from [122]. Copyright 2018, AIP Publishing.</figcaption></figure>

**Current use** These models are widely used in practical applications such as fluvial or aeolian transport. In fact, various formulations have been proposed to express the Shields parameter as a function of various variables [122] (such as the particle size, density, etc.). A comparison between the various predictions is displayed in Fig. 39.

Another example concerns river morphology where morphodynamic models often rely on the Exner equation [64], which expresses the conservation of mass between deposited sediments within a bed and suspended sediments. This equation relates the changes in the bed elevation to the amount of sediment that are resuspended (i.e., the divergence of the sediment flux) and is usually formulated as:

$$ \frac{\partial z}{\partial t} = \frac{-1}{\rho_{\text{bulk}}} \frac{\partial q}{\partial x} \qquad (64) $$

with z the bed elevation, ρ<sub>bed</sub> the bulk density of particles, and q the sediment flux. This has led to the development of various expressions for this particle flux (interested readers are referred to dedicated reviews on this topic such as [64, 22]). Such formulations have also been used to predict the resuspension of sand
due to wind around obstacle (see for instance [149], which relies on RANS simulations of the fluid flow to estimate the local fluid velocity required as an input to these models).

### 4.3.6. Empirical formulas

**General description** Empirical formulas correspond to models based on empirical observations instead of relationships based on theoretical descriptions of the phenomena at play. These models are usually obtained by careful analysis of data coming from experimental observations or numerical experiments. From this analysis, simple formulas are used to express a relationship between the desired quantity (often called the observable) and different input variables known to affect this observable.

**Model formulation** When designing an empirical formula, the first step is to choose the observable of interest, and the second step is to select the input variables to be included in the formula. To be more specific, let us consider the case where we are interested in describing the time-averaged resuspension rate $\tau_r^{\Delta t}$ of aerosols in the atmosphere at short times (here less than 24 h). We further assume that the empirical formula should capture the evolution of the time-averaged resuspension rate $\tau_r^{\Delta t}$ as a function of the friction velocity $u_*$ and time $t$ only. The best fit obtained using experimental data from three experiments gives [268]:

$$ \tau_r^{\Delta t} = 0.01 \frac{(u_*)^{1.43}}{t^{1.03}}. \qquad (65) $$

Another fit can be derived if we also include information on the particle diameter $d_p$, the particle density $\rho_p$, and the average substrate roughness height $\langle r_a \rangle$. This gives [268]:

$$ \tau_r^{\Delta t} = 0.42 \frac{(u_*)^{2.13} d_p^{0.17}}{t^{0.92} \langle d_a \rangle^{0.32} \rho_p^{0.76}}. \qquad (66) $$

Experiments (real or numerical) used to derive an empirical equation must therefore have a form of internal consistency regarding fluid flow conditions (laminar vs. turbulent), particle size (small vs. large with respect to the boundary layer), particle shape (quasi-spherical vs. one dominant dimension, like in rod-like particles), surface type (bare vs. with vegetation), surface roughness (small vs. large compared to the particle diameter), wettability (of both the particle and the surface). Note that it is possible to capture multiple regimes with empirical formula, for instance with asymptotic equations, but this is usually outside their scope.

**Typical output and characteristics** The main interest of empirical formulas is that they provide quick and relatively reliable information on an observable of interest. For that reason, empirical formulas are often used in practical applications (e.g., dispersion of pollen or allergens). Empirical formulas can also be used in numerical simulations of particle resuspension and transport, mostly as source terms for the amount of resuspended particles during a given amount of time knowing the local fluid velocity.

As pointed out, their main drawback is their range of validity. In fact, these equations can only capture the physics and regimes that were present in the experiments used for model tuning, and they depend on the relative amount of data in each experiment as well as the estimation of model parameters that may have not been measured at the time of the experiment. Consequently, should we be interested in applying one such formula in a different context (where, for instance, the relative humidity of air plays a role), a new empirical formula would have to be designed to include the relative humidity as one of the input parameters (provided that information on this variable is available within the initial dataset). As a result, before applying an empirical formula to a specific problem, it is imperative to understand the range of applicability and the model limitations. For instance, the aforementioned empirical models fail to predict the short-term prompt resuspension which can account for the instantaneous aerosolization of a large portion of the deposit. This aspect is particularly important in the evaluation of health-effects.
**Current use** In the past, empirical formulas were based on dimensional analysis (see the review [130] and references therein). This has led to a number of linear regression or power laws for the time-averaged resuspension rate [268, 269]. Nowadays, new machine learning (ML) type approaches have emerged to derive data-driven models. While the traditional black-box approach is not very appealing due to its lack of transparency, new methods allow for the derivation of equations and the imposition of physics-based constraints. For instance, symbolic regression is a type of regression analysis that combines mathematical building blocks like mathematical operators and functions to derive the simplest equations that best fit a data set. Compared to other regression techniques, the user does not have to know the functional relationship among the variables a priori and the method uncovers both the model structure and the model parameters. Prior knowledge can be incorporated by using shape-constrained symbolic regression, in which an expected behavior can be prescribed (e.g., monotonicity of the function with respect to certain inputs). In general, these functions will be more complex than those obtained by dimensional analysis and power laws and have the potential to reveal more descriptive equations of the underlying physical system.

To illustrate the limitation of empirical approaches, various empirical models derived from various at-
mospheric resuspension measurements (given by Eqs. (65) and (66)) have been compared to additional
experimental data on the resuspension of lycopodium spores. As illustrated in Fig. 40, some of the models
are not able to reproduce the additional measurements due to the lack of accurate considerations of surface
properties (especially roughness) [268].

<figure><img src="image_52.png" /><figcaption>Figure 40: Evaluation of the performance of empirical models (see Eqs. (65) and (66)) against the data from Wu and Braaten experiment on the resuspension of lycopodium spores. Reprinted with permission from [268]. Copyright 2003, Elsevier.</figcaption></figure>

Nevertheless, due to their simplicity, these models have been extensively used and coupled to various approaches. For instance, empirical models are frequently applied to estimate particle resuspension from trees [270]: the principle is to express directly the fraction of deposit being resuspended as a function of the wind speed, the leaf area density, and particle concentration, which are then taken as source terms (i.e., new particles entering the physical domain) in CFD simulations of particle dispersion in atmospheric environments. Similar empirical expressions are introduced as source terms in CFD simulations of radionuclide dispersion (see for instance [271]). Expressions of the Shields parameter are also frequently used as source terms even in relatively high-level CFD simulations. For example, it has been used in combination to Large-Eddy Simulation (LES) of the fluid flow to study creep motion over a ripple [272] or even resuspension induced by helicopter hoover over a sandy surface [273].

4.3.7. *A hierarchical view on modeling approaches*

The various modeling approaches outlined in this section can be organized as a function of their information content. This is shown in Fig. 41 in which the horizontal arrow indicates a reduction of information in the treatment of the particle resuspension process. A few remarks can be made:
(a) Only two methods, namely the exact dynamical and the dynamical PDF, consider the complete resuspension process in the sense that they distinguish between detachment and actual resuspension. Both introduce specific effects for deposited particles, such as adhesive forces, but retain the same formulation as for particles transported by the fluid flow away from walls. Consequently, they are naturally consistent with corresponding simulations of particle transport in the fluid near-wall boundary layers. This makes them attractive candidates if one wishes to address the overall process in which deposition needs also to be considered (for instance, for saltating particles);

(b) In contrast, the kinetic and force-balance approaches appear as static methods since they assume that resuspension is equivalent to detachment and, therefore, to the rupture of an initial equilibrium, or static, state of deposited particles. Regardless of the differences in how adhesive and fluid entrainment forces are handled, these two models can be seen as equivalent, with one (the kinetic approach) being formulated in sample space and the second (the force-balance approach) written directly in physical space;

(c) The last two methods, the system-scale model and empirical formulas, bypass the treatment of the physical process at play in particle resuspension altogether and attempt directly to provide analytical expressions for either some key parameter or a time-averaged particle resuspension rate.

<figure><img src="image_53.png" /><figcaption>Figure 41: Sketch summarizing the various resuspension models and their level of description (arrow).</figcaption></figure>

## 5. The art of modeling: new insights into existing models

Two remarks can be drawn from the description of the physics involved in particle resuspension and current attempts at capturing its key aspects. First, there is no first-principle formulation acting as the reference theory and from which simplified versions would be derived. Second, there is no lack of models but rather of profusion of them. What is however puzzling is that not only do these models operate at different levels of description but they involve sometimes completely different treatments of the same physical phenomena. For instance, adhesion forces are regarded as input parameters in some models while they are considered as characteristic features to be predicted in others. This situation indicates that very different conceptions of a modeling approach are present but without being explicitly stated. To chart the modeling landscape, it is therefore useful to step back and address specifically what constitutes a model.

In practice, the formulation of a model implies building a consistent theoretical framework under known hypotheses which have direct impact on the quantities being predicted and the range of applicability. Some of these choices follow from theoretical considerations (like choosing the relevant mechanisms and forces at play), some from the applications which are studied (with the expected range of validity and applicability), while others refer to practical issues (e.g., the tractability of a model in real-life situations or its desired efficiency). Thus, modeling requires to combine and balance widely different aspects often manifested as
constraints to meet. In that sense, developing a model is akin to an art and it is henceforth referred to as
the “art of modeling”.

The objective of this section is to introduce, or recall, the underlying principles which makes up the reading grid with which the various modeling approaches presented in Section 4.3 will be assessed in Section 6. To do so, we consider the following questions:

(1) What is a model?

(2) How is a model conceived?

(3) How is a model used and evaluated?

5.1. *Defining a model*

Etymologically, the word model is derived from Latin *modulus* meaning small measure. In the field of physical sciences, a model corresponds to “a system of postulates, data, and inferences presented as a mathematical description of an entity or state of affairs”. By extension, it can also refer to “a computer simulation based on such a system”<fnref n="1" />. This definition shows that there are several notions that need to be addressed separately:

* The model, which corresponds to a simplified representation used to explain the behavior of a physical system or the workings of a physical process.

* Its implementation in a software, which involves the development/use of dedicated numerical algorithms to solve the set of equations defining the model (or provide accurate approximations).

* The realization of one or several numerical simulations using the model implemented in a software together with the analysis of the results obtained.

In this article, we essentially focus our attention on the first notion, i.e., the conception of a model as a description aiming at reproducing some key features of a physical system of interest. The numerical implementation of such models is left out of the scope of the present paper because it relates more to the fields of computational physics, numerics, and computer sciences (interested readers are referred to publications dedicated to this topic, such as [274, 275]). This does not mean that the development of accurate and efficient algorithms should be overlooked, as it has direct consequences on the precision of the results obtained and on their usefulness to the scientific community.

5.2. *Conceiving a model*

Modeling refers to the action of building models. However, beneath this deceitfully trivial definition lies a much richer and complex endeavor. For instance, Bonate mentioned in his introduction to the book entitled *The art of modeling* [276]: “there is more to modeling than the act of modeling. To be an effective modeler requires understanding the different types of models and when one type of model is more appropriate than another, how to select a model from a family of similar models, how to evaluate a model’s goodness of fit, how to present modeling results both verbally and graphically to others, and to do all these tasks within an ethical framework”.

To put it differently, we suggest here that modeling requires to address the following aspects, which are indicated in Fig. 42 (some of these notions are discussed in [9, Section 3]).

* The context and objectives of a study are one of the main factors influencing the choice of a model. Indeed, once the specific problem at hand is outlined, one has to start by identifying the quantities of interest that are relevant to fulfill the selected objectives. Given that it is unfortunately often overlooked, it is worth emphasizing that this step is an essential one in the construction of a model.

<footnote n="1">¹Merriam-Webster. (n.d.). Model. In Merriam-Webster.com dictionary. Retrieved in June, 2022, from https://www.merriam-webster.com/dictionary/model (see sense 3 and etymology)</footnote>
<figure><img src="image_54.png" /><figcaption>Figure 42: Sketch summarizing the various notions involved when conceiving a model including: aspects related to the context, theoretical considerations and specific constraints.</figcaption></figure>

For instance, in the context of multiphase flows, one could be interested in various aspects, including: identifying where particles are coming from (especially to distinguish between natural and anthropogenic sources), characterizing these particles (particularly their nature and toxicity for living organisms), measuring their dynamics (with fluxes and rates of transport, deposition and/or resuspension), or examining their fate (such as the location where they accumulate in the environment or if they can be inhaled by humans). The choice of a quantity of interest is directly derived from such questions: chemical analysis of the particles is paramount to characterize their composition, while measurements of rates/fluxes is more relevant to investigate the dynamics of particles in a flow or their fate (e.g. by placing sensors at specific locations, such as the average human breathing height in medical environments).

Once the factors of interest are determined, the minimum information required to derive these factors of interest needs to be clarified. For example, in the frame of statistical formulations of turbulent reactive single-phase flows, the objective is usually to obtain the average value of chemical source terms. As a result, the minimum information corresponds to the one-point one-time PDF of the concentration of scalars involved in the chemical reactions. This example is developed in [9, Section 3.2] in more details. What should be remembered is that the models thereby developed should provide at least this minimum amount of information. Of course, more complete models providing additional information are still viable (strictly speaking, this extra information is not required with respect to the objectives).

Another item to specify is the range of applicability. For instance, a model designed to capture the transport of sand in the atmospheric boundary layer (where canopy effects are important) cannot usually be directly applied to reproduce the motion of plankton in the ocean. In fact, it requires several adaptations: first, the model has to be adjusted/customized to capture the key features of the fluid flow oceans (possibly with stratification effects as those described in [67]); second, the model needs to be extended to treat the motion of plankton (which can be self-propelled and even display collective motion [277]). It is therefore important to define right from the outset whether a model is intended only for one specific application or for a range of situations.
* Having defined the context and objectives, the next step consists in analyzing the physics involved to identify the key mechanisms and the fundamental physical phenomena at play. For particle resuspension, these issues were already touched upon in Section 2.2. As indicated in Fig. 42, this provides answers to the following questions: What are the essential phenomena? What are the phenomena that can be neglected? And what are the phenomena that can be approximately accounted for but without too many details? These answers can be used as guidelines for a model construction since we look for approaches that allow to treat the essential phenomena without approximation while still containing, or resolving, enough information to model those phenomena that can be approximated.

For instance, if someone is interested in the motion of large suspended particles, it is reasonable to neglect Brownian motion while gravity must be properly accounted for. However, this means that the model cannot be directly applied to study the motion of small and light colloid particles for which Brownian motion is paramount while gravity can be neglected. Depending on which type of particles is considered, the same physical force is then regarded as either essential or negligible. Note, however, that if the objective is to treat any group of particles with sizes ranging from colloids to sediments, both Brownian motion and gravity become essential phenomena and, therefore, must be duly included (even though they can have little impact for a specific subsets of particles). This is in line with the remark above on the chosen range of applicability of the model to develop/select. In a similar situation, if we consider particles embedded and carried by turbulent flows (think about aerosols emitted from a chimney), particle transport is obviously a key phenomenon. If, furthermore, we wish to handle without loss of accuracy any set of particles having different diameters or even chemical properties, then this clearly tips the scales in favor of particle tracking methods. More detailed discussions can be found in [8, 9] on distinctions between slow variables representing degrees of freedom that need to be treated without approximation and fast variables representing other degrees of freedom whose effects can be mimicked with stochastic models.

* The conception of a model can be further restrained by specific constraints, such as those related to the model precision, flexibility, tractability, and/or modularity. More specifically, let us consider again the case of particle transport by turbulent flows for which, as pointed out above, particle tracking approaches are attractive candidates. Yet, depending on the level of precision required, the choice of the expressions for the hydrodynamic forces and surface forces can vary. For instance, if particles are much larger than the Kolmogorov scale, one can choose between very fine computations of hydrodynamic forces (based on PR-DNS simulations) or approximate expressions (using only the fluid velocity from PP-DNS methods, see Section 4.2.1). Of course, this has consequences in terms of the level of information contained in each formulation. Yet, it also implies that the computational costs associated with these two models are not comparable since the level of information is not the same. If the first approach is chosen, this limits ourselves to academic situations involving simple geometries and low Reynolds-number flows whereas the second one is more tractable and applicable to a much wider range of practical and real-life situations (see a similar discussion in [9, section 2.4]).

A complementary issue is related to the extensibility of the model, which concerns how easily a model can be modified to account for additional phenomena and mechanisms. The issue is to ascertain whether a chosen formulation is able to accommodate future developments when more complex physics is considered. For instance, we may wish to have formulations that can be extended in a straightforward manner by adding extra variables attached to each particle (representing additional degrees of freedom, such as rotation or orientation, etc.) without having to revisit the whole construction when we go, for example, from spherical particles to particles with more complex shapes. This is indeed one of the main issue to be addressed later on in Section 7.

In short, it is seen that modeling requires a deep understanding of the issues at stake, the specific objectives, as well as the limitations behind current fundamentals/theories. To rephrase Bonate, one could say that a modeler has to contextualize all the notions required. This contextualization is a long-lasting process which requires patience and efforts to avoid over-simplifying or overlooking some of the information.
### 5.3. Evaluating a model

Once a model is formulated, its accuracy has to be assessed: this corresponds to the so-called validation phase of a model. For that purpose, the model is usually implemented in a software. Then, either analytic formulas derived from this model or numerical simulations obtained by running the corresponding software are compared with available experimental data. This three-step process is illustrated in Fig. 43, which shows that there is actually a cycle between this triptych (modeling, implementation, analysis). In fact, comparison of a model predictions to new experimental data can bring out that some features are not properly captured. This can be related either to new properties that have been unveiled by more recent observations or to properties already known but that were previously deemed unimportant. Regardless of the origin of such discrepancies between observations and model predictions, the result of this cycle is that models are not frozen objects but are constantly evolving.

<figure><img src="image_55.png" /><figcaption>Figure 43: Illustration of the life cycle of scientific models: models constantly evolve as they are further tested in new configurations/applications and, on singular occasions, a sudden shift in the paradigm can lead to the development of a whole new field/category of models.</figcaption></figure>

For instance, in the community of particle-laden flows, accurate predictions of the transport and deposition of particles have been a long-standing topic [97]. Yet, even a cursory look at the evolution of these models reveals how they have been continuously enriched. In fact, numerical simulations first took into account particle transport by a flow (often turbulent) while particle-wall interactions were considerably simplified or downright neglected. Indeed, particles reaching the wall were often considered either as destroyed (numerically speaking) or as irreversibly deposited on the surface [97]. This gave first practical evaluations of the amount of particles deposited in a pipe as a function of time. Later, new experimental measurements provided more precise information on the morphology of the multilayer deposits thus formed, as well as on their effect on the flow (up to the complete clogging of the pipe). In order to reproduce these phenomena, questions related to the modeling of particle-wall interactions quickly came up since they play a key role in the morphology and cohesion of the deposits formed. Deposition models were therefore extended with new descriptions taking into account these particle-wall (and similarly particle-particle) interactions. At the same time, progress in the field of surface science (especially with AFM techniques) highlighted the role played by surface roughness on inter-surface forces. As a result, models for surface forces were also further enriched to account for the role of surface roughness.

### 5.4. Summary

In this section, we have introduced the “art of modeling”, which corresponds to a set of guidelines to rely on when establishing a model. The list of criteria used in this analysis does not pretend to be the only one, since additional criteria could be selected and their ordering changed. We believe, however, that these criteria already constitute a detailed enough framework allowing for a comprehensive investigation. Applied to the resuspension models, they turn out to be helpful in revealing the complexity of the models themselves and how different physical notions are deeply entwined (this is detailed in the next Section 6).
These guidelines are intended to help scientists and engineers alike in developing future models or in choosing the appropriate model with respect to the context of their study. This is indeed a task that should not be overlooked, especially due to the richness of the phenomena at play in particle resuspension, which involve various disciplines and a range of spatio-temporal scales. As it follows from the previous discussions, the choice of a model requires to weigh the desired range of validity and applicability of the model as well as how accurately the key phenomena should be described (e.g., mechanisms and forces at play). Therefore, as also claimed by Ancey [65], modeling “does not rely solely on explicit knowledge - that is a collection of equations, rules, diagrams - but also on the implicit knowledge gained through personal experience and interaction with highly practiced peers”.

# 6. Analyzing resuspension models

Drawing on the guiding principles presented in Section 5, we propose now to shed new light on existing resuspension models. This requires to adopt a leading thread. Indeed, the steps sketched in Fig. 42 should not be regarded as marking a one-way street leading from one context to the details of the models. Most of the time, the construction of a model implies to address the various aspects discussed in Section 5.2 in a different order than the one sketched in Fig. 42.

To avoid a cumbersome presentation with too many entries, we choose here to concentrate on the physics involved (the mechanisms, fundamental interactions, forces, etc.) and to introduce considerations about various contexts in the course of the discussion. This can lead to a variety of reading grids depicted in Fig. 44: they differ by the order of the criteria used to assess modeling approaches. In the following, we rely on the one displayed on the left, with the different criteria ordered as:

6.1 - Identifying, selecting and modeling the relevant mechanisms and phenomena considered;

6.2 - Determining the fundamental interactions at play and choosing how they are described;

6.3 - Clarifying the required and resolved information (including its physical meaning and level of description);

6.4 - Determining if the model is coupled to other approaches, and assessing the consistency / coherency between these descriptions (including the level of information);

6.5 - Evaluating the uncertainty and sensitivity of the outputs, thereby claryfying the range of applicability and precision of the model.

<figure><img src="image_56.png" /><figcaption>Figure 44: Illustration of 3 different reading grids to analyze resuspension models. The grids differ mostly in the order through which the criteria are treated (with slight variations in the content of each criterion). The left one is used throughout Section 6.</figcaption></figure>
## 6.1. Accounting for resuspension mechanisms

The first criterion is related to the identification and selection of the relevant mechanisms. In particular, the objective is to answer the following questions: What are the main physical mechanisms? Which ones should be treated without approximation? Which ones can be approximated or even neglected?

**Selection of the resuspension scenario** Various mechanisms have been highlighted in Section 2.2.3, including: rolling, sliding and lifting motion as well as collision-induced events. These mechanisms are not necessarily exclusive since a particle can move through a combination of several of them (see Fig. 9b).

It is then tempting to sort existing models with respect to the mechanisms that are accounted for. However, this would not reflect the assumptions that have motivated such choices. For instance, it appears reasonable to neglect inter-particle collisions in the case of sparse monolayer deposits, i.e., with only a few particles on a large surface. Yet, this assumption is erroneous when applied to dense monolayer deposits (i.e., where the inter-particle distance is smaller or equal to the particle size) or to multilayer deposits. Hence, it is not surprising to find out that models established in the context of fluvial and aeolian transport account for inter-particle collisions [3] (e.g., impulse criterion and some system-scale models) since they are by essence designed to be applied to complex multilayer beds. Meanwhile, many models developed in the multiphase flow community aimed at reproducing the measured resuspension from sparse monolayer deposits. As a result, these models often neglect inter-particle collisions while focusing on rolling/sliding/lifting motion. Yet, these models often focus on a single mechanism: some account only for rolling motion (e.g., the dynamic PDF approach [252] or the Rock'n'Roll model [24]), others only for lifting motion (like the RRH model [124]), or even a combination of rolling/sliding/lifting (as in force- and torque-balance approaches [103]). Another noteworthy consequence of neglecting inter-particle collisions is that it allows to treat resuspension events as independent from one another.

<figure><img src="image_57.png" /><figcaption>Figure 45: Sketch summarizing the various notions involved when choosing the relevant mechanisms at play in a model for particle resuspension. They are sorted according to the type of model (generic or customized), the type of surface (sparse/dense monolayer or multilayer) and the type of particles. The consequences in terms of the key mechanisms at play are highlighted (red) together with the possible models that can be used (dark green).</figcaption></figure>

From this discussion (see Fig. 45), it appears that the main issue is whether the model should be universal or specialized/customized for a given task/application. This choice is not trivial and has many consequences in terms of the range of applicability and validity of a given model. For instance, general-purpose resuspension models present the advantage of working for any size of particles and different type of deposits (i.e., monolayer or multilayer). This is the case for particle tracking formulations in general, and especially for the recent dynamic N-body approach since it reproduces any type of motion for spherical particles within a complex sediment bed [215]. On the contrary, customized models are not necessarily
applicable to other context/applications. This is obvious for empirical formulations, which can be viewed as
ultimate customized models: empirical approaches can indeed only be used in a context that is representative
of the data used to generate the empirical equation.

**Description of the mechanisms** Having chosen the mechanisms to be included in the model, the next step consists in selecting the appropriate level of description. For the sake of clarity, let us consider the case where a customized model for the rolling motion of colloidal particles is developed. One has now to decide how to treat rolling motion in the model. In line with the discussions in Section 5, there are three options. First, rolling can be treated without approximation by resorting to a dynamic model, which explicitly tracks the motion of particles as they roll on the surface (like in the dynamic PDF approach [252]). Second, rolling can be modeled to capture only its effect on the resuspension by resorting, for example, to impulse criterion [259] or to kinetic PDF approaches [24]. Third, rolling can be simplified to a rupture of balance for the incipient motion, as done in static torque-balance approaches.

**Modularity of a model** Another noteworthy aspect is the flexibility of a model. In the context of the mechanisms at play in particle resuspension, a model is considered flexible if the mechanisms can be changed without having to rework on all the foundations of the model. In practice, adding/removing/changing one of the mechanisms considered is a relatively straightforward operation when using Lagrangian tracking formulations (either dynamic ones or static ones like force- and torque-balance models). It basically amounts to adding/removing/changing one of the equation of motion of the particles. It can turn out to be a little bit more complicated when using integrated kinetic approaches, since a specific formula is obtained considering a set of predefined mechanisms. Hence, a new formula can still be obtained when the set of mechanisms is changed, but it has to be rederived (if it has not yet been done). Similarly, formulas obtained with system-scale approaches can also be extended to account for a different set of mechanisms since they are based on force- and torque-balance approaches. Empirical models are the only category that does not allow any flexibility on the mechanisms involved since each expression is fully determined by the data set and the corresponding conditions.

## 6.2. Fundamental interactions at play

The second criterion is related to the selection of the relevant interactions. More precisely, the aim is to
provide answers to the following questions: What are the fundamental interactions in particle resuspension?
Which ones should be treated explicitly? Which ones are to be modeled? Which ones can be neglected?

To address these issues, we shed light on the reasoning that drives the choice of the relevant fundamental
interactions and on the subsequent selection of an appropriate description for each force.

**Selection of the fundamental interactions** Various forces can play a role in resuspension. As detailed previously in Section 4.1, this includes contributions from: interactions with external fields (like gravity), hydrodynamic forces (due to particle-fluid interactions), contact and non-contact forces (related to particle-surface or particle-particle forces).

As for the resuspension mechanisms, the issue is whether one aims at obtaining a general-purpose or a
customized model. To be more specific, let us imagine that we wish to develop a general-purpose model for
the resuspension of hard spherical particles with arbitrary properties (size, density, and composition), for
any type of deposit (i.e., monolayer or multilayer), and exposed to any kind of flow (such as gases or fluids,
laminar or turbulent). In that case, it is natural to take into account most (if not all) of the forces listed
in Section 4.1. On the contrary, customized models allow to neglect some of the forces. For example, in
the case of sparse monolayer deposits, inter-particle collisions effects as well as particle-particle interactions
can be neglected compared to hydrodynamic effects. Hence, when choosing whether to include a given force
within the equations of motion, the issue is related to the range of applicability of the model. Drawing on
this notion, we propose to review in the following the motivations that allow to neglect some of the forces,
instead of comparing each model separately.
* The most common external force included in resuspension models is gravity (since particles studied are usually deposited on surfaces within a given planet). The case of electromagnetic forces is easier to discern, since these forces are typically included in a model only when particles having electro-magnetic properties are actually exposed to an electromagnetic field.

* Among the various hydrodynamic forces, the most important ones are the drag and lift forces. These forces are always occurring for particles deposited on a surface due to the strong shear flow in the vicinity of the surface. Yet, it is possible to neglect one of the two forces depending on the case studied. In fact, models for monolayer deposits are often customized to capture either the resuspension of small colloidal particles through rolling motion (which is induced by drag forces) [218] or of large sediments through direct lift-off (with negligible drag forces) [124]. Apart from these specific cases, other hydrodynamic forces can also enter the particle equation of motion, such as: the added-mass force (which is negligible when particles are much heavier than the fluid) and Basset history forces (which are often neglected since precise expressions are rather complex). As already discussed, Brownian motion can also be included when studying the motion of small colloidal particles.

* Surface forces are often simplified to adhesive and friction forces. Indeed, other surface forces only occur in specific situations, making it relatively easy to know if these forces should be accounted for. To be more specific, capillary forces are expected to occur when dealing with hydrophilic particles exposed to an airflow with a high relative humidity and for low surface roughness. An exception may be when dealing with organic particles, like spores, whose outer layer can change properties in humid environments. Hence, most situations are handled with only two major contributions: adhesion forces and friction forces. Since friction is a force that resists the translational motion between two surfaces in contact, it has to be included when sliding is a dominant mechanism (and can be neglected if sliding is negligible). Meanwhile, adhesive forces are usually preventing the motion of small particles (like colloids) and have to be included in such cases. However, when dealing with large and heavy particles (such as gravels in rivers [259]), adhesion forces can be neglected since they are several orders of magnitude smaller than gravity.

At this point, it is also worth noting that the previous discussion applies to all resuspension models. While this seems obvious for dynamic models (which include a number of forces in the equation of particle motion), it remains valid for empirical formulas too. In fact, empirical formulas can be adapted to describe the evolution of a quantity of interest as a function of the forces involved (similarly to system-scale models, see [122]). This would lead to empirical expressions of the form: $\tau_r^{\Delta t} = G(F_{hydro}, F_{surf}, F_{ext})$. Yet, such empirical formulas would be more intricate to use since they would require information on the forces rather than (or complementary to) information on the fluid/particle/surface properties. Hence, in terms of applicability, it is harder to keep track of the range of validity of such formulas. A more complete discussion on the required/resolved information from each model will be given later in Section 6.3.

**Description of the selected interactions** Once the relevant fundamental interactions are selected, one has to decide how to treat each of these forces in the model. In practice, two choices can be made: a force can be treated without approximation or it can be modeled. In the following, we illustrate these issues for the hydrodynamic and surface forces, which pose the greatest modeling challenge (see also Fig. 46).

* Since hydrodynamic drag and lift forces are proportional to the local instantaneous fluid velocity around each particle, they fluctuate due to near-wall turbulence (see [204, 205]). As a result, models for drag and lift forces involve two questions: how is the force evaluated and how is the relative velocity estimated? As illustrated in Fig. 46, hydrodynamic forces can be obtained with two approaches. First, simple models consider that the distribution of hydrodynamic forces is predetermined right from the onset (e.g., using Gaussian distributions as in kinetic PDF approaches [24]). The second approach relies on the use of dedicated models to obtain the distribution of hydrodynamic forces as a result. This includes both fine-scale descriptions (based on surface integration) or approximate expressions (see Section 4.1.2). Yet, when approximate expressions are chosen, information on the instantaneous
<figure><img src="image_58.png" /><figcaption>Figure 46: Sketch summarizing the various notions involved when choosing the fundamental forces in models for particle resuspension. It is presented in terms of the various forces involved (hydrodynamic, surface, external, collisions) together with the expressions of these forces and the information required on additional variables that enter these expressions (when relevant).</figcaption></figure>

fluid velocity at the particle position is needed. As described in Section 4.2.1, this instantaneous flow velocity can be accessed using three methods: it can be computed directly (using DNS), it can be approximated using turbulence models (such as RANS/PDF models), or one can consider a predefined distribution of these velocities.

* Similarly, adhesive forces can fluctuate very rapidly due to surface roughness (see Section 4.2.2). Information on adhesion forces can be obtained by solving numerically specific models or by assuming predetermined distributions (see also Fig. 46):

- There are two main difficulties when trying to come up with a model for surface forces between rough surfaces: (1) how to compute the forces? and, (2) how to account for the presence of roughness? As mentioned in the analysis of existing models for contact forces (see Section 4.1.4), these forces can be evaluated either with contact mechanics theories (that include deformation effects) or with a Hamaker approach for van der Waals forces (i.e., neglecting deformations). Hence, the choice between these two descriptions of surface forces is driven by the particle properties (here, their ability to deform). Yet, rough surfaces are very chaotic in nature and they display a range of roughness sizes and non-trivial distributions on the surface (see Fig. 28). Nevertheless, methods for the precise calculation of the force between rough surfaces are available when information on the exact surface topology is accessible. When models are expected to remain tractable in complex realistic situations, roughness is approximated with known functions (e.g. resorting to Fourier transforms, fractal representations, or simple hemispherical asperities placed on a smooth plane, see Fig. 30). This choice results therefore from a careful balance between the precision required, the minimum amount of information required/available, and the computational costs
associated to the method.

- Models can be even further simplified by imposing the distribution of adhesion forces right from the onset (as in kinetic PDF approaches [24]). The interest of such models is that they do not require any sort of complex computation. But one of their main drawbacks is that information on adhesion force distributions has to be accessible and provided beforehand, either from experimental measurements or fine numerical simulations. This is an important point, which questions the universal nature of such formulations (see the discussion in Section 6.3).

**Flexibility of the model** In the context of fundamental interactions, model flexibility is related to inter-operability and modularity. Inter-operability means that a force can be replaced by a similar force without having to rework on all the foundations of the model. Modularity means that the force can be extended by another formulation (including possibly extra effects). As for the mechanisms, model flexibility is inherent to particle tracking approaches, which allows to change not only any of the fundamental interactions entering the equations of motion but also their expressions. Meanwhile, integrated kinetic approaches can be adapted to account for additional forces but also for different expressions of the forces entering these equations. The only category of models that suffer from a lack of flexibility is the empirical models (since they rely on a predetermined set of data).

### 6.3. Nature of resolved/required information

The third criterion consists in clarifying the required and resolved information within each model (including its physical meaning and its level of information). More precisely, this step aims at answering the following questions: What are the input variables required by a model? What are the quantities expected to be predicted by a model? What is the physical nature of these variables? What is the minimum level of information required for these variables?

By the term “variable”, we mean here a number of possible physical entities that enter the formulation of a model. As displayed in Fig. 47, this covers three types of variables: the properties related to fluid/surface/particle (e.g., the fluid density or the surface roughness), the values of the fundamental forces at play (like the drag force acting on a particle), and the quantities describing particle resuspension (such as the fraction remaining).

<figure><img src="image_59.png" /><figcaption>Figure 47: Illustration of how the various variables can be used either as required information or resolved information. These variables/physical entities are sorted according to three categories (properties related to fluid/particle/surface, fundamental forces or resuspension quantities).</figcaption></figure>

As illustrated in Fig. 47, all these variables can be viewed as forming a sort of “pool of possible variables” and, depending on the context of the study, one may decide which are the required/resolved quantities.
Here, we wish to emphasize that the choices made during the first two steps (namely the mechanisms and the fundamental forces) have profound consequences on the information predicted by a model and the information required to use it. To be more specific, let us consider the case of a dynamic model coupled to a PR-DNS (as in [215]). The use of an individual particle tracking approach implies that detailed information can be obtained about resuspension quantities (such as velocities of individual particles, surface fluxes/rates like the Shields parameter or even volume concentrations). In addition, such a model yields predictions on the forces at play, especially on hydrodynamic drag and lift forces which are here explicitly treated (note that no information is available on adhesion forces since they are neglected for such large particles but they can easily be included if the model is extended to treat small particles). Meanwhile, this model requires detailed information on particle properties (e.g., shape and size) and on fluid properties (density, viscosity, velocity).

As it transpires from Fig. 47, values of the fundamental interactions at play can also appear as required information in other models. One example is the kinetic PDF approaches (like the Rock'n'Roll model [24, 278]), which computes information on resuspension rates using predetermined distributions for adhesion and hydrodynamic forces. In practice, the RnR model usually relies on log-normal and Gaussian distributions for adhesion and hydrodynamic forces [24], respectively. More recently, the model was extended to account for other distributions (like a Weibull distribution for adhesion forces, as in [278]). However, recent measurements have shown that adhesion forces can display complex multimodal distributions when using realistic outdoor surfaces contaminated by dirt instead of clean laboratory surfaces [125]. Not only can the distribution of adhesion forces take complex shapes, but it is also highly dependent upon the characteristics of each surface. In other words, unless we limit ourselves to the same well-characterized surfaces, the distribution of adhesion forces cannot be regarded as following a universal form. This means that, if we expect to apply the same model to any kind of particle-surface situations, adhesion forces should be resolved by the chosen modeling approach. In that sense, a universal formulation based on kinetic PDF can only be obtained provided that the distribution of adhesion forces are computed (and not prescribed). In fact, the same remark holds for hydrodynamic forces in kinetic PDF approaches.

Another important notion is related to the level of information contained in a model inputs/outputs. In fact, the level of information contained in the required variables has consequences on the level of information which can be resolved for the predicted quantities. This is especially the case when treating the adhesion forces between rough surfaces (provided that these are relevant interactions). If surface roughness is described using reduced statistical information (like the rms roughness size and the average peak-to-peak distance), only statistical information on the adhesion forces can be predicted by the model. Hence, unlike with detailed representations of surface roughness, the adhesive force of an individual particle at a precise location over the surface cannot be accessed. While this seems obvious for the computation of adhesion forces, it can be more ambiguous in the overall resuspension models (which involves descriptions of roughness, hydrodynamic forces, etc.). For instance, even if hydrodynamic forces are treated without approximation, the use of such a statistical approach for roughness implies that only statistical information on resuspension quantities can be predicted. This aspect resurfaces in Section 6.4 about consistency issues.

## 6.4. Consistency issues

The fourth criterion consists in assessing the consistency of the overall resuspension model, especially when it is composed of multiple sub-models (e.g., with description of hydrodynamic forces and the corresponding turbulent flow together with description of adhesion forces and the corresponding surface roughness). In particular, the objective is to evaluate the coherency and compatibility of the expressions entering the model with respect to the following issues: What is the application domain considered? What is the expected range of applicability? What is the minimum level of information expected/required?

Due to the multidisciplinary and multiscale aspects involved in particle resuspension, it is hard to keep track of all the assumptions involved in each formulation. However, checking consistency between these formulations is paramount in order to assess the range of validity and applicability of a model. In fact, the precision of a model is fully determined by the minimum level of information contained across all the variables entering a model (see the previous section). We illustrate these issues with a few examples:
* First, resuspension models can be inconsistent in terms of the descriptions of the variables entering the model. For instance, let us consider the case of a model based on a direct calculation of the adhesion force using information on the exact surface topology. We further assume that this formulation is combined with approximate expressions for hydrodynamic forces (like Eq. (18)) where the fluid velocity is not fully resolved in the near-wall region (e.g., through RANS simulations based on so-called wall functions or, in other words, with no down-to-the-wall integration). This would lead to a discrepancy between the level of description of the adhesion and hydrodynamic forces. In fact, such a model would amount to solving exactly physical aspects at the nanoscopic scale (where adhesion forces matter) while simplifying other effects acting at the hydrodynamic scale (i.e., over millimeter-size scales). Hence, the error made at the hydrodynamic scale would completely offset the precision of the methods for adhesion forces. Strictly speaking, such adhesion models would only be consistent with fully-resolved information on the near-wall instantaneous fluid velocity (e.g., using PR-DNS as in [215]).

While such inconsistencies are relatively obvious since adhesion forces act at scales much smaller than hydrodynamic forces, the opposite is not necessarily straightforward. In fact, due to the wide separation of scales, it is tempting to simplify or even neglect adhesion forces compared to hydrodynamic forces. This would be reasonable when studying large suspended particles for which adhesion forces are much smaller than the gravity force. However, when studying smaller particles, it is preferable to have a coherent level of information for both forces, as well as for their corresponding required physical entities (near-wall fluid velocity and surface roughness). Such issues have been largely explored in the context of particle transport and deposition [97].

<figure><img src="image_60.png" /><figcaption>Figure 48: (a) Visualization of the wake of a rotor in ground effect. Reprinted with permission from [279]. Copyright 2010, Vertical Flight Society. (b) Instantaneous flow structures in an excited, round impinging jet with embedded azimuthal vortices. Isosurfaces of the second invariant of the velocity-gradient tensor are colored according to their distance from the bottom surface. This shows the complex flow induced by a rotor which directly affects re-entrainment of detached particles in the near-wall region. Reprinted with permission from [273]. Copyright 2017, Cambridge University Press.</figcaption></figure>

* Second, resuspension models are often coupled to other approaches to compute the fluid and particle dynamics. For instance, it would be consistent and coherent to couple a dynamic resuspension model to a PR-DNS combined to N-particle tracking (as in [215]). In fact, this amounts to using a N-body Lagrangian formulation to describe both dynamics and resuspension. In contrast, associating a dynamic resuspension models to a reduced description of a two-phase flow (like RANS) would be inconsistent, since the only information provided by the flow solver is the average one-point velocity field and not the instantaneous velocity field (this point is detailed in Section 6.4.1). The reverse also holds: resuspension models that require information on the fluid velocity should be coherent with the level of information available in the chosen approach for the fluid phase. For example, an empirical expression for the resuspension rate can be coupled to fine simulations of a turbulent flow (like LES) and to Lagrangian tracking of detached particles (as in [273]). As displayed in Fig. 48, this has allowed the authors to investigate the role of near-wall turbulent structures induced by helicopter rotors on particle dynamics in the near-wall region. While this does not pose any difficulty in terms of numerics,
it can be questionable in terms of physical descriptions. In fact, the choice of an empirical formula implies that the average number of resuspended particles is evaluated using only a one-point one-time information on the instantaneous velocity. Yet, LES does provides more detailed information about the instantaneous but filtered velocities, which is not used for the resuspension model (i.e. attached particles) but is taken into account as soon as they get detached from the surface.

### 6.4.1. Avoiding typical statistical pitfalls

Coupling inconsistent descriptions of the fluid and particle phases does not just limit the overall precision but can also result in deeply flawed formulations. This typically occurs when the nature of the required statistical information is not properly analyzed and when what can - or cannot - be extracted from a given statistical model is not well understood.

To exemplify such pitfalls, let us consider a small volume, say $\delta V(\mathbf{x})$, around a given point $\mathbf{x}$ in which $N(\mathbf{x})$ particles are located at a given time $t$. In the spirit of N-body approaches (like DEM), we may wish to determine all the possible particle-particle interactions between these $N(\mathbf{x})$ particles, including inter-particle collisions or whether one particle is rolling on another one (see Fig. 49). Such a direct treatment of particle-particle interactions requires, however, to know the relative position of each particle within the small volume. In a statistical description, this means that we must have access to the PDF $p(t; \mathbf{x}_p^{(1)}(t), \mathbf{x}_p^{(2)}(t), \dots, \mathbf{x}_p^{(N(\mathbf{x}))}(t))$ of the $N(\mathbf{x})$ particle locations, $\mathbf{x}_p^{(i)}(t)$ being the position at time $t$ of the particle labeled $(i)$ with $i = 1, \dots, N(\mathbf{x})$ such that $\mathbf{x}_p^{(i)}(t) \in \delta V(\mathbf{x})$. This information is only accessible from a N-particle PDF method.

However, if we are using a one-particle Lagrangian PDF approach, it is essential to be aware that we have only access to the one-particle PDF of particle locations, that is $p(t; \mathbf{x}^{(i)}(t)) (\forall i = 1, \dots, N)$, but not to the $N(\mathbf{x})$-PDF. Actually, by resorting to a locally homogeneity hypothesis, through which probabilistic expectation at point $\mathbf{x}$ is replaced by an ensemble averaging over the $N(\mathbf{x})$ particles located in the small volume $\delta V(\mathbf{x})$, these $N(\mathbf{x})$ particles correspond to $N(\mathbf{x})$ samples of the one-particle PDF $p(t; \mathbf{x})$. Since such one-particle PDF methods are developed in the spirit of mean-field formulations, this means that the only information that can be extracted corresponds to one-point statistics (like the particle volume fraction $\alpha_p(t, \mathbf{x})$ or the mean particle velocity $\langle \mathbf{U}_p \rangle(t, \mathbf{x})$). In a weak sense, all we can write for the particle position PDF is

$$p(t, \mathbf{x}) \approx \frac{1}{N} \sum_{i=1}^{N} \delta(\mathbf{x} - \mathbf{x}^{(i)}(t)), \qquad (67)$$

which, using the local homogeneity assumption mentioned above, translates into $p(t, \mathbf{x})\delta V(\mathbf{x}) \simeq N(\mathbf{x})/N$. Therefore, all the possible ways through which the $N(\mathbf{x})$ particle locations are distributed within the small volume $\delta V(\mathbf{x})$, such as the two ones displayed in Fig. 49, cannot be distinguished. In the context of one-particle PDF methods, they are all equivalent since they yield the same weak approximation of the one-particle PDF, and correspondingly of the one-point PDF (e.g., mean concentration, mean velocity). Consequently, all attempts at calculating direct particle-particle interactions between the subset of $N(\mathbf{x})$ particles are bound to be meaningless.

Note that the same reasoning indicates that in one-particle method, even if the number of particles making up a deposit is correctly predicted, we still cannot deduce exactly from the simulations the resulting morphology or structure of the deposit. This means that the different ways according to which deposited particles are organized (see, for instance, the left panel in Fig. 12 in Section 2.2.3) cannot be distinguished. Capturing deposit morphology is, of course, possible but requires additional hypothesis to be introduced and a specific model to be developed.

A similar pitfall appears if we try to account directly for actions at distance between particles when information is incomplete. To emphasize this point, let us now consider two small volumes, $\delta V(\mathbf{x}_1)$ around point $\mathbf{x}_1$ and $\delta V(\mathbf{x}_2)$ around point $\mathbf{x}_2$ which contain $N(\mathbf{x}_1)$ and $N(\mathbf{x}_2)$ particles, respectively. If we wish, for instance, to extract the spatial correlation between particle velocities (so as to deduce the length correlation or the energy spectrum) at points $\mathbf{x}_1$ and $\mathbf{x}_2$, we are then considering

$$\langle \mathbf{U}_p^{(1)} \mathbf{U}_p^{(2)} | (\mathbf{x}_p^{(1)} \in \delta V(\mathbf{x}_1); \mathbf{x}_p^{(2)} \in \delta V(\mathbf{x}_2)) \rangle = \int \mathbf{V}_p^{(1)} \mathbf{V}_p^{(2)} p(t; \mathbf{x}_1, \mathbf{V}_p^{(1)}, \mathbf{x}_2, \mathbf{V}_p^{(2)}) d\mathbf{V}_p^{(1)} d\mathbf{V}_p^{(2)}. \qquad (68)$$
<figure><img src="image_61.png" /><figcaption>Figure 49: Illustration of consistency issues, with the amount of information available in N-point PDF formulations (middle) and one-point PDF approaches (bottom).</figcaption></figure>

Similarly, if we consider two-body particle interactions deriving from a potential $U_p^{\text{int}}(\mathbf{r})$, then the total force acting on a particle located at $\mathbf{x}_1$ is

$$F_p(\mathbf{x}_1) = \int \nabla U_p^{\text{int}}(\mathbf{x}_p^{(2)} - \mathbf{x}_1) p(t; \mathbf{x}_1, \mathbf{x}_p^{(2)}) d\mathbf{x}_p^{(2)}, \quad (69)$$

which involves the two-particle location PDF (obtained as the marginal of the two-particle position-velocity PDF in Eq. (68)). Using the principles introduced in Section 5.2, this means that, if our objective is to handle without approximation quantities such as the ones in Eqs. (68)-(69), then the minimal information required is the two-particle PDF. However, if we have first selected a one-particle PDF model, this information has been averaged out and we have only access to the one-particle marginals. In practice, this means that we have $N(\mathbf{x}_1)$ and $N(\mathbf{x}_2)$ particles in the two small volumes $\delta V(\mathbf{x}_1)$ and $\delta V(\mathbf{x}_2)$, but for each particle labeled (i) in the first volume (with $i = 1, \dots, N(\mathbf{x}_1)$), we cannot determine anymore its corresponding partner (j) in the second volume (with $j = 1, \dots, N(\mathbf{x}_2)$). The $N(\mathbf{x}_1)$ samples can be used to derive mean quantities at point $\mathbf{x}_1$ by Monte Carlo estimates while the $N(\mathbf{x}_2)$ samples can be used to derive mean quantities at point $\mathbf{x}_2$ by the same procedure, but no correlation between the two particle sets can be extracted from the one-particle PDF simulation. To be able to do so implies to track not a large number of particles but a large number of particle pairs.

The conclusion is that, if the objective is to treat without approximation quantities involving particle-particle interactions at a distance or at contact, we must be using a N-particle PDF approach. This requires that the velocity field of the carrier fluid (the continuous phase) be exactly known. In consequence, this rules out typical turbulence models, such as the RANS and even the LES ones, since they correspond to one-point statistical formulations. Applying a DEM approach for the particle phase is therefore only consistent with performing a DNS for the fluid phase. Phrased differently, coupling directly a DEM formulation to a one-point turbulence model leads to a deeply flawed and ill-based overall formulation. This is, unfortunately, a rather common mistake due to a poor analysis of the statistical requirement in terms of available information.
### 6.5. Uncertainty quantification

The fifth (optional) criteria is related to the statistical analysis of the model. The aim is to provide answers to the following queries: what are the overall performances of a model? In particular, how sensitive are the results to the variables entering a model? What are the uncertainties on the required variables and their effect on the resolved variables?

<figure><img src="image_62.png" /><figcaption>Figure 50: Sketch summarizing the key features in the analysis of results obtained from simulations: by identifying the requested variables (together with the associated uncertainty) and the resolved variables entering a model, one can evaluate the uncertainty propagation.</figcaption></figure>

As displayed in Fig. 50, this step requires to have clarified all the notions introduced in the previous sections (see also Fig. 44), especially the identification of the required/resolved variables. Hence, if performed, this analysis should be carried out after choosing/establishing a model. It usually relies on existing statistical tools, which allow two kinds of examinations:

* **Sensitivity analysis (SA):**

Sensitivity analysis consists in assessing how the outputs of a given model respond to variations in its inputs. In practice, sensitivity analysis ranks the inputs with respect to their importance in the variability of the output. Hence, this allows to distinguish the requested variables that mostly influence the resolved quantities (often in non-linear/exponential ways) and those that are less relevant. In addition, sensitivity analysis can reveal interactions among input variables and correlations among output quantities.

* **Uncertainty quantification (UQ):**

Uncertainty quantification is the science of quantitative characterization and reduction of uncertainties in models (and more generally any computational or practical application). It tries to determine how likely certain outcomes are if some aspects of the system are not exactly known. For that purpose, one usually relies on uncertainty propagation to estimate how uncertainty in the input variables is propagated through a model to the resolved variables.

From the previous definitions, it appears that uncertainty analysis and sensitivity analysis are closely related to each other. In fact, a good practice is to run both of them in tandem.

In the following, we would like to insist on the interest for scientists to rely more on such statistical analyses, which have been seldom used in the context of particle resuspension [280, 281, 282]. For that purpose, we do not delve into the existing statistical tools developed for SA and UQ techniques (interested readers are referred to dedicated literature on this topic, such as [161, 283]). Instead, we wish to underline what can be gained from such analyses. In particular, these analyses can prove very helpful in:

* **Model calibration:** Since a large number of variables enter resuspension models, sensitivity analysis can facilitate the calibration stage by focusing on the most sensitive parameters. Without such knowledge on the ranking of input variables, the risk is to lose time calibrating non-sensitive parameters.

* **Model simplification:** SA and UQ techniques allow to identify redundant parts within a model or inputs that have no effect on the desired output quantity. Hence, such information is helpful to potentially remove such parts/inputs from existing models in order to simplify them.
* Model improvement: SA and UQ techniques can reveal previously unsuspected connections/correlations between experimental measures, model inputs, and model predictions [280]. This helps designing updated or higher-order models in which the depiction of the observed phenomena is improved. Such methods can even help inching towards unified models that actually bridge the gap between monolayer and multilayer deposits as well as between aeolian and fluvial resuspension (see also [284, 242]).

At this stage, it is worth pointing out that gaining such precious knowledge is not free. In fact, one of the drawbacks of statistical analyses is that they require extensive amounts of data, each dataset corresponding to a slight variation in the input variables. In addition, replicates of the same input conditions are needed to assess the intrinsic stochasticity on a given output (like the time-averaged resuspension rate $\tau_r^{\Delta t}$). Hence, SA and UQ are delicate and time-consuming tasks, especially when relying on long experimental measurements and/or on computationally expensive models [282].

This points to the importance of understanding first these tools carefully in order to make the best out of them, before actually applying them in practical situations. This relates to one comment in Gigerenzer's book entitled *Reckoning with risk: learning to live with uncertainty* [285], where he criticized that “data analysis is typically taught as a set of statistical rituals rather than a set of methods for statistical thinking”.

## 6.6. Summary

At first sight, the presentation of the guiding principles related to the ‘art of modeling’ in Section 5 may appear like a step into a domain more akin to the philosophy of sciences than physics itself. However, not only is it sometimes sound to revisit the basics behind modeling efforts, but addressing the complex maze of the modeling proposals introduced in Section 4 from such a standpoint has proven fruitful in several aspects. First, this has allowed to bring out a reading grid according to which these model proposals have been organized and assessed. Second, this has provided us with a general framework that is helpful when attempting to go beyond present model limitations, in particular, by keeping us clear of two false impressions.

Indeed, new researchers in particle resuspension might derive the impression that models are bound to be incomplete, one way or another. Yet, once present shortcomings are clearly brought out, this incompleteness can be regarded as an incentive to pursue research work so as to push back present limits. Another possible impression from the profusion of available models is the misleading belief that particle resuspension has been thoroughly investigated, leaving no areas on which to shed new light. As recalled in Section 5, models are, however, regularly evolving through the triptych cycle (modeling/implementation/analysis). Therefore, thanks to this analysis, the stage for a proper assessment of model formulation is set, thereby paving the way for future extensions/developments of these models. This is addressed in the following Section 7.

# 7. Paving the way for new models and experiments

The objective of this section is to introduce future works needed to improve our understanding of particle resuspension as well as its modeling. Although still relatively unexplored at the moment, we believe that two directions in particular have the potential to become avenues of progress. They concern:

i. More complete and realistic models accounting for particles with complex shapes (e.g., non-spherical, deformable, rough particles) and intricate surface forces (related to surface heterogeneities);

ii. More general-purpose and unified approaches bridging the gap between two separate existing class of models (dedicated either to monolayer or multilayer resuspension) and covering a larger range of applications (i.e., not necessarily customized for either fluvial or aeolian resuspension).

In the following, these two roads are described in more details by addressing, at the same time, new experimental measures and refined models accounting for these new phenomena.
7.1. *Beyond spherical particles: towards more realistic descriptions*

From the discussions on existing models in Sections 6.1 and 6.2, it appears that a complete description of the resuspension process relies not only on the knowledge of the mechanisms at play but also on the specific forces entering these formulations. In turn, information on particle properties (e.g., size, density) as well as on surface characteristics (like roughness, composition) is required for adhesion forces (when relevant), while information on the fluid flow is needed for hydrodynamic forces. This implies that more complete descriptions of particle resuspension in realistic situations require further investigations of two aspects detailed below: the properties of particles and those of surfaces.

7.1.1. Resuspension of particles of complex shapes

Up to now, a large majority of particle-tracking approaches have focused on the case of spherical particles. This might appear as a somewhat strong assumption since many particles are not spherical (like pollen, fibers, or sand). Nevertheless, it has allowed to analyze precisely the physical phenomena involved and to develop satisfactory models that accurately reproduce the key resuspension features of spherical particles. However, we need now to go beyond this first approximation and consider more realistic particle shapes. In particular, two categories of particles are worth exploring:

• Non-spherical rigid particles:

**Physical origin:** Particles in the environment often display complex shapes. For example, many biological agents have rod-shaped geometries [176] (like bacillus spores shown in Fig. 51a), while fibers are usually thin and elongated objects [286] (as seen in Fig. 51b). Other particles have irregular surface features [287], like sand grains shown in Fig. 51c or pollens with spiky surfaces visible in Fig. 28d.

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

(a) Electron micrographic image of
*Bacillus anthracis* (Pasteur) spores.
Reprinted with permission from [288].
Copyright 2006, John Wiley & Sons.

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

(b) Electron micrographic image of 6 mm
long Basalt-fiber reinforcing a soil.
Reprinted with permission from [289].
Copyright 2021, Springer Nature.

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

(c) Electron micrographic image of quartz
sand. Reprinted with permission from
[290]. Copyright 2021, Multidisciplinary
Digital Publishing Institute.

Figure 51: Images of particles with a non-trivial geometry, ranging from rod-like shapes (left), thin elongated objects( middle) to multifaceted particles (right).

**Role in resuspension:** Particle shape anisotropy can impact resuspension in the following ways:

- First, hydrodynamic and adhesion forces differ from the case of spherical particles. In particular, the strength of both forces depend on the orientation of the particle with respect to the underlying surface and to the fluid streamlines. For instance, as depicted in Fig 52, elongated particles undergo much lower adhesion forces when only the tip of the particle is in contact with the surface (left figure) compared to the adhesion force obtained when the particle has its major axis oriented parallel to the surface (middle figure). Meanwhile, hydrodynamic drag and lift forces are much lower in the case of a streamwise orientation (middle figure) than a wall-normal or tangential orientation (left or right figures). This comes from Eq. (18), which shows that hydrodynamic forces are directly proportional to the area exposed to the fluid flow.

- Second, as a result of these forces, incipient motion depends also on the orientation of the particle with respect to the surface and to the direction of the flow. In particular, the relative importance
<figure><img src="image_66.png" /><figcaption>Figure 52: Illustration of how forces and probable motion change with respect to a spheroid orientation in contact with a surface.</figcaption></figure>

between the resuspension modes (namely sliding, rolling, and lifting) varies with the particle orientation. In particular, it has been shown that sliding occurs sooner than rolling for hemispherical or irregular particles while rolling usually occurs sooner than sliding for spheres [291]. To give another example, let us consider the case of a spheroid (as the one displayed in Fig 52). When the spheroid is in the tangential orientation case (i.e. when the flow is perpendicular to the particle rotation axis), the rod is expected to have a higher chance of rolling on the surface compared to a sphere with a radius equal to the spheroid second axis. On the contrary, when the flow is perfectly aligned with the particle major axis (streamwise orientation), the rod is expected to be harder to move than a sphere with a radius equal to the rod thickness due to lower drag force and the fact that this orientation somewhat impedes rolling.

**Experimental needs:** This brings out new challenges in terms of measurements. In particular, new experiments are needed to determine the type of motion of such non-spherical particles, both at the onset of motion and after detachment. For that purpose, one might consider techniques that have recently proven useful to measure the rolling motion of large millimeter-size particles [129]: it consists in using transparent particles with tiny black spots drawn on their surfaces while monitoring their motion using optical detection methods (see Fig. 53a). This technique provides access to both translational and rotational motion of particles rolling on the surface (see Fig. 53b). Hence, it can deepen our understanding of the resuspension modes for particles with complex shapes while providing valuable information for the development and validation of models for such particles. A parameter to carefully control is the boundary layer depth, which can be modified by the choice of fluid (e.g., air vs. water), to capture both particles immersed in the boundary layer and particles sticking out of it.

**Modeling challenges:** Following the methodology introduced in the description of current models in Section 4.3, the first step consists in selecting the degrees of freedom that characterize non-spherical particles as mechanical objects. To illustrate the additional modeling challenges compared to the spherical-particle situation, it is interesting to start with what seems a simple case of anisotropic particles, namely rigid ellipsoids. Indeed, ellipsoidal particles correspond to basic geometrical forms obtained by deforming a sphere by means of affine transformations. They are characterized by their extension ($r_1, r_2, r_3$) in each of the three directions, leaving a number of possible shapes such as: spheroids (i.e., ellipsoids of revolution $r_2 = r_3$), prolate shapes (extended ellipsoids like rods or pencils $r_1 > \{r_2, r_3\}$), oblate shapes (flattened ellipsoids like disks $r_1 < \{r_2, r_3\}$), and Zingg ellipsoids ($r_2/r_1 = 2/3$ and $r_3/r_2 = 2/3$, representative of sediment shapes [216]). Let us even consider the specific case of rod-shaped particles, which corresponds to a prolate spheroid (as shown in Figs. 52 or 56). Such rod-shaped particles can be characterized by a single parameter, namely the spheroid aspect ratio $\lambda = r_{p,L}/r_{p,S}$, which measures the ratio between the semi-major length $r_{p,L}$ (i.e. along the symmetry axis) and the semi-minor length $r_{p,S}$. From the sketches in Fig. 52, it is however clear that the basic forces acting on a rod-shaped particle depend on its orientation since the projected areas are now function of the rod orientation with respect to the flow streamlines. In a kinetic description
<figure><img src="image_67.png" /><figcaption>Figure 53: Images showing how markers on the surface of a transparent particle can be used to extract information on the rolling motion and, as a result, on the translational and rotational velocities.</figcaption></figure>

(the extension to non-fully resolved flow fields is addressed later), this means that the particle state vector retained for spheres, $Z_p = (X_p, U_p)$, must be extended to account, at least, for the particle orientation $P_p$, so that the minimum choice is to have $Z_p = (X_p, U_p, P_p)$. The rod orientation evolves as a function of the particle rotational velocity $\Omega_p$, since we have

$$ \frac{d\mathbf{P}_p}{dt} = \mathbf{\Omega}_p \times \mathbf{P}_p. \qquad (70) $$

It is worth pointing out that keeping only the rod orientation implies that some (rotational) information is lost. Indeed, one can capture tumbling motions (i.e., rotation around directions orthogonal to $P_p$), but not spinning motion (i.e., rotation around the direction aligned with $P_p$). For this reason, a slightly extended particle state vector is obtained by retaining the particle rotational velocity, which means that we are now considering $Z_p = (X_p, U_p, \Omega_p)$. For this state vector, the evolution equations which govern the dynamics of rods are given by [292]:

$$ \frac{d\mathbf{X}_p}{dt} = \mathbf{U}_p, \qquad (71a) $$

$$ \frac{d\mathbf{U}_p}{dt} = \frac{\nu_f}{m_p} \mathbb{R}^{-1} \hat{\mathbb{K}} \mathbb{R} \cdot (\mathbf{U}_s - \mathbf{U}_p) + \mathbf{F}_{\text{non-drag}}, \qquad (71b) $$

$$ \frac{d(\mathbb{I}_p \cdot \hat{\mathbf{\Omega}}_p)}{dt} = -\hat{\mathbf{\Omega}}_p \times (\mathbb{I}_p \cdot \hat{\mathbf{\Omega}}_p) + \hat{\mathbf{T}}. \qquad (71c) $$

In the equation for the translational velocity, Eq. (71b), the drag force has been singled out in the first term on the rhs while $\mathbf{F}_{\text{non-drag}}$ refers to all other relevant forces (lift, adhesion, etc.). It is seen that the expression of the drag force for rod-shaped particles involves a non-isotropic resistance tensor $\mathbb{K}$ since this drag force depends on the spheroid orientation with respect to the fluid. This resistance tensor can be expressed by a purely diagonal tensor $\hat{\mathbb{K}}$ in the local coordinate system attached to each spheroid $(\hat{x}, \hat{y}, \hat{z})$ (see Fig. 56). The tensor $\mathbb{R}$ hence corresponds to the rotation tensor that allows to go from the global frame of reference to the local one. The exact expressions of this resistance tensor can be expressed in terms of the Euler parameters and are detailed elsewhere (as in [293, 294]). Meanwhile, the equation for the rotational velocity involves the rotational inertia tensor $\mathbb{I}_p$ and the torque $\hat{\mathbf{T}}$ of all forces acting on a spheroid (more details in the pioneering work of Jeffery [295]). Note
that, in order to keep the rotational inertia tensor constant in time, Eq. (71c) is not written in the original inertial frame of reference, as in Eq. (71b), but in the particle-attached frame of reference (this is represented, for instance, by the notation $\tilde{M}$ which indicates that the torque $M$ is to be expressed in the frame of reference shown in the right panel of Fig. 56). Regardless of the details of the Jeffery model, it is important to point out that $\tilde{M}$ is function of the fluid velocity gradients at the rod center of mass, which we can express by $\tilde{M} = \mathcal{F}(G_{ij})$ where $G_{ij} = \partial U_{f,i}/\partial x_j$ is the fluid velocity gradient. Therefore, the Jeffery description relies on both the fluid velocity seen $U_s$ and the velocity gradient seen $G_{ij}$ (note that while $U_s$ is a vector, $G$ is a second-order tensor).

As it transpires from Eqs. (71), models for the dynamics of rigid spheroids require to have proper expressions for the forces acting on the particles, among which $F_{non-drag}$ or expressions of the drag force that go beyond the limit of validity of the Stokes regime (low particle Reynolds numbers).

i - Hydrodynamic drag and lift forces: As mentioned in Section 4.1.2, hydrodynamic drag and lift forces can be computed with microscopic descriptions, based on a direct integration of the extra fluid stress tensor along the entire surface of the particle. However, when developing models that remain tractable in realistic applications, one has to rely on approximate expressions. In the case of rigid spheroids, the drag force can be correlated to the angle $\phi$ between the direction of the fluid and the spheroid main axis [296]:

$$C_D = C_{D,\phi=0^\circ} + (C_{D,\phi=90^\circ} - C_{D,\phi=0^\circ}) \sin^{a_0}\phi \\ \text{with } C_{D,\phi=0^\circ} = \frac{a_1}{\mathcal{R}e_{\mathrm{p}}^{a_2}} - \frac{a_3}{\mathcal{R}e_{\mathrm{p}}^{a_4}} \text{ and } C_{D,\phi=90^\circ} = \frac{a_5}{\mathcal{R}e_{\mathrm{p}}^{a_6}} - \frac{a_7}{\mathcal{R}e_{\mathrm{p}}^{a_8}} \quad (72)$$

where $a_i \in [0...8]$ are nine shape-dependent coefficients fitted to PR-DNS integrations of the drag force on spheroids (see [296, table 2]).

Similar expressions have been suggested for the lift force:

$$C_L = \left( \frac{b_1}{\mathcal{R}e_{\mathrm{p}}^{b_2}} + \frac{b_3}{\mathcal{R}e_{\mathrm{p}}^{b_4}} \right) (\sin\phi)^{b_5+b_6} \mathcal{R}e_{\mathrm{p}}^{b_7} (\cos\phi)^{b_8+b_9} \mathcal{R}e_{\mathrm{p}}^{b_1 0} \quad (73)$$

with $b_i \in [1...10]$ ten shape-dependent coefficients fitted to PR-DNS integrations of the drag force on spheroids (see [296, table 3]).

These expressions were tested and validated for different particle Reynolds numbers $1 \le \mathcal{R}e_{\mathrm{p}} \le 300$, several aspect ratios $\lambda = (1.25, 2.5, 5)$ and various angles $0 \le \phi \le 90$. To the authors knowledge, no generic expressions valid for a larger range of shapes/Reynolds number have been proposed.

ii - Adhesion forces: As mentioned in Section 4.1.4, the adhesion between rigid bodies is well described by van der Waals forces. Hence, it is possible to use a microscopic description based on SEI methods (to integrate molecular-molecular forces over the exact geometry). Yet, one can also resort to analytical formulas using the Derjaguin approximation (for short separation distances). This gives a formula for VDW potential energy between a spheroid and a sphere which is similar to Eq. (28) for two spheres, but with a different geometrical prefactor [176]:

$$E_{\text{VDW,Spheroid1-Sph2}} = -\frac{A_{\text{Ham}}}{6h} \frac{r_{\text{p1,L}} r_{\text{p,2}}}{\sqrt{\left(\frac{r_{\text{p1,L}}^2}{r_{\text{p1,S}}} + r_{\text{p,2}}\right) (r_{\text{p1,S}} + r_{\text{p,2}})}} \quad (74)$$

Using the same arguments on the prefactor, the formulas that include the effect of surface roughness can also be generalized for spheroidal particles (see [176]).

iii - Other forces: To the authors' knowledge, no general formulation has been proposed for the capillary force or for the electrostatic forces between a spheroid and a plate [176].
At present, models for rigid non-spherical particles are usually limited to ellipsoids and follow the Jeffery formulation, Eqs. (71). These geometries are indeed easy to describe with mathematical expressions, making them compatible with fine simulations of the underlying fluid flow. For instance, PR-DNS [297] or LES [298] were performed in the context of large spheroids in water. These simulations have already confirmed the fact that ellipsoidal particles are more prone to sliding motion (or mixed combined motion) than purely rolling motion [298]. With respect to the models considered in Section 4.3, this means that only the first category, namely the exact dynamical formulation, has started to be extended to address rigid non-spherical particles. In this framework, complete information on the fluid flow is accessible and the particle state vector is indeed limited to $Z_p = (X_p, U_p, \Omega_p)$. In contrast, the formulation of extended dynamic PDF models based on the Jeffery description raises considerable challenges since the particle state vector is likely to be considerably augmented to include both the fluid velocity as well as the velocity gradients seen by rod-shaped particles. First suggestions have been made [9, Section 10.2] but work remains to be done to clarify the minimum amount of information on the fluid velocity gradients which is needed and to formulate new stochastic models able to capture their key statistical characteristics for rod dynamics. At the other side of the model spectrum, empirical models require ad-hoc experiments and parameters describing non-sphericity will become input variables. For intermediate formulations like kinetic PDF models, non-sphericity could be captured in the input PDF of adhesion forces. The issue is then to prescribe the resuspension mechanism, or their relative importance, depending on the amount of particles for each orientation (if not all equally possible).

**Summary:** The key challenge for non-spherical rigid particles is to come up with general-purpose models that can be applied regardless of particle size. Indeed, particle anisotropy has always an impact on hydrodynamic and adhesion forces, however different for low- or large-inertia particles. Tracer particles (i.e., with a small inertia) were recently shown to follow the fluid streamlines and align their orientation with the velocity gradient [299]. As a result, their anisotropy might not severely change their resuspension. The same does not hold for inertial particles, whose dynamics is not aligned with the flow streamlines. This suggests that a universal model can only be achieved provided that, at the very least, the particle orientation is included in the description.

* Deformable and/or flexible objects

**Physical origin:** Apart from having complex shapes, particles can also deform. For example, plastic debris can undergo plastic deformations when trapped in a sediment bed, leading to changes in shape near the contact area. Meanwhile, elongated particles (like fibers [286] or polymers) as well as quasi-2D objects (such as leaves) can deform under the action of external forces, leading to stretching, bending or twisting (as illustrated in Fig. 54). Another example is related to biological particles, which can change their size depending on the local environment (like spores shrinking in low relative humidity [48]).

**Role in resuspension:** Resuspension of deformable particles can differ significantly from the one of rigid particles. This is due to a combination of two effects: first, plastic deformations of the surfaces in contact can make adhesive forces much more complex; second, the flexibility of elongated particles (like fibers or polymers) can lead to a partial contact between the fiber and the surface (the remaining parts forming arches above the surface, as shown also Fig. 55a).

**Experimental needs:** New measurements of the adhesion between flexible elongated particles and surfaces are required to characterize how adhesion between them depends on the fiber orientation, the fiber properties (flexibility), the contact configuration (i.e., the size of the regions actually in contact) and the pulling angle. In fact, a flexible polymer has a certain orientation on the surface which induces a higher resistance to lateral sliding/rolling motion along the major axis. By grafting the end of a polymer to an AFM and pulling it away from the surface [302], recent experimental measurements have confirmed that desorption changes according to the pulling angle (see Fig. 55b). Ideally, these techniques should be coupled with real-time optical measurements of the particle shape to record the contact configuration between the flexible particle and the surface. They should help designing fine
<figure><img src="image_68.png" /><figcaption>(a) Example of a 3D reconstruction of a fiber conformation from experimental data (inset showing the side and top view, the color reflects the curvature). Reprinted with permission from [300]. Copyright 2018, IOP Science.</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>(b) 3D reconstruction of a DAS10 leaf with digital longitudinal (i.e., from leaf base to leaf tip) and transverse cross sections shown as black lines (units are in µm). Reprinted with permission from [301]. Copyright 2014, Oxford University Press.</figcaption></figure>

<caption>Figure 54: Images of deformable particles (such as elongated 1D fibers or flat 2D surfaces).</caption>

models to capture such intricate adhesion regimes, using for instance discretization techniques where a flexible particle is represented as a chain of ellipsoids or other tractable shapes [303].

<figure><img src="image_70.png" /><figcaption>(a) Schematic of the experiment, with a polymer partially adsorbed on the surface and pulled by AFM tip.</figcaption></figure>

<figure><img src="image_71.png" /><figcaption>(b) Average desorption force values plotted as a function of the pulling angle.</figcaption></figure>

<caption>Figure 55: Adhesion of a deformable polymer on a surface. Reprinted with permission from [304]. Copyright 2018, American Chemical Society.</caption>

**Modeling challenges:** Apart from the extension of adhesion models to account for plastic deformations (briefly addressed in Section 4.1.4), the key challenge is to develop new models for deformable 3D particles. Following again the methodology introduced in the description of current models in Section 4.3, a two-step process if used: (a) the first step consists in defining a set of additional degrees of freedom to capture all the features/phenomena observed experimentally; (b) the second step consists in developing the corresponding evolution laws for these extra particle-attached variables. To illustrate the additional modeling challenges compared to rigid anisotropic particles, it is interesting to start with what seems to be a simple case of 3D deformable objects, namely flexible elongated fibers. Compared to rod-shaped particles, a flexible fiber is a thin and elongated object that can hardly be described by the sole end-to-end distance (as displayed in Fig. 56). Instead, depending on the forces acting on the fiber, it can form complex sinuous shapes (somewhat like a snake). Describing such objects requires therefore to handle additional degrees of freedom compared to rigid rod-shaped particles. In the frame of particle tracking approaches, this would be directly implemented by extending the state vector and adding the corresponding equations. In practice, one can resort to three levels of description shown in Fig. 56:
<figure><img src="image_72.png" /><figcaption>Figure 56: Illustration of the possible descriptions of flexible elongated particles, from the most detailed one (on the left) to the reduced models (like chain-spring, dumbbell or ellipsoids, on the right).</figcaption></figure>

i - *Slender body theories* consist in describing a fiber as a slender elastic filament whose curvature is driven by the stresses acting on the fiber. In that case, its deformability can be described using beam theory, i.e. as an elongated cylindrical object whose deformation is given by the Euler-Bernoulli theory [305]. In the case of an inextensible and inertialess fiber in a viscous flow, the dynamics of a fiber of length $l$ and circular cross-section $a$ can be approximated using the following equation [306, 307, 308]:

$$ \partial_t \mathbf{X} = \mathbf{U}_f(\mathbf{X}, t) + \frac{1}{\alpha} \mathcal{D} [\partial_s(T \partial_s \mathbf{X}) - F \partial_s^4 \mathbf{X}], \quad \text{with } |\partial_s \mathbf{X}|^2 = 1 \quad (75) $$

$$ \text{where } \alpha = \frac{8\pi \rho_f \nu_f}{b} \quad \text{and} \quad \mathcal{D} = \mathbb{I} + \partial_s \mathbf{X} \partial_s \mathbf{X}^T, $$

with $s \mapsto \mathbf{X}(s, t)$ the curvilinear axis (parametrized by the arc-length coordinate $s \in [-\ell/2, \ell/2]$), $F$ the fiber's bending modulus (also called flexural rigidity, defined as $F = EI$ where $E$ is the Young modulus and $I$ the fiber moment of inertia), $b = -\ln(\ell^2/a^2e)$ the shape parameter and $\mathcal{D}$ the mobility matrix (whose form arises from the anisotropic drag exerted by the fluid on the slender fiber). This equation has to be complemented with boundary conditions at both ends of the fiber (free-end or fixed-end). This equation shows that such fibers deform due to the action of a velocity gradient along the fiber's length. Solving such equations numerically requires to discretize the curvilinear axis as an ensemble of inextensible and undeformable segments, i.e. with a length that is much smaller than the flexibility length (which measures the minimum length at which a fiber deforms under a given stress).

ii - *Chain-spring models* are more simplified models, where an elongated deformable fiber is represented by a chain of beads connected together through spring-like segments, or through chain-rod models that use rigid rods [9, 309, 310]. Hence, for each bead $\mathbf{X}_i$ composing the fiber, one has to solve the following equation:

$$ m_i \frac{d^2 \mathbf{X}_i}{dt^2} = -\zeta \left[ \frac{d\mathbf{X}_i}{dt} - \mathbf{U}_f(\mathbf{X}_i, t) \right] + \lambda_i (\mathbf{X}_i - \mathbf{X}_{i-1}) - \lambda_{i+1} (\mathbf{X}_{i+1} - \mathbf{X}_i), \quad (76) $$

where $\zeta$ denotes the individual drag coefficient of each bead $i$ and $\lambda_i$ is the tension between the $i$-th and the $(i-1)$-th beads. This tension can be described using simple spring models or fixed distance (as in Kramer chains). Such formulations allow to capture bending and twisting of elongated flexible particles with reduced costs. Even more reduced descriptions can be obtained using three beads connected with spring-like segments (trumbbells) or two connected beads (dumbbells).

iii - *Ellipsoids* can also be used as an extermely simplified representation of 3D elongated and flexible fibers. When dealing with inertia-less fibers with a size smaller than the Kolmogorov scale, it has indeed been shown recently that the dynamics of such objects is similar to the one of inertia-less
ellipsoids. This is due to the fact that they follow the fluid streamlines and fibers are always in a stretched state, except when interacting with rare but strong turbulent vortices [299]. Yet, by resorting to such reduced descriptions, details about the exact curvature of an elongated object are definitely lost.

Coming up with a reduced formulation for the orientation and deformation of a flexible fiber requires to clarify first the resolved information that is required. For instance, when one is interested in accurately reproducing the possible entanglement of fibers deposited on the wall, resorting to reduced descriptions based on ellipsoids or dumbbells is irrelevant since they only provide information on an end-to-end length and not on the curvature. In such cases, models should rely at least on formulations based on three or more degrees of freedom (i.e., trumbbells or any N-bbells with $N \ge 3$).

Things become even more complex when dealing with quasi-2D surfaces like leaves. In that case, few studies have explored their deformation due to the complexity of the fluid-structure interactions involved and the fact that each result can hardly be generalized to other geometries. For example, a recent study on 2D sheets of papers has shown that bending such surfaces can induce stronger stiffness than unbent sheets of paper [311]. In turn, material stiffness can affect particle resuspension. For instance, similarly to pole vault, a semi-flexible fiber can have only one of its ends in contact with the surface and be bent due to local shear. As the fiber resuspend, the energy that has built-up through elastic deformation can be converted into a wall-normal velocity, thereby leading to its motion away from the surface. In addition, tractable models require developing new reduced descriptions for the dynamics and deformation of such complex 2D surfaces, if relevant.

**Summary:** The key challenges for deformable 3D particles is to come up with general-purpose models that can be applied regardless of the particle shape and mechanical properties (flexibility, extensibility). This poses difficulties both in terms of experimental characterizations and numerical descriptions of the deformation of such particles across their whole surface (which can be a simple 2D line as in fibers or more complex 3D surfaces like leaves). Coming up with reduced formulations for the orientation and deformation of such complex 3D objects is thus paramount but requires first to clarify the resolved information that is needed.

### 7.1.2. Resuspension from complex surfaces

Current resuspension models already account for surface roughness. Nevertheless, several issues remain to be addressed and are described below.

**Experimental needs:** Current measurement techniques do not allow for simultaneous measurement of both adhesion forces and particle resuspension. This is due to the procedures currently used to measure the adhesion force: all existing methodologies rely on detecting the moment at which the contact between the particle and the surface is ruptured. In practice, various methods can be used to measure adhesion forces, which are summarized below and in Fig. 57:

* The centrifuge method consists in placing a surface with particles deposited on it in a centrifuge and measuring the velocity at which the particles are detached. Provided that the properties of each particle are known (size, density), one can reconstruct the detachment force using information on the velocity at which detachment occurs [312]. This can be used to measure either the wall-normal or tangential forces required to set particles in motion [43]. In addition, recent advancements have allowed to track the motion of particles as they move close to the surface (using continuous image acquisition [313]).

One of the advantages of this technique is that it can easily accommodate particles immersed in liquids. Hence, this method is popular in the pharmaceutical and food industry for evaluating the adhesive properties of polydisperse powders. Yet, its drawback is that it requires information on particle properties (mass, volume) and any uncertainty directly impacts the precision of the measurement.

* The drop impact test method consists in attaching particles upside down to a support stub inside a tube. The stub is then dropped from a given height and comes to a sudden stop when it hits a stopper.
<figure><img src="image_73.png" /><figcaption>Figure 57: Illustration of the existing techniques to measure adhesion forces. They are based on measuring particle detachment through various methods (such as centrifuge, drop impact, vibration, electric field detachment or AFM).</figcaption></figure>

Particle detaching from the surface fall through a hole in the stopper. The number of particles detached is determined by microscope image analysis of the surface before and after the fall or by weighting the sample. Video analysis of the impact provides impact velocity and contact time, from which adhesion force can be computed. By dropping the surface from different height, the fraction detached and the corresponding adhesion force can be determined.

This method has been used for instance in the food industry, to determine the adhesiveness of seasoning to tortilla chips [314]. This method has the same disadvantages as the centrifuge one.

* The *vibration method* consists in inducing vibrations at a given frequency and intensity until the particles deposited on a surface are detached. Particle detachment events are continuously monitored by image analysis and correlated with acceleration and surface displacement, which can be measured via a laser-scanning-vibrometer [315].

This method provides information on adhesion force distributions (see for instance [316]). Another advantage is that the apparatus can be enclosed in a humidity chamber to control relative humidity and, therefore, study capillary forces [317]. However, this method has drawbacks similar to the centrifuge one: the particle mass needs to be known to infer adhesion. In addition, the use of optical acquisition techniques restricts the applicability of this method to particles with a size close to 1 µm [315, 317]. Another important issue is that it is restricted to mechanically stable substrates [318].

* The *electric field detachment method* characterizes particle adhesion by using an electric field to remove particles from a surface placed between two parallel planar electrodes [319]. This method requires to charge particles using induction up to the saturation charge [320]. Once DC voltage is applied, particles move from the bottom electrode to the top electrode. Particle detachment is observed by monitoring the current flowing between the electrodes. The adhesion force is computed from the particle charge and the electrostatic field strength.

This technique allows to obtain the adhesion force distribution by using different applied voltages (as in [319]). In addition, the adhesion of particles of various sizes can be measured simultaneously. However, its main limitation is that it can only measure adhesion of conducting particles. In addition, the particle diameter is needed to compute the saturation charge.

* The *Atomic Force Microscopy* (AFM) is a form of scanning probe microscope that can be operated in various medium (air, gases, vacuum, or liquid) [321]. It can be used directly to measure surface
roughness (with a resolution of a couple of nm in the transverse direction and 0.1 nm in the wall-normal direction). Alternatively, a particle can be glued to the tip of a cantilever and, by bringing the probe particle in contact with a fixed surface and then retracting the probe, one can measure adhesion forces (simply by recording the deflection of the cantilever due to adhesion). To obtain statistical information on adhesion force distributions, the same procedure is repeated at different locations by moving the probe over the surface in a raster pattern.

The key advantage of AFM is that it allows to measure adhesion forces as small as $10^{-18}$ N [322]. In addition, the force is determined without any assumption concerning the particle size or shape. However, one of the limitations is that the technique can only determine the adhesion of one particle at a time, hence, it could be time-consuming for polydisperse particle collections or irregular particles. Another limitation is the fact that the direction of contact is predetermined by the way the particle is attached to the cantilever and the particle cannot adjust freely on the surface, like when real powders are dispersed on a surface.

This overview shows that, in a sense, measuring adhesion seems like measuring detachment events, except that the forces used to measure adhesion may be different from those that induce resuspension in the system of interest. For instance, surface vibration can be used to determine the adhesion between sand grains while, in practical applications, detachment is mostly induced by the wind. In addition, no information is provided on the re-entrainment of the particles once they detach from the surface. As a result, the two phenomena of detachment and resuspension are not measured concurrently. New experimental protocols or techniques might allow for simultaneous measurements in the near future. However, this might prove to be a very intricate tasks especially if plastic deformations occur (since the surface topology will be definitely changed after each contact/detachment series, making subsequent measurements unrelated to the previous one).

**Modeling challenges:** The absence of simultaneous measures of adhesion forces and resuspension events has profound consequences for the models of small colloidal particles. In fact, these models were validated by comparing only statistical information, like the time-averaged resuspension rate $\tau_r^{\Delta t}$, and not information on individual events, like the instantaneous resuspension rate of a single particle and its velocity. While models with a reduced description have proven useful to predict such averaged quantities (especially the median critical velocity $u_{r=50\%}$), more advanced formulations are needed to better capture rare events like resuspension at very low/high velocities [323] and the influence of seasonal variations [26]. This brings out two issues for reduced models:

*   **Refined models for rough surfaces:** To capture both short- and long-time resuspension rates, one cannot rely anymore on simplified descriptions of surface roughness (or adhesion forces). In fact, predictions of rare events (like resuspension at low or high velocities) can only be obtained provided that the tails in the distribution of adhesion and hydrodynamic forces are properly captured. This has implications in the current models for adhesion forces or surface roughness.

When the adhesion force is computed directly from a reduced description of surface roughness (see Section 4.2.2), roughness is often described with a statistical approach (like fractals, wavy sinusoids or random hemispherical asperities on a plate). In particular, approximate formulas for the adhesion force often rely on simple average parameters that are extracted from surface topologies (namely, the average roughness, the rms roughness, the peak-to-peak distance, or the highest/lowest peaks displayed in Fig. 58a). Yet, these parameters do not contain enough information on surface features, especially on the whole range of sizes and their corresponding density on the surface [218]. Since large roughness features are expected to induce higher adhesion forces than smaller ones and since they can provide shelter for particles in their wake, such large features can be responsible for the presence of particles that are very hard to resuspend (and only do so at high velocities). Hence, to go beyond current limitations, more systematic measurements of surface roughness are needed to describe the whole distribution of roughness sizes and the corresponding density on the surface.

When the adhesion forces are presumed (as in kinetic PDF approaches), the main difficulty comes from the use of prescribed laws (like Gaussian or log-normal distributions) to fit the measured distribution.
In fact, such fit are expected to filter the information on the tails of the distributions, possibly leading to erroneous predictions of rare events. This can be easily circumvented either by using more complex fit (e.g., based on multi-peak distributions, as shown recently in Fig. 29) or using directly the distribution of forces measured experimentally (provided it is detailed enough to contain relevant statistics on the extreme values of the adhesion forces).

<figure><img src="image_74.png" /><figcaption>(a) Sketch of a rough surface and the corresponding roughness parameters ($R_a$ and $R_{rms}$). Reprinted with permission from [218]. Copyright 2018, Elsevier.</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>(b) SEM images of a glass substrate contaminated by inorganic debris. Reprinted with permission from [224]. Copyright 2018, Elsevier.</figcaption></figure>

<caption>Figure 58: Illustrations of the complexity in describing rough surfaces: (a) average roughness parameters only provide limited information on roughness; (b) substrates can display very different roughness characteristics due to local inhomogeneities (e.g., when contaminated particles are present on the substrate).</caption>

*   **Models for inhomogeneous rough surfaces:** An additional difficulty to model real surfaces is that they can become contaminated by other objects that can bound to the surface. Such a case is displayed in Fig. 58b, which shows inorganic debris contaminating a glass substrate left outdoors for a certain duration [224]. Similar contamination can occur when studying the resuspension of biological particles from surfaces in outdoor environments [48]. The consequence of contamination is that surface roughness and adhesion forces can significantly change when a particle is in contact with regions of the surface that are not contaminated compared to regions that are contaminated. As a result, the adhesion force can display very complex multimodal distributions [224]. In addition, if contaminants are present on either the particle or the surface, more measurements are needed to characterize the spectrum of particle/surface interactions [224]. This implies more rigorous and time-consuming measurements so that the surface area used to “scan” the surface roughness and/or adhesion forces is large enough to be representative. To describe such surfaces, one can consider describing independently each region with uniform roughness properties, leading to the identification of small patches (i.e., isolated regions with irregular roughness features).

*   **Models for other surface heterogeneities:** Drawing on the previous comment, more realistic descriptions of particle resuspension should also include information on any type of surface heterogeneities. In the present paper, we have extensively discussed the role played by surface roughness, which corresponds to heterogeneities in the topology. However, surfaces can also display heterogeneities in their composition (as for composite materials like concrete) as well as in their chemical properties (e.g., charge). This can result in significant local variations in the resuspension of particles attached to such complex materials. To address such issues, the question is whether one is interested in capturing the variations of the resuspension rate on each local region (depending on its properties) or if only macroscopic information has to be resolved (e.g., on the average resuspension over the whole surface).

**Summary:** The limitations discussed here illustrate that improved understanding of the adhesion and resuspension processes are expected to emerge only from tight coupling between new experiments and fine-
scale simulations. This can then lead to more advanced models that are based on reduced descriptions of
these phenomena.

7.2. *Beyond individual resuspension events and collection motion: towards unified descriptions*

Broadly speaking, most of the models presented in Section 4.3 fall into two main but separate categories. On the one hand, many models have been designed to capture resuspension from sparse monolayer deposits where particles are hardly affected by the presence of other particles. On the other hand, some models have been formulated to address directly multilayer resuspension, therefore considering that a bed of particles is initially present on the surface and accounting for collective particles effects. Apart from the extensions related to more complex individual particle characteristics discussed in Section 7.1, it is important to address what takes place in between these two asymptotic cases and, more precisely, the new phenomena induced when particle-particle interactions interplay with particle-fluid and particle-surface forces. The development of unified approaches (i.e., descriptions that can span all the range of possible applications and types of deposits) hence requires to investigate two important aspects detailed below: the transition between individual/collective events and the role of collective motion on resuspension.

7.2.1. *Bridging the gap between monolayer and multilayer deposits*

Thanks to the use of dense monolayer deposits (i.e., deposits where the initial surface concentration of particles is high enough so that the average inter-particle distance becomes comparable to their diameter), studies on the influence of particle-particle interactions on resuspension have recently emerged [324, 111, 110]. These studies have revealed several phenomena that can help bridging the gap between the descriptions of sparse monolayer and multilayer resuspension:

• Collision propagation:

As displayed in Fig. 10b in Section 2.2.3, when the initial distance between deposited particles is small enough, collision propagation can occur due to the motion induced by a particle colliding with one of its neighbors. As a result, when dealing with dense monolayers, particle collision effects cannot be neglected as usually done in models for sparse monolayer deposits.

The role of particle collision in dense monolayers is two-fold:

- First, the collision frequency is directly driven by the initial particle layout on the surface. Indeed, it is obvious that two particles located next to each other collide if upstream particles are set in motion due to hydrodynamic actions. Hence, the initial presence of clusters of particles on the surface can lead to a large number of particles being set in motion, especially if one of the upstream particles starts moving, which is more probable since upstream particles are more prone to interact with a flow that is unaffected by the presence of other particles. On the contrary, a particle initially deposited further downstream is going to be hit by a moving particle only if it is within its actual path and provided that it is not too far (at which point, the moving particle might be too distant from the surface to interact with other deposited particles, see Fig. 59).

This means that new experimental studies should provide detailed information on the initial lay-
out of particles on the surface, their gaps and patches, especially to identify the regions where
the mean particle distance is lower than the average value derived by assuming a uniform particle
distribution over the surface (see also Fig. 59). For that purpose, two types of measurements are
relevant depending on the level of information required: first, the exact location of each particle
(using optical techniques and algorithm detecting the presence and position of each particles);
second, reduced statistical information, for example by resorting to radial distribution functions
(which can prove very helpful to characterize each region where particles are homogeneously dis-
tributed). These experimental data can then be used either in deterministic particle-tracking
approaches (where each particle is placed according to the exact positions measured) or in PDF
approaches (where only statistical information on the initial particle placement is required). Re-
gardless of the level of description chosen, the corresponding resuspension models should be
<figure><img src="image_76.png" /><figcaption>Figure 59: Sketch showing how inter-particle collisions can lead to multiple resuspension events depending on the initial layout of the particles: only few collisions are expected in regions with low surface concentrations while numerous collisions are expected in regions with high surface concentration (e.g., due to clustering).</figcaption></figure>

extended to include specific algorithms to detect or to account for inter-particle collision events: collisions can naturally be detected in N-body tracking approaches whereas statistical models have to be devised in PDF formulations. These statistical models can be built upon the notion of a “residence time” [110] which corresponds to the amount of time spent by a moving particle in the vicinity of the surface, where it can interact with its neighbors. As displayed in Fig. 59, this allows to determine the average area covered by a particle as it moves near the surface and, using the statistics on the particle layout, to estimate the collision probability with other particles within this area.

— Second, information on the outcome of each collision is needed. In fact, as displayed in Fig. 60, four scenarios can take place for the outcome of a binary collision depending on the motion of each particle after the collision: (i) only the impacting particle moves; (ii) only the initially fixed particle moves; (iii) both particles move; (iv) both particles are immobile. Drawing on existing results for the collision of particles embedded in a fluid [134], one can expect that the velocity and direction of motion of each particle after the collision depend on a number of physical variables, including the impact angle, the impact velocity, the particle properties (deformation) and the adhesive forces (especially those acting on the fixed particle which prevent motion).

New experiments should be carried out to investigate the outcome of inter-particle collisions, especially to identify the energy and momentum exchanges during such events. In turn, new models are expected to be developed in the near future to account for these possible outcomes (possibly including theories for plastic deformations which can lead to significant dissipation of energy during such collisions, when relevant).

In the previous paragraph, we have insisted on the role of collision propagation assuming that only bi-
nary collisions take place. However, when the initial surface concentration of particles is high enough,
more complex collisions involving three (ternary) or more particles can occur more frequently. There-
fore, future studies should investigate whether such complex collisions have an impact on collision
propagation and, when relevant, develop accurate models to reproduce their effect. For that purpose,
it might prove useful to get inspiration from existing works on avalanche [325] and more generally
granular media [326]. We will come back to the role of such collective motion later in Section 7.2.2.

• Complex surface motion:

Studies on dense monolayer deposits can also shed light on more complex surface motion. As we
have seen, descriptions of sparse monolayer deposits are typically based on three resuspension modes
<figure><img src="image_77.png" /><figcaption>Figure 60: Sketch showing the possible outcomes of each inter-particle collision event: the moving particle can either come to a halt or continue its motion, while the initially fixed particle can remain deposit or start moving. This leads to 4 possible outcomes (with variations in each one due to velocity differences for example).</figcaption></figure>

(namely rolling, sliding and lifting). However, as particles collide with each other, some of them can stop moving (especially when the impacted particle is tightly bound to the surface) [327, 118]. In multilayer resuspension, the cessation of motion has been studied due to its important role in sediment motion in rivers [75]. It was shown that, once particles start moving, the fluid velocity can be decreased below the threshold of incipient motion without resulting in an immediate halt of the particle. This is related to the inertia of sediments, which can continue their motion for a certain amount of time even when the flow stops.

Note that such complex motions with possible arrest are also possible in sparse monolayer resuspension. For example, rolling/sliding colloids can encounter varying adhesion forces as they migrate on the surface. As a result, if they reach a region with high enough adhesion forces/torques which lasts for a sufficient duration, these particles can eventually stop moving as displayed in the right panel of Fig. 61 (see for instance [328, 329, 330]). Unfortunately, there is at the moment little information available on the frequency and importance of such events on resuspension. There is however a difference in the physical interactions involved. In sparse monolayer deposits, complex particle motion along a wall surface is only due to the competition between hydrodynamic and adhesion forces. In denser monolayer deposits, particle-particle collisions start to play a role and this leads to a much more intricate interplay between hydrodynamic, adhesion, and collision effects.

As a result, more systematic characterization of particle dynamics near the surface is needed, including the distinction between their mode of motion (rolling/sliding) and potential arrest of the particles. This brings out interesting questions about migrating particles that come to a stop: can these particles start moving again later? Or do they remain immobile due to the local strong forces preventing their motion? Such information can be key to model long-term resuspension. In fact, as depicted in Fig. 61, if particles come to rest and reattach on the surface in regions with very strong adhesion forces (and/or possibly low hydrodynamic forces due to sheltering/shielding effects), these particles are likely to remain on the surface for a long time.

Hence, a noteworthy consequence of such complex motion is that it can lead to more pronounced extreme events: adhesion measurements might show the rare occurrence of strong adhesive forces but, due to this complex dynamics, many more particles can end up trapped in these regions. This is all the more important when dealing with dense monolayers where, once one particle is trapped in such a region, a cluster can quickly be formed due to particles colliding with this already deposited particle (that can act as an attractor for this region). Consequently, the need for more detailed investigations of
<figure><img src="image_78.png" /><figcaption>Figure 61: Illustration of the complex motion of particle migrating on a surface.</figcaption></figure>

adhesion forces on rough surfaces (mentioned previously in Section 7.1.2) should be supplemented with careful analysis of particle dynamics. This implies that previous measurements of particle resuspension by placing particles and simply counting the number of particles remaining on consecutive images after being exposed to a given flow are not sufficient anymore.

We believe that precise measurements of such particle dynamics are now within our grasp thanks to recent advances in terms of experimental techniques (especially to measure rolling motion [129], as displayed in Fig. 53). In fact, some insights have already been provided regarding the rotational/translational velocity of migrating particles as well as the time/distance traveled between their incipient motion and their detachment. For instance, micrometer-size particles were shown to be able to travel a distance roughly equal to hundred times their diameter before detaching from engineered surfaces [330]. Another noteworthy issue is that the distinction between detached particles and particles moving on a surface raises new questions regarding the accuracy of previous measurements of particle resuspension. Actually, previous experiments often relied on counting the number of particles remaining on a certain region of interest. Yet, there is no distinction made between detached particles and particles rolling/sliding out of the observation area. This would be especially true if the mean distance traveled by moving particles before actual detachment is large.

<figure><img src="image_79.png" /><figcaption>Figure 62: Simulations showing the relative importance of the migration (i.e., rolling, sliding), detachment and arrest as a function of the friction velocity. Simulations have been carried out using a recent dynamic PDF model [83] and considering the case of 40 µm particles on a rough substrate (25 nm asperities covering 2% of the surface) exposed to an airflow.</figcaption></figure>

As far as modeling approaches are concerned, such complex motion would be naturally accounted for by particle tracking approaches. To illustrate this aspect, a numerical simulation has been performed using the recent dynamic PDF approach for 40 µm spherical particles [83]. The model has been extended to account for sliding and lifting motion (arrest is captured when the adhesion forces/torques exceed the hydrodynamic forces). As seen in Fig. 62, the relative importance of these mechanisms varies with the fluid velocity. Migration arrest happens to reach its highest value at intermediate velocities. This can be interpreted as follows: the number of particles reaching halting regions initially increases with the fluid velocity since particles are able to travel longer distance (thereby increasing their probability to reach such regions); yet, as the fluid velocity goes beyond a certain threshold, the adhesion forces in