Royston Fernandes
email: royston.fernandes@inra.fr

Sylvain Dupont

Eric Lamaud

Investigating the role of deposition on the size distribution of near-surface dust flux during erosion events

Keywords: Dust flux, Size distribution, Dust emission, Dust deposition, Deposition velocity Dust flux, Size distribution, Dust emission, Dust deposition, Deposition velocity

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Mineral dust aerosols are a fundamental component of the Earth system. They refer to fine suspended soil particles (diameter < 20 µm) released by strong winds in regions with erodible soils [START_REF] Knippertz | Mineral Dust: A Key Player in the Earth System[END_REF]. The life-cycle of airborne dust involves three stages: (1) the release from the surface, (2) the transport in the atmosphere over long distances by turbulent diffusion and against gravitational sedimentation, and (3) the deposition back on to the surface. These three stages shape the particle size distribution (PSD) of the near-surface dust flux. Interest in understanding the PSD of dust flux is fueled by the impact of the aerosol size on different Earth processes (e.g., [START_REF] Mahowald | The size distribution of desert dust aerosols and its impact on the earth system[END_REF] such as weather and climate (e.g., [START_REF] Claquin | Radiative forcing of climate by ice-age atmospheric dust[END_REF][START_REF] Otto | Atmospheric radiative effects of an in situ measured saharan dust plume and the role of large particles[END_REF], ecosystem nutrient supply (e.g., [START_REF] Swap | Saharan dust in the amazon basin[END_REF][START_REF] Knippertz | Mineral Dust: A Key Player in the Earth System[END_REF], and human health (e.g., [START_REF] Weuve | Exposure to particulate air pollution and cognitive decline in older women[END_REF][START_REF] Derbyshire | Natural minerogenic dust and human health[END_REF].

The entrainment of dust particles from the surface is driven by sandblasting (e.g., [START_REF] Shao | Effect of saltation bombardment on the entrainment of dust by wind[END_REF][START_REF] Marticorena | Modeling the atmospheric dust cycle: 1. design of a soil-derived dust emission scheme[END_REF][START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF][START_REF] Shao | A simple expression for wind erosion threshold friction velocity[END_REF][START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distributions in source areas[END_REF][START_REF] Shao | A model for mineral dust emission[END_REF]Kok, 2011b), and, to a lesser extent, by aerodynamic forces (e.g.,

Preprint submitted to Aeolian Research

February 7, 2019

from an impacting saltator is available to break surface inter-particle cohesive bonds (E). In dust emission schemes, saltation is considered as the necessary process for releasing dust. The dust flux is, therefore, often a function of the saltation flux (e.g., [START_REF] Shao | Effect of saltation bombardment on the entrainment of dust by wind[END_REF][START_REF] Marticorena | Modeling the atmospheric dust cycle: 1. design of a soil-derived dust emission scheme[END_REF][START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distributions in source areas[END_REF][START_REF] Shao | A model for mineral dust emission[END_REF]Kok, 2011b;[START_REF] Kok | An improved dust emission model -part 1: Model description and comparison against measurements[END_REF]. While the link between saltation and dust fluxes is well established, the processes shaping the PSD of the dust flux remain unclear.

In dust emission schemes, the emitted dust PSD is usually influenced by (1) the surface interparticle cohesive bond, (2) the soil PSD, and/or (3) the wind intensity, as reviewed below.

The surface inter-particle cohesive bond (E). The parametrization of E differs from one scheme to the other. In [START_REF] Shao | A model for predicting aeolian sand drift and dust entrainment on scales from paddock to region[END_REF], E was assumed proportional to the drag force acting on the particle multiplied by some length scale, leading to E ∝ d 3 p , where d p is the diameter of the released dust. Later, [START_REF] Shao | A model for mineral dust emission[END_REF] and [START_REF] Kok | An improved dust emission model -part 1: Model description and comparison against measurements[END_REF] related E to the van der Waals bond at the surface, leading to E ∝ d 2 p . Then [START_REF] Shao | Physics and Modelling of Wind Erosion[END_REF] suggested an ensemble surface cohesive bond proportional to the particle size, leading to E ∝ d p ::: The ::::: same ::: was ::::::::: considered ::: by :::::::::::::: [START_REF] Kok | An improved dust emission model -part 1: Model description and comparison against measurements[END_REF] .

Differently, [START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distributions in source areas[END_REF] proposed an inverse relationship between E and the particle size, E ∝ d -0.018 p (deduced from their Table 1). This allowed them to reproduce [START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF] wind-tunnel experiment where finer particles were released under higher wind conditions, i.e., higher energy of the impacting saltators.

The soil particle size distribution. The influence of the soil PSD on the dust flux changes according to the dust emission schemes. [START_REF] Marticorena | Modeling the atmospheric dust cycle: 1. design of a soil-derived dust emission scheme[END_REF] related the dust flux PSD to the percentage of clay present in the soil, i.e., percentage of particles lower than 3.9 µm.

In [START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distributions in source areas[END_REF], the emission of dust is derived from a distribution with three lognormal modes, where the contribution of each mode depends on the strength of the inter-particle cohesive bond and on the wind intensity. In [START_REF] Shao | A model for mineral dust emission[END_REF], the emission of dust is derived from a distribution with two modes ::: two ::::: PSDs, minimally and fully disaggregated modes :::::::::: distributions, with an adjustable weight parameter to determine their relative contribution. Kok (2011b) suggested that aggregates behave like brittle materials where the emitted dust PSD is determined by patterns in which cracks nucleate and propagate. A compilation of published flux measurements in [START_REF] Kok | The physics of wind-blown sand and dust[END_REF] indicated a limited effect of soil texture on the observed dust flux PSD.

This insensitivity of the dust flux PSD to the soil granulometry was also suggested in [START_REF] Reid | Dynamics of southwest asian dust particle size characteristics with implications for global dust research[END_REF].

The wind intensity. The influence of the wind intensity on the emitted dust PSD remains unclear. On one hand, [START_REF] Shao | A model for mineral dust emission[END_REF] and [START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distributions in source areas[END_REF] considered an enrichment of the dust flux in smaller particles with increasing wind intensity, due to the stronger disintegration of aggregates :: by :::::::: adjusting :::: their :::::: weight :::::::::: parameters :: as :: a ::::::: function ::: of :::::::::: atmospheric ::::::: forcing. On the other hand, Kok (2011b) proposed the independence of the emitted dust PSD to the wind intensity, deduced from a compilation of field measurements showing negligible differences of the dust flux PSD under various wind conditions (Kok, 2011a).

The definition of the dust flux simulated by dust emission schemes remained also vague in regard to the measured fluxes used to evaluate-fit them (e.g., [START_REF] Shao | Effect of saltation bombardment on the entrainment of dust by wind[END_REF][START_REF] Marticorena | Modeling the atmospheric dust cycle: 1. design of a soil-derived dust emission scheme[END_REF][START_REF] Shao | A model for mineral dust emission[END_REF][START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distributions in source areas[END_REF][START_REF] Kok | An improved dust emission model -part 1: Model description and comparison against measurements[END_REF]. The dust fluxes simulated by these schemes are most likely the :: the ::::::: surface emitted dust fluxes at the surface since these schemes do not explicitly simulate turbulence nor consider the effect of deposition during emission. However, these simulated fluxes have most often been evaluated-fitted against field measurements performed well above the surface (several meters) and at various distances from the upwind border of the dust source area, i.e. fetch sizes ranging from 100 m to more than 10 km (e.g., [START_REF] Gillette | Fine particulate emissions due to wind erosion[END_REF][START_REF] Gomes | Validation of a dust production model from measurements performed in semi-arid agricultural areas of spain and niger[END_REF][START_REF] Gillies | Eolian suspension above the saltation layer, the concentration profile[END_REF][START_REF] Zobeck | Wind-induced dust generation and transport mechanics on a bare agricultural field[END_REF][START_REF] Fratini | Size-segregated fluxes of mineral dust from a desert area of northern china by eddy covariance[END_REF][START_REF] Sow | Size resolved dust emission fluxes measured in niger during 3 dust storms of the amma experiment[END_REF]. This implies an hypothesis of PSD similarity between near-surface and surface fluxes, and thus a negligible impact of the In this study, the near-surface dust flux will refer to the vertical turbulent-diffusive flux (F wc ) at several meters above the surface. It corresponds to the amount of particles transported away from the surface by the flow turbulence. This flux results from a balance between the emission (F emi ) and deposition (F dep ) fluxes, the storing rate of particles in the air, and the gravitational settling flux (F sed ) due to particle weight. Importantly, F wc does not include F sed and it is not equivalent to the dust emission flux at the surface (F emi ). See Table 1 for a summary of the flux notations used in this paper.

One-dimensional dust-dispersal model

A simple 1D dust-dispersal model was developed to simulate the whole dust life-cycle in a column of air representative of the atmospheric surface boundary layer above an infinite bare erodible surface under neutral thermal stratification (Fig. 1). The model accounts for (1) dust emission through sandblasting, based on an energy budget linking the surface cohesive forces to the dust particle size, (2) dust turbulent transport within the surface boundary layer, and (3) dust deposition at the surface.

Wind

The wind intensity is quantified through the friction velocity (u * ), which is related to the turbulent flow shear stress (τ) induced by the surface:

u * = τ/ρ a , (1)
where ρ a is the air density. The shear stress quantifies the momentum flux absorbed by the surface. Above an homogeneous surface, τ is constant with height within the surface boundary layer. The similarity theory leads to the well-known logarithmic velocity profile expressed as:

u (z) = u * κ ln z z 0 , ( 2
)
where z is the vertical coordinate, z 0 the surface roughness length, and κ the von Karman constant (= 0.40).

In the presence of saltation, the total shear stress (τ) is partitioned within the saltation layer between the momentum flux aborbed by the surface and that absorbed by the saltators. The absorption of momentum by the saltators modifies the wind profile. This was accounted for in [START_REF] Raupach | Saltation layers, vegetation canopies and roughness lengths[END_REF] as:

u (z) = u * κ ln z z 0 + 1 - √ r γ + ln z 0 H s + E z H s , (3)
where H s = 0.3969u [START_REF] Raupach | Saltation layers, vegetation canopies and roughness lengths[END_REF].

Unlike saltating particles, the relatively small volume fraction of dust particles has negligible effect on the flow field.

Dust transport

An Eulerian approach was used to describe the turbulent transport of dust particles in the atmospheric surface layer. Here, dust particles were assumed spherical, non magnetic, and electrically neutral. The size range of dust particles is divided into n b bins, each bin being characterized by a mean particle diameter d b . Hence, the conservation equation of the dust concentration c b of the bth bin is given by:

∂c b (z, t) ∂t = - ∂F wc,b (z, t) ∂z + v s,b ∂c b (z, t) ∂z , (4)
F wc,b (z, t) = -K d ∂c b (z, t) ∂z , (5)
where

K d = κzu * 1 -1 - √ r exp (-z/H s )
is the dust eddy diffusivity assuming similarity in turbulent transport between momentum and dust particles [START_REF] Gillette | Measurements of aerosol size distributions and vertical fluxes of aerosols on land subject to wind erosion[END_REF]. This expression of K d converges to its surface layer value κzu * above the saltation layer [START_REF] Raupach | Saltation layers, vegetation canopies and roughness lengths[END_REF].

At the surface, the total dust flux (F tot ) accounting for both turbulent-diffusive and gravitational settling fluxes, is expressed as:

F sur f tot,b (t) = F emi,b (t) -F dep,b (t) , (6)
where F emi,b and F dep,b are the emission and deposition fluxes at the surface of dust particles of the bth bin, respectively (Fig. 1).

Saltation

Saltating particles are characterized by a size distribution p s whose diameters range from

D S 1 to D S 2 .
The total saltation flux is the integration of the saltation flux of all particle sizes according to p s :

Q tot = D S 2 D S 1 Q(D)p s (D)δD, (7)
where the saltation flux Q(D) of particles of diameter D is estimated from [START_REF] White | Soil transport by winds on mars[END_REF] as follows:

Q(D) = 2.61 ρ a g 1 + u * ,t (D) u *        1 - u * ,t (D) u * 2        u 3 * , (8)
where u * ,t (D) is the threshold friction velocity deduced from [START_REF] Marticorena | Modeling the atmospheric dust cycle: 1. design of a soil-derived dust emission scheme[END_REF].

On average, the kinetic energy of the saltating particles of diameter D impacting the surface is:

E imp,D = 1 2 m imp v 2 imp , (9)
where m imp = ρ p πD 3 /6 is the particle mass, and v imp is the velocity of the impacting saltator. In a first order approximation, v imp is estimated as v imp = 5u * ,t (D) [START_REF] Kok | An improved dust emission model -part 1: Model description and comparison against measurements[END_REF].

Using a probabilistic approach as often used in splash schemes, the kinetic energy of impacting saltators (E imp ) is distributed between the energy used for saltator rebound (fraction reb ), the energy used for ejecting new saltators (fraction e j ), and the energy lost to the surface (fraction f r ), such as:

reb + e j + f r = 1. ( 10
)
Here, reb = 2γ 2 P reb , where P reb is the rebound probability [START_REF] Anderson | Sediment transport by wind: Toward a general model[END_REF]) and γ = 0.55 [START_REF] Rice | An experimental study of multiple grain-size ejecta produced by collisions of saltating grains with a flat bed[END_REF]; and f r = 0.96 (1reb ) [START_REF] Ammi | Three-dimensional analysis of the collision process of a bead on a granular packing[END_REF].

This energy distribution at the surface during the saltator impact does not account for the energy fraction used for releasing dust particles from sandblasting. We hypothesized that the energy for dust emission is a fraction d of the energy lost to the surface ( f r E imp,D ).

Dust emission

At the surface, dust particles lie on larger sand grains, surrounded by other dust grains, forming aggregates. A particle in this system experiences short-range dielectric attractive forces known as van der Waals forces due to inter-particle interactions (dust-dust particles or dustsand particles). Estimating the van der Waals forces between two spherical particles in contact is simple but it becomes much more complex when it comes to estimate the cohesive forces of particles in an aggregate [START_REF] Shao | A model for mineral dust emission[END_REF]. This explains the various parametrizations proposed in the literature for the inter-particle cohesive forces as mentioned in the introduction. Overall, these parametrizations express the inter-particle cohesive bond as a function of the dust particle diameter:

E coh,b = Ad β b , (11)
where β can be referred to as the inter-particle cohesive bond exponent defining the relationship between the dust particle diameter and the cohesive bond, and A is a constant. In existing dust emission schemes, β varies from -0.018 to +3 as reviewed in the introduction section.

The number of dust particles of bin b released by an impacting saltator of size D is defined as the ratio between the available energy for releasing dust and the energy required to eject one dust particle:

N emi,b = d f r E imp,D Ad β b , (12)
where d is the fraction of the energy lost to the surface by the impacting saltator ( f r E imp,D ) that is used to release dust. Equation 12 can be rewritten as:

N emi,b = α f r E imp,D d β b , ( 13
)
where α is the dust emission coefficient resuming the two unknown constants d and A.

Hence, the emission dust flux at the surface of the bth bin is:

F emi,b = ζ d (d b , u * )p d (d b ) D S 2 D S 1 N emi,b (D) Q(D) l(D) p s (D)δD, ( 14
)
where p d (d b ) is the proportion of dust particles from the bth bin available at the surface, Q/l represents the vertical flux of saltating particles or in other words the flux of impacting saltators at the surface, with l the average hop length of the saltating particles. Here, l was simply deduced as the maximum horizontal distance traveled by a projectile launched with a lift-off velocity V e j and a lift-off angle θ e j : l = V 2 e j sin 2θ e j /g. We chose V e j = 0.63u * and θ e j = 50

Dust Deposition

Dust deposition on soil occurs only through dry deposition accounting for gravitational settling, turbulent mixing, and brownian diffusion (e.g., [START_REF] Seinfeld | Atmospheric chemistry and physics: from air pollution to climate change[END_REF]. Hence, the deposition flux of dust particles from the bth bin is:

F dep,b = v dep,b c sur f b , (15)
where c sur f b is the dust concentration close to the surface (middle of the first grid cell, see Fig. 1), and v dep,b is the particle dry deposition velocity that is classically parameterized as a set of resistances:

v dep,b = 1 R a + R b + R a R b v s,b + v s,b , (16)
where R a = log (z s /z 0s ) / (κu * ) is the aerodynamic resistance accounting for turbulent transfer near the surface (z s is the middle of the first grid cell) with z 0s the saltation-layer roughness length estimated as per [START_REF] Raupach | Saltation layers, vegetation canopies and roughness lengths[END_REF]; Fig. 2 presents the variation of the deposition velocity with the dust particle size for different wind intensities. The deposition velocity exhibits a minimum for particles around 0.5 to 1 µm in diameter. With increasing wind intensity, the deposition velocity increases, especially for the largest particles, due to the larger decrease of R b with increasing particle size.

R b = u * S -2/3 c + 10 -3/S t -

Simulation configurations

The computational domain was meshed along the vertical using a stretched grid from the surface to the top of the domain. The grid size was varying from ∆z min = 0.01 m at the surface to ∆z max = 0.24 m at the top (Fig. 1). This mesh ensures a grid resolution fine enough to simulate exchanges close to the surface (where particle concentration is high), while reducing the computational time by considering larger grid cells for exchanges farther away. Unless otherwise specified, the total domain height was 200 m corresponding to the approximate depth of the atmospheric surface layer under neutral thermal stratification.

The dust conservation equation (Eq. 4) was solved using the Crank-Nicolson method [START_REF] Crank | A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type[END_REF]:

c t+∆t b -c t b ∆t = 1 2 F t+∆t i c b , z, t + ∆t, ∂c b ∂z , ∂ 2 c b ∂z 2 + F t i c b , z, t, ∂c b ∂z , ∂ 2 c b ∂z 2 , ( 17
)
where ∆t = 0.01 s is the time step and F t i is the right-hand side term of Eq. 4 at time t. A Neuman boundary condition was applied at the top of the domain (∂c b /∂z = 0). This resolution method has the advantage of being implicit and of ensuring second order convergence in time. The spatial derivatives were discretized following a finite volume approach, ensuring the conservation of the amount of dust particles within the computational domain. In our simulation analysis, the time variations of dust concentration and dust flux profiles were converted into spatial variations from an upwind virtual border of an erodible surface, corresponding to t = 0 s, to distances x = U int t, referring to the fetch length, where U int is the integral wind speed (U int = z re f 0 u (z) dz/z re f , where z re f was chosen equal to the height of the investigated dust flux, 3 m here, unless otherwise specified). This time-space conversion allows us to evaluate the sensitivity of the dust flux PSD to the deposition process according to the fetch length.

Model evaluation

Before using our model to investigate the sensitivity of the near-surface dust flux PSD, we find it important to evaluate its ability to reproduce published dust flux behaviors.

Dust concentration vertical profile compared to analytical solution

At steady state over an extended homogenous surface (F emi,b = F dep,b ), the conservation equation of dust concentration (Eq. 4) leads to the following analytical solution [START_REF] Shao | Physics and Modelling of Wind Erosion[END_REF]:

c b (z) = c b (z r ) (z/z r ) v s,b /(κu * ) , (18)
where z r is a reference height.

The concentration profile simulated by our model at equilibrium between F emi,b and F dep,b is consistent with this analytical solution (Fig. 5). This was verified for an erosion event with u * = 0.40 ms -1 and 16.0 µm dust particles, in a 35 m high domain. This equilibrium was reached for a fetch x eq = 36 km. With a larger domain, x eq increases as the emitted dust particles have a larger volume to disperse, reducing the near-surface dust concentration, and thus the deposition.

For typical surface atmospheric boundary layer of 200 m high, x eq ≈ 200 km. The equilibrium fetch also increases with decreasing particle size. For 5 and 10 µm dust particles and a 200 m high domain, x eq ≈ 320 and 640 km, respectively. Simulations with particles smaller than 5 µm led to an equilibrium fetch larger than a few thousand kilometers. within a large-eddy simulation (LES) airflow model that simulated instantaneous wind. They simulated 20-minute long erosion events, with three dust particle bins 1.5, 6.7, 14.2 µm, under three wind conditions (u * = 0.47, 0.63, 0.77 ms -1 ). By assuming an equiprobable emission of the three size bins, they observed an enrichment of the 2-m high F wc in small particles with time.

This enrichment was explained by the nonstationarity of the erosion process, due to the lower deposition velocity of the smaller dust particles. Here, we reproduced these simulations with our simple 1D model using similar erosion configurations as in [START_REF] Dupont | Near-surface dust flux enrichment in small particles during erosion events[END_REF]. µm) and the impoverishment in 6.7 µm particles with time and with increasing wind intensities.

The small differences between both models may be explained by the difference in modelling the turbulence in the saltation layer between both approaches. (1997).

Here, our 1D model is used to reproduce the [START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF] experiment by simulating erosion events in a domain with the same vertical size (70 cm) as their wind tunnel. The interparticle cohesive bond exponent β was set to -0.018 and the surface dust size distribution (Fig.

3

) was chosen so as the F emi PSD for u * = 0.40 ms -1 equates that observed by [START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF]. Two simulations were conducted: one where F emi PSD exhibited no dependence to the wind intensity (ζ(d b , u * ) = 1 in Equation 14), and the other one where F emi PSD followed the dependence to the wind intensity proposed by [START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF] 3

). The 3-m high F wc PSD in number simulated by our model (Fig. 8a) is close to that of Kok (2011b) for small fetch lengths (< 5 m). As the fetch increases, the fraction of the largest particles in the dust flux decreases due to their surface deposition, and thus the flux enriches in small particles. This enrichment in small particles (impoverishment in large particles) is accentuated with u * as the deposition velocity increases (Fig. 2). The differences between the simulated nearsurface dust flux PSD and the emitted flux PSD proposed by Kok (2011b) are amplified when the PSDs are expressed in mass as it emphasizes the role of the largest particles (Fig. 8b).

Sensitivity of the dust flux PSD

To investigate the sensitivity of F wc PSD, 10 simulations were conducted varying either (1) the inter-particle cohesive bond exponent (β), ( 2) the surface dust-size distribution (p d ), or (3) the wind intensity (u * ). Each simulation started from an air clean of dust. This sensitivity of F wc PSD is presented according to the fetch size, with values ranging from 5 m to 10 km. To focus solely on the possible impact of the deposition process on F wc PSD, the F emi PSDs were considered independent of the wind intensity, i.e. ζ(d b , u * ) = 1 in Equation 14.

Table 3 summarizes the range of each parameter-forcing considered in this analysis. For p d , four surface dust-size distributions (D1, D2, D3 and D4) were considered (Fig. 3). They were chosen so as to obtain for β = 2 the following F emi PSDs (Fig. 9d): equinumber emissions across bins (D1), strong emission around 0.7 µm corresponding to particles with the lowest deposition velocity (D2), two emission peaks at 0.7 and 8 µm (D3), and strong emission around 8 µm corresponding to particles with large deposition velocity (D4). Fig. 9 presents the F emi PSDs obtained for the different combinations of β and p d considered in our analysis. These PSDs cover the usual distributions encountered in existing emission schemes [START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distributions in source areas[END_REF][START_REF] Shao | A model for mineral dust emission[END_REF][START_REF] Kok | An improved dust emission model -part 1: Model description and comparison against measurements[END_REF].

The variation of the 3-m high F wc PSDs is presented in number in presents the fraction in number and in mass of F wc partitioned in four bins (Z1 to Z4) on the total dust flux F wc,tot , according to the fetch size. Here, these four bins cover the whole dust size range.

They were defined based on the variations of the particle deposition velocity with the particle size (Fig. 2) such as Z1 bin (< 0.3 µm) covers the region of predominant brownian deposition, Z2 bin (0.3 -1.0 µm) covers the region of minimum deposition velocity, Z3 bin (1.0 -8.0 µm) covers the region of sharp rise in deposition velocity, and Z4 bin (> 8.0 µm) covers the region of predominant gravitational settling.

An enrichment of F wc in particles with the lowest deposition velocity (Z2 bin) is observed with the fetch when the emitted dust at the surface (F emi ) is deprived in dust from the Z2 bin (Fig. 10a; D3, D4 in Fig. 10b; β = 0, 1 in Fig. 10c). This enrichment slightly extents to particles of the Z1 bin (Fig. 11a; D3, D4 in Fig. 11b; β = 0, 1 in Fig. 11c), and implies an impoverishment in particles from Z3 and Z4 bins. When F emi is dominated by particles of the Z2 bin, no enrichment (impoverishment) is observed as F emi is already dominated by particles the least likely to deposit (D2 in Figs. 10b, 11b and β = 3 in Figs. 10c,11c).

The resulting difference from this enrichment between F wc and F emi PSDs represents a few percent of the flux for fetch lengths around 100 m, but can rise to about 10% for fetch longer than 1 km (Fig. 11). These percentages are observed for fluxes expressed either in number or in mass.

This modification of F wc PSD with the fetch occurs mainly within the first few hundred meters from the upwind border of the source area, then it evolves slowly with the fetch (Fig. 10a). This is visible, for example, in the simulation with the D1 surface PSD and u * = 0.50 ms -1 (Fig. 11a)

where the fraction of F wc from bins Z3 and Z4 decreases by nearly 5% within the first 200 m, followed by only a 2% decrease up to 10 km.

The enrichment of F wc in particles of the Z2 bin is accentuated with the wind speed (Figs.

10a and 11a). For a D1 surface PSD and β = 2, the fraction of F wc from bin Z2 increases by about 2% at a fetch around 100 m between u * = 0.30 and 0.50 ms -1 , and by about 5% at fetch longer than 1 km (Fig. 11a). In a different way, for u * = 0.50 ms -1 , F wc from bin Z2 already reaches in number 39% of the total F wc for a fetch of 100 m while a fetch of 10 km is needed to reach the same fraction for u * = 0.30 ms -1 (Fig. 11a). Furthermore, the particle size below which such enrichment occurs decreases from 5 µm to 3 µm (Fig. 10a) when u * increases from 0.30 to 0.50 ms -1 , emphasizing the role of the wind speed in shaping F wc PSD through its impact on the deposition velocity.

The strength of F wc enrichment in small particles depends on the F emi PSD. This latter one is shaped by the combination of the surface PSD and the inter-particle cohesive bond exponent (β) (Fig. 9), β controlling the difference in PSD between the dust available at the surface and the dust emitted from the surface. Regardless of the surface PSD, a high positive value of β (≈ 3) results in a strong emission of the smallest particles, and conversely a strong emission of the largest particles occurs for a small value of β (≈ -1) (Figs. 10c and11c). The enrichment of F wc in dust of the Z2 bin is minimal for these extreme cases due to either the insignificant concentration of small particles (β ≈ -1) or the predominant concentration of particles of the Z1-Z2 bins (β ≈ 3) (Fig. 10c).

Discussion and conclusion

The particle size distribution (PSD) of the near-surface dust flux (F wc ) was found to be mainly influenced by both the deposition process and the size distribution of the emitted dust (F emi ).

The F emi PSD appeared shaped from the mutual choice of the size distribution of available dust at the surface (p d ) and the inter-particle cohesive bond exponent (β) (Fig. 9). Different combinations of both quantities can lead to the same F emi PSD. This explains why [START_REF] Shao | A model for mineral dust emission[END_REF] and [START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distributions in source areas[END_REF] were able to reproduce from their schemes the F wc PSD observed by [START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF] while using different values of β, 2 and -0.018, respectively, and different surface dust size distributions (Fig. 3). Hence, the uncertainty of both quantities in existing dust emission schemes can be simply resolved by defining directly the F emi PSD. Defining F emi PSD becomes, however, all the more complex if p d is wind speed dependent as suggested by [START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF].

For constant wind and dust emission conditions, our simulations suggested a distinction between the PSDs of F emi and F wc . This PSD difference increases with the length of the fetch, i.e.

the distance from the upwind border of the source area. A difference of a few percent in the flux fraction in number of small dust particles (0.3 to 1 µm) between F emi and F wc PSDs is observed for 100-m long fetch (Fig. 11). This difference can rise to 10% for fetch longer than 1 km. In mass, this difference affects rather larger particles (> 1 µm). This modification of F emi and F wc PSDs is explained by the deposition process, as previously found by [START_REF] Dupont | Near-surface dust flux enrichment in small particles during erosion events[END_REF] from LES. The deposition of dust particles slowly increases with the fetch as the concentration of dust in the air is enhanced. This process is particle size dependent due to the dependence of the deposition velocity on the particle size (Fig. 2).

Theoretically, the dust deposition should increase with the fetch until an equilibrium is reached between deposition and emission (stationary state). Our simulations showed that this steady state occurs for fetch lengths larger than a few thousand kilometers for particles lower than 5 µm (section 3.2). This suggests a low probability of reaching a dust steady state during an erosion event. However, this result was obtained for an air clean of dust at the beginning of the simulation, i.e. a clean air upwind from the source area. For a background air already charged in dust particles passing through a source area, the deposition is expected to be larger and the equilibrium between deposition and emission should be reached for smaller fetch lengths. In that configuration of a background dust concentration, the magnitude of the difference in PSD between F emi and F wc should be larger than that observed from our simulations.

Emission and deposition appeared in our simulations as parallel inseparable processes, the impact of the deposition on the F wc PSD being already visible (a few percent) for short fetch around 100 m (Fig. 11). These parallel emission and deposition processes may question the evaluation of existing dust emission schemes. As mentioned in the introduction, these schemes simulate the PSD of emitted dust at the surface (F emi ) but have been usually evaluated against measurements performed well above the surface (F wc ) where deposition should probably have started to sort particles. Since dust emission schemes do not account for deposition, it was assumed that the PSD of F emi is identical to that of F wc , and time-space independent for stationary wind conditions. Our results suggest that this evaluation procedure is not appropriate when comparing with field experiments, for which fetch lengths usually range from 100 m to several kilometers. However, the small fetch length of wind-tunnel experiments may ensure a PSD similarity between near-surface and surface emitted dust fluxes, as shown in our simulation of the wind-tunnel experiment of [START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF] (section 3.4).

The wind intensity modifies the F wc PSD by increasing the surface deposition, especially for large particles (Fig. 2), and thereby accelerating the attainment of equilibrium between emission and deposition. A rise of the friction velocity (u * ) from 0.30 to 0.50 ms -1 leads to a 10% increase (decrease) of the dust flux fraction of small (large) particles for 1-km long fetch (Fig. 11a).

Hence, an observed enrichment of F wc in small particles with wind intensity could be explained by the enhancement of the deposition of larger particles for fetch lengths larger than 100 m, adding to the higher release of small dust through a stronger disintegration of aggregates as proposed by [START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF].

The magnitude of the difference between F emi and F wc PSDs depends also on the difference in diameter between the particles the most emitted at the surface and the particles with the lowest deposition velocity (around 0.7-2.0 µm following Fig. 2). The farther the peak of F emi PSD is from the minima of the deposition velocity, the greater is the difference in PSDs between F wc and F emi (Fig. 10b,c). Again, this is explained by the dependence of the deposition velocity to the particle size. The near-surface dust flux enriches in particles with the lowest deposition velocity.

For F emi dominated by particles with the lowest deposition velocity, the PSDs between F emi and F wc remain similar regardless of the fetch (D2 in Fig. 10b, β=3 in Fig. 10c). This similarity in PSD does not mean an equilibrium between the emission and deposition processes (steady state).

In conclusion, the role of the deposition process when evaluating existing dust emission schemes against field experiments should be considered, especially for large fetch lengths and/or with a large-scale background dust concentration. Improving our prediction of the PSD of F wc requires not only a better understanding of the size distribution of F emi but also a better understanding of the deposition velocity according to the particle size. This is all the more complex as both, emission and deposition processes, are difficult to disentangle in experiments.

Flux

Parameter

Range of values Friction velocity (u * ) 0.30, 0.40, 0.50 ms -1 PSD of available dust at the surface (p d ) D1, D2, D3, D4 (see, Fig. 3) Inter-particle cohesive bond exponent (β) -1, 0, 1, 2, 3

  the :::::::::: gravitational ::::::: settling :: of ::: the :::::: largest :::::::: particles ::: and ::: the ::::::: particle :::::: surface ::::::::: deposition :::::: across ::all ::::: sizes. ::::: Only :::::::::::::::::::::[START_REF] Shao | Parameterization of size-resolved dust emission and validation with measurements[END_REF] applied :: a ::::::::: correction :: to ::: the :::::::: measured :::: dust :::: flux :: to ::::::: account ::: for ::the :::::::::: graviational ::::::: settling ::: but :::::: without :::::::::: considering ::: the ::: full : deposition process.[START_REF] Dupont | Near-surface dust flux enrichment in small particles during erosion events[END_REF] obtained from large-eddy simulations (LES) of aeolian erosion events an increasing difference in PSD with time between near-surface dust flux and emission flux as a result of a first particle sorting through turbulence diffusion, gravitational settling, and more importantly surface deposition. Their results were, however, limited to three particle sizes (1.40, 6.70 and 14.20 µm) and to one soil PSD with an equiprobable emission of the three dust particle sizes, limiting a more general view of the impact of dust deposition on the near-surface dust flux PSD according to the soil PSD. Interestingly, by converting the erosion duration of their simulations into a fetch size, their results mean that the difference in PSD between near-surface dust flux and emission flux should increase with the fetch. A quantification of this difference according to the fetch size and the soil PSD would benefit to the erosion community.The goal of the present study is twofold: (1) investigate the role of deposition in shaping the PSD of near-surface dust flux during erosion events, and, more generally, (2) investigate the sensitivity of the PSD of near-surface dust flux to the surface inter-particle cohesive bond, the PSD of available dust at the surface, the wind intensity, and the fetch length. To that purpose we extend the initial study of[START_REF] Dupont | Near-surface dust flux enrichment in small particles during erosion events[END_REF] by considering (1) dust particles ranging from 0.1 to 16 µm, (2) various soil PSDs, and (3) the influence of the inter-particle cohesive bond exponent (β) on the emitted dust flux PSD. Because the LES approach used in Dupont et al. (2015) was too time consuming for performing such a detailed sensitivity analysis, we developed for this analysis a simple one-dimensional (1D) physically-based dust dispersal model, accounting for dust emission, transport, and deposition in the atmospheric surface boundary layer. After evaluating the model against published dust flux behaviors, the model is used to study the sensitivity of the PSD of the near-surface dust flux.

  where t is time, F wc,b is the dust turbulent-diffusive flux of the bth size bin, v s,b is the settling velocity of the particles of the bth size bin given as v s,b = ρ p gd 2 b C c / (18ρ a ν), where ρ p is the dust particle density, C c is the Cunningham slip correction factor, and ν is the kinematic viscosity of air (Seinfeld and Pandis, 1998). Equation 4 means that the dust concentration variation in a grid cell results from a balance between incoming and outgoing turbulent-diffusive (F wc,b ) and gravitational settling (F sed,b = v s,b c b ) fluxes due to turbulent motions and particle weight, respectively.

  [START_REF] Shao | Physics and Modelling of Wind Erosion[END_REF]. The coefficient ζ d (d b , u * ) accounts for the modification of the size distribution of emitted dust due to the enhancement of aggregate disintegration with wind speed as proposed by[START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF]. Unless otherwise specified, ζ d (d b , u * ) = 1 in our simulations, implying the size distribution of the emission dust flux is independent of the wind intensity.

  1 is the quasi-laminar resistance accounting for brownian diffusion and inertial impaction on particles. Here, S c = ν/D g is the Schmidt number, D g = k b TC c / 3πρ p νd b is the brownian diffusivity, with k b the Boltzman constant, and T the air temperature, set to 27 • C in this study. S t = u 2 * v s,b / (gν) is the Stokes number.

3. 1 .

 1 Dust flux versus wind intensity compared to existing dust emission schemesExisting dust emission schemes show a clear increase of the dust flux with increasing wind intensity. Here, our model was used to simulate 15-minute erosion events (fetch up to 20 km) with different stationary wind intensities (u * ), starting with an air clean of dust. The surface dust size distribution (p d ) was considered enriched in small particles (D1 distribution in Fig.3). The cohesive bond exponent β was set to 2. This combination of β = 2 and distribution D1 leads to an equiprobable PSD of emitted dust. The dust emission coefficient (α) that modulates the amplitude of the dust flux in our model (Eq. 13), was chosen so as the simulated dust flux at u * = 0.40 ms -1 fits the one fromShao et al.'s (1993) scheme.The 3-m high F wc,tot , including all particle sizes, as simulated by our model after 15-minute erosion, exhibits a similar trend with u * as the fluxes obtained from the emission schemes of[START_REF] Shao | Effect of saltation bombardment on the entrainment of dust by wind[END_REF] and[START_REF] Marticorena | Modeling the atmospheric dust cycle: 1. design of a soil-derived dust emission scheme[END_REF] (Fig.4). This agreement confirms the accurate sensitivity to the wind intensity of the magnitude of the dust flux simulated by our model. The same was verified for different combinations of β and p d by fitting different values of α.

3. 3 .

 3 Dust flux enrichment in small particles compared to Dupont et al. (2015) large-eddy simulation Dupont et al. (2015) performed detailed three-dimensional simulations of soil erosion of a bare surface by representing the main erosion processes (saltation, sandblasting, dust suspension)

Fig. 6

 6 Fig. 6 presents the time variation of the fractions in number of F wc,b on F wc,tot , for the three size bins and three wind intensities, as estimated from our model and compared to Dupont et al. (2015). Our 1D model reproduces quite efficiently the flux enrichment in the smallest bin (1.5

3. 4 .

 4 scheme, equivalent to β = -0.018. Later,[START_REF] Shao | A model for mineral dust emission[END_REF] validated his dust emission scheme by reproducing similar time-independent PSD of the dust flux as a function of the wind intensity. To obtain his result, he chose β = 2 with a different surface dust size distribution than Alfaro et al.

  Fig. 7 compares the mass size distributions of F wc simulated by the current model and observed by Alfaro et al. (1997), for a 5-m long fetch and three friction velocities. For u * = 0.40 ms -1 , our model agrees with the observations of Alfaro et al. (1997) regardless of ζ(d b ). At higher wind speeds, our results only agree with Alfaro et al.'s ones when F emi PSD is a function of the wind speed (ζ(d b , u * ) 1). This result implies a negligible effect of deposition, sedimentation and turbulence, in shaping the PSD of the near-surface dust flux and, thus, in differentiating the PSDs of F wc and F emi , for the conditions of the wind tunnel experiment. This result leaves open the possibility of the dust flux enrichment in small dust particles with wind intensity due to an enhancement of aggregate disintegration[START_REF] Alfaro | Modeling the size distribution of a soil aerosol produced by sandblasting[END_REF].

  Fig. 10 according to the fetch size. Fig. 10 includes results for (a) three wind conditions (u * = 0.30, 0.40 and 0.50 ms -1 ) with β = 2 and the surface dust-size distribution D1, (b) three surface dust-size distributions (D2, D3 and D4) with β = 2 and u * = 0.40 ms -1 , and (c) four surface cohesive bond exponents (β = -1,0,1 and 3) with the surface dust-size distribution D1 and u * = 0.40 ms -1 . Additionally, Fig. 11

Figure 1 :Figure 2 :Figure 3 :Figure 5 :Figure 6 :

 12356 Figure 1: Schematic representation of the 1D computational domain used in this study.

Figure 7 :Figure 8 :Figure 9 :Figure 10 :

 78910 Figure 7: The mass size distribution of the vertical dust fluxes (F wc ) as simulated by our model (line-dots) for three wind conditions (u * ), and compared with the wind-tunnel observations from Alfaro et al. (1997) (filled squares). In one simulation, the PSD of available dust at the surface (p d ) is independent of u * (ζ d = 1) while in the other simulation p d depends on u * as suggested by Alfaro et al. (1997) (ζ d 1)

Figure 11 :

 11 Figure 11: Variation with the fetch length of the fraction in number (left figures) and in mass (right figures) of the 3-m high F wc on the total turbulent-diffusive dust flux F wc,tot according to (a) the wind intensity (u * = 0.30, 0.40 and 0.50 ms -1 ), (b) the surface dust-size distributions (D2, D3 and D4, see Fig. 3), and (c) the inter-particle cohesive bond exponents (β = -1,0,1 and 3). F wc is divided into 4 bins: Z1 (d p ≤ 0.3 µm) -particles with dominant brownian deposition, Z2 (0.3 µm < d p ≤ 2 µm) -particles the least likely to deposit, Z3 (2 µm < d p ≤ 8 µm) -particles with the deposition velocity the most sensitive to wind intensity, and Z4 (d p > 8 µm) -particles the most likely to deposit. The fractions at the smallest fetch lengths are equivalent to those of the surface emitted dust flux (F emi ).

Table 2

 2

summarizes the main input parameters required by our model to simulate dust dispersal during an erosion event. Simulations started with an air clean of dust. Saltating particles had a geometric mean diameter of 210 µm and a geometric standard deviation of 1.8. Dust particle diameters ranged from 0.1 to 16 µm, divided into 15 log-normal bins.

Table 1 :

 1 Dust flux notation.

	Parameter

Table 2 :

 2 Main input parameters of our 1D dust dispersal model.

Table 3 :

 3 Range of values of the parameters considered in our sensitivity analysis of the dust flux PSD.

), including all particle sizes, as a function of the friction velocity (u * ) obtained from our 1D model and compared with the empirical parametrizations of (1) [START_REF] Marticorena | Modeling the atmospheric dust cycle: 1. design of a soil-derived dust emission scheme[END_REF]: F wc,tot = γQ tot , with γ = 10 -2 , and (2) [START_REF] Shao | Effect of saltation bombardment on the entrainment of dust by wind[END_REF]:

), with C = 4 and u * t = 0.2 ms -1 . Note that the dust emission coefficient (α) of our model (Eq. 13) was chosen so as F wc,tot fits the value obtained from [START_REF] Shao | Effect of saltation bombardment on the entrainment of dust by wind[END_REF] scheme at u * = 0.40 ms -1


Mostafa Minh Quan Nguyen

Safdari Shadloo

Abdellah Hadjadj

Benoit Lebon

Jorge Peixinho
email: jorge.peixinho@ensam.eu

Perturbation threshold and hysteresis associated with the transition to turbulence in sudden expansion pipe flow

Keywords: Transition to turbulence, Sudden expansion pipe flow, Instability control, Direct numerical simulation (DNS)

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The flow through an axisymmetric sudden expansion in a circular pipe is a basic configuration, which occurs in many industrial applications, such as heat exchanger, mixing chamber, combustion chamber, etc. This basic geometry is also used as a building block to model more complex flows such as those occurring in arterial stenoses [START_REF] Pollard | A contribution on the effects of inlet conditions when modelling stenoses using sudden expansions[END_REF], pistons [START_REF] Boughamoura | Numerical study of a piston-driven laminar flow and heat transfer in a pipe with a sudden expansion[END_REF], and transportation pipes [START_REF] Koronaki | Numerical study of turbulent diesel flow in a pipe with sudden expansion[END_REF], among others. In these applications, the capacity of predicting when the flow will become turbulent is crucial. In the literature, there are many efforts in theoretical analysis [START_REF] Teyssandiert | An analysis of flow through sudden enlargements in pipes[END_REF], experimental explorations [START_REF] Back | Shear-layer flow regimes and wave instabilities and reattachment lengths downstream of an abrupt circular channel expansion[END_REF][START_REF] Latornell | Some observations on the evolution of shear layer instabilities in laminar flow through a sudden expansion[END_REF] and numerical simulations [START_REF] Macagno | Computational and experimental study of a captive annular eddy[END_REF][START_REF] Varghese | Direct numerical simulation of stenotic flows. Part 1. Steady flow[END_REF] focusing on this problem, or a similar geometries such as the planar abrupt expansion [START_REF] Fearn | Nonlinear phenomena in a symmetric sudden expansion[END_REF][START_REF] Xia | Imaging velocity profiles: Flow through an abrupt contraction and expansion[END_REF][START_REF] Baloch | On two-and threedimensional expansion flows[END_REF][START_REF] Bertolotti | Numerical and experimental models of post-operative realistic flows in stenosed coronary bypasses[END_REF][START_REF] Varghese | Direct numerical simulation of stenotic flows. Part 1. Steady flow[END_REF][START_REF] Tsukahara | DNS of viscoelastic turbulent channel flow with rectangular orifice at low Reynolds number[END_REF]. More recently, Lebon et al. (2018a) found a new mechanism for periodic bursting of the recirculation region in the flow of a circular pipe with the expansion ratio of 1:2. Yet, a consensus about the sequence of events in the transition from laminar to turbulence seems to be relatively well-established, but the exact value of critical Reynolds number is still not firmly determined, and the different transition scenarios are not yet fully elucidated.

The flow is mainly controlled by the inlet Reynolds number, Re = U d/ν, where U is the bulk velocity at the inlet, d is the inlet diameter and ν is the fluid kinematic viscosity. In the laminar state, the flow is axisymmetric. As Re increases, the flow starts to break its symmetrical properties but remains steady until Re = 1139 ± 10 as reported in the experiments by [START_REF] Mullin | Bifurcation phenomena in the flow through a sudden expansion in a circular pipe[END_REF]. It should be noted that the onset of symmetry breakdown is at much lower value for the case of channel sudden expansion, as reported to be Re ≈ 40 by [START_REF] Fearn | Nonlinear phenomena in a symmetric sudden expansion[END_REF] and Re = 216 by [START_REF] Drikakis | Bifurcation phenomena in incompressible sudden expansion flows[END_REF]. In the sudden circular pipe expansion flows, oscillatory and intermittent bursts were reported to appear at Re ≈ 1567 ± 16 by [START_REF] Mullin | Bifurcation phenomena in the flow through a sudden expansion in a circular pipe[END_REF], 1500 < Re < 1700 by [START_REF] Sreenivasan | An instability associated with a sudden expansion in a pipe flow[END_REF], and Re = 750 by [START_REF] Latornell | Some observations on the evolution of shear layer instabilities in laminar flow through a sudden expansion[END_REF]. Then, the flow develops into localised turbulence at even higher Reynolds numbers (Lebon et al., 2018b).

Recently, a global stability analysis was performed by [START_REF] Sanmiguel-Rojas | Global mode analysis of a pipe flow through a 1:2 axisymmetric sudden expansion[END_REF]. They showed that the flow can remain axisymmetric up to Re ≈ 3273, which is much larger than the value found experimentally. Subsequent simulations by [START_REF] Cliffe | Adaptivity and a posteriori error control for bifurcation problems II: Incompressible fluid flow in open systems with Z 2 symmetry[END_REF] also indicated that the steady supercritical bifurcation point lies at even higher Reynolds numbers, i.e. Re ≈ 5000. Moreover, a detailed study of transient growth stability was performed by [START_REF] Cantwell | Transient growth analysis of flow through a sudden expansion in a circular pipe[END_REF]. They showed that the sudden expansion amplifies the energy of infinitesimal perturbations up to six orders of magnitude in the inlet and then decay. The difference between the critical Re in experiments and simulations may be explained by the fact that, in experimental studies, the imperfections of the apparatus have a strong impact on the measured critical Reynolds number. Therefore, the values of critical Reynolds number seem to be dependant both on the perturbation nature and its amplitude. In this case, a numerical simulation with a well-defined finite amplitude perturbation is required to better understand the underlying physics.

Many direct numerical simulations (DNS) on 1:2 expansions were performed before, but with an initial arbitrary velocity field. For example, [START_REF] Tsukahara | DNS of viscoelastic turbulent channel flow with rectangular orifice at low Reynolds number[END_REF] showed in a orifice configuration how the turbulent kinetic energy evolves into localised turbulence. [START_REF] Moallemi | Instability and localized turbulence associated with flow through an axisymmetric sudden expansion[END_REF] used a steady parabolic inflow in a circular pipe expansion and showed the emergence of an instantaneous fluctuating vorticity fields with a maximum near the reattachment point of the recirculation region. The first DNS of finite amplitude perturbation in sudden expansion flow was performed by Sanmiguel-Rojas and [START_REF] Sanmiguel-Rojas | Finite-amplitude solutions in the flow through a sudden expansion in a circular pipe[END_REF], where many interesting results are reported. They showed that, for a range of Re, the flow in laminar state can be forced to disordered motion when a transverse velocity perturbation is added at the inlet. The minimal amplitude of the transverse velocity perturbation required to initiate disordered motion scaled with Re -0.006 . This result has an important role in the passive flow control. Now, it is natural to ask if this result is universal? Will another kind of perturbation have the same behaviour? Is the scaling law still valid? The question is crucial especially by knowing that the previous perturbation scheme used by Sanmiguel-Rojas and Mullin ( 2012) is simply a transverse velocity disturbance in the velocity field. Therefore, it does not satisfy the no-slip boundary condition at the inlet wall. Another interesting result is: when Sanmiguel-Rojas and Mullin (2012) increased and then decreased Re, a hysteresis loop was found for 1450 < Re < 1850. This point deserves more attention, since in the original work, the process of variation of Reynolds number as well as the physical time of the reported flow state are not specified. Recently, [START_REF] Selvam | Flow in a circular expansion pipe flow: effect of a vortex perturbation on localised turbulence[END_REF] presented results on a vortex perturbation but no systematic study of the amplitude threshold and the hysteresis was carried out. One may ask, will the hysteresis occur if Re varies in quasi-static manner? Or will the results change if the observation time is different? In the present study, the authors propose a similar study, but with the vortex disturbance in order to revisit the threshold scaling and the hysteresis loop behaviour.

Numerical set-up

The present work focuses on a circular pipe flow with a sudden expansion with the expansion ratio of 1:2. The fluid flow system is solved using Nek5000 [START_REF] Fischer | nek5000 Web page[END_REF], a well-validated high-order spectral element code for transitional and fully turbulent flows [START_REF] Selvam | Localised turbulence in a circular pipe flow with gradual expansion[END_REF][START_REF] Ducoin | Direct numerical simulation of flow instabilities over Savonius style wind turbine blades[END_REF]. The governing equations are mass and momentum conservations in an isothermal incompressible limit:

∇ • u = 0 (1) ∂ ∂t u + u • ∇u = -∇p + 1 Re ∆u ( 2
)
where u is the velocity field, p is the pressure, and t is the physical time.

The density of flow is set to unity for simplicity. The computational domain is axisymmetric as depicted in Fig. 1(b). The region upstream of the expansion is called inlet and has a diameter, d, and a 5d length. The downstream region, after the expansion, has a diameter D = 2d and a length L = 150d. The whole domain contains 62 300 spectral elements where each element consists of P 3 Gauss-Legendre-Lobatto (GLL) points, P being the polynomial order. P = 5, with total number of 7.9 million calculating points, is used for the reported results. Moreover, additional simulations, performed with P = 6 (13.5 million GLL points) for Re = 1000, showed the laminar recirculation length and the drag coefficient remain essentially the same. A classical set up for the inlet velocity, located at z = -5d, is the Hagen-Poiseuille profile, which satisfies the non-slip condition at the inlet wall:

u o inlet (r = d/2, z = -5d) = 0 when r = √ x 2 + y 2 = d/2 and u o inlet (x, y, z = -5d) = 2U 1 -4(x 2 + y 2 )/d 2 e z , (3)
where (e x , e y , e z ) are the three unit vectors of the cartesian base.

In order to initiate turbulence in subcritical Re, a vortex perturbation is added to the inlet parabolic profile [START_REF] Selvam | Flow in a circular expansion pipe flow: effect of a vortex perturbation on localised turbulence[END_REF]. This modifies the expression of velocity inlet as:

u inlet (x, y, z = -5d) = u o inlet + AΩ    -y x 0    =         -AΩy AΩx 2   1 -4(x 2 + y 2 r 2 )            , (4)
where A and Ω are the amplitude and the intensity of the vortex perturbation, respectively. Ω is defined by (i) the vortex radius, R Ω = d/4, (ii) the center (x Ω , y Ω ) and (iii) the distance from the centerline to the center of the vortex:

r Ω = (x -x Ω ) 2 + (y -y Ω ) 2 .
Then, the expression of Ω is given by:

Ω =      1, r Ω ≤ R Ω /2 8 (R Ω -r Ω ) /d, R Ω /2 ≤ r Ω ≤ R Ω 0, r Ω > R Ω (5)
The position of the vortex is fixed with x Ω = -d/4 and y Ω = 0, such as the inlet profile obeys the incompressible, non-slip and non-penetrate boundary conditions. The radius and the position of the vortex, r Ω , is chosen such that Ω vanishes at the boundary. By injecting Ω = 0 into Eq. ( 4), one can easily verify that u inlet vanishes at the boundary as well.

Preliminary tests with a large vortex perturbation on the pipe centreline indicate the required amplitude to initiate turbulence is relatively large. It corresponds to the limit case of the rotating Hagen-Poiseuille flow discharging into a sudden expansion [START_REF] Miranda-Barea | Experimental study of rotating hagen-poiseuille flow discharging into a 1:8 sudden expansion[END_REF]. By positioning the vortex at (x Ω , y Ω ) = (-d/4, 0), the vortex breaks the axial symmetry [START_REF] Wu | Osborne Reynolds pipe flow: Direct simulation from laminar through gradual transition to fully developed turbulence[END_REF] of the flow and naturally deflects the recirculation region. The distribution of Ω(x, y) can be visualised in Fig. 1(a).

Distinction between different states of instabilities

The laminar or turbulent state of a flow can be monitored from the time evolution of the drag coefficient:

C z (t) = ν πdLU 2 L z=0 2π θ=0 ∂u z ∂r r=d r dθ dz, (6)
where d and L = 150d are the radius and the length of the downstream pipe, respectively. Here r = √ x 2 + y 2 and θ = arctan (y/x) are the positions in the radial and azimuthal direction in the cylindrical coordinate system, respectively.

Table 1 presents a summary of the selected simulations along with parameters: Re, A, initial conditions and observed states. While the laminar state is noted as LS, there exist two unsteady states, US1 and US2, with distinct amplitudes and axial position for the disordered turbulent patch. For LS, C z is steady, i.e. dC z /dt ≈ 0. However, for US, C z is unsteady and larger than the laminar C z , i.e. C z,LS < C z,U S .

In order to describe the unsteady flow patterns, the axial velocity fluctuations, u z , at the pipe centreline are obtained by subtracting the instantaneous axial velocity, u z , to the steady laminar case u 0 z at the same spatial location:

u z (0, 0, z, t) = u z (0, 0, z, t) -u 0 z (0, 0, z), (7)
and plotted in a space-time diagram. The reference value u 0 z is the laminar state with no perturbation, i.e. A = 0. The diagram is constructed from horizontal lines at given times. The colour within the line is the value of u z described by 1500 points along the centreline in the z direction at each 100 time-step. The final diagram is the superposition all computed times with the resolutions of 0.1d in space. This diagram is completed by time evolution of C z graduated at the top of the diagram. Fig. 2 presents u z space-time diagrams and the evolution of drag coefficient, C z , for four different perturbation amplitudes, A, at Re = 2000. The colour intensity in the diagrams increases with the deceleration of the streamwise velocity. The centreline position is typical of this flow configuration and other radial positions close to the centreline would lead to essentially the same qualitative behavior as the localized turbulence fills up the pipe radially. There are two (c), which can move axially along the downstream section.

In Fig. 4, a summary of the simulation cases is presented in the form of a threshold curve where the boundaries between LS and US2 are shown. Systematically, the US1 is between the LS and the US2 indicating the border between laminar and disordered motion is sensitive to initial conditions. The threshold can be described by a power-law fit: A ∝ Re -3 , which is much steeper than the value of -0.006 that is reported by Sanmiguel-Rojas and [START_REF] Sanmiguel-Rojas | Finite-amplitude solutions in the flow through a sudden expansion in a circular pipe[END_REF] indicating that the flow is more and more sensitive as Re increases. Note that the perturbation used by Sanmiguel-Rojas and Mullin ( 2012) is a transverse velocity (tilt), i.e. an addition of a y-transverse velocity component to the parabolic inlet flow, which creates a velocity discontinuity in the inlet section and at the wall. The effects of the generated shear could depend on the mesh resolution as well as on how the solver interpolates the discontinuity. A recent experimental study from Lebon et al. (2018a,b) has shown a power-law of Re -2.3 , which is in closer agreement to the present numerical results. In the experiments, several disturbances mechanism were tested: (i) single, (ii) suction and (iii) periodic inout or synthetic jet from a hole in the wall. Clearly, a direct comparison between the experiments and the DNS is not straightforward for the reason that the vortex perturbations and the jet disturbances are different. Yet, the effect of the vortex perturbation or the jet disturbance is to break the flow symmetry and distort the recirculation region.

Transient growth of unsteady flow patterns

Taking advantage of the time-accurate and 3D opportunities of our DNS results, the unsteady flow patterns are examined. Figure 5(a) presents the space-time diagram for Re = 1360 and A = 0.2, corresponding to case 5a in Table 1. The flow experiences several sequences of laminar, oscillatory and disordered motion, as well as relaminarisation of the turbulent localised patches. In Fig. 5(b), two profiles of u z along the centreline at t = 1000 s, before the emergence of a turbulent patch, and at t = 2500 s, during a turbulent patch, are presented. The profiles can be divided into three regions. In the first region, z < 20d, the flow experiences a gradual deceleration. For 20d < z < 50d, weak wavy unsteady oscillations are observed. The wavelength is ≈ 2.3d. In the third region, z > 50d, the oscillations decay for t = 1000 s, whereas they are considerably amplified for t = 2500 s. The and (d) at z = 80d have a different behaviour, correspond to LS and US2, respectively. Thus, at early stages of development, the differences between the two cases are not particularly striking, so the transition process involves subtle effects. These observations are in disagreement with the transient growth analysis [START_REF] Cantwell | Transient growth analysis of flow through a sudden expansion in a circular pipe[END_REF] in the sense that the disturbance generated at the inlet grows while traveling downstream, mainly because of its interactions with the recirculation region. In the current DNS, the present methodology uses finite amplitude perturbation, which can either decay or grow. Furthermore, the perturbation can remain silent or bounded for a long time and then suddenly grow because of the nonlinear interactions between the unsteady flow and the downstream of the recirculation region, where the reattachment point is located.

The unsteady flow pattern can be also analysed by using a Reynolds averaging technique based on a LS domain. The base flow is approximated by taking the average fields during the LS from t = 1500 to 1800 s: u(x, y, z) = u(x, y, z, t) 1500<t<1800 .

(8)

Then, the unsteady pattern can be expressed using:

u (x, y, z, t) = u(x, y, z, t) -u(x, y, z). (9)
The amplitude of u (x, y, z, t) is small and stays bounded in laminar phase.

For convenience, the coordinate system is converted from Cartesian, u (x, y, z, t), to cylindrical, u (r, θ, z, t). Next, the Fourier transform of u (r, θ, z, t) in azimuthal direction, noted as u , is computed:

u (r, k θ , z, t) = 1 2π 2π 0 u (r, θ, z, t)e -2iπθk θ dθ, (10)
where, the azimuthal modes energy are computed point-wise as:

   e z (r, k θ , z, t) e y (r, k θ , z, t) e z (r, k θ , z, t)    = 1 2    | u x (r, k θ , z, t)| 2 | u y (r, k θ , z, t)| 2 | u z (r, k θ , z, t)| 2    , (11)
and the total energy of each azimuthal modes is the sum of all three directions.

From the linear stability point of view, the perturbation u could be decomposed into growing and decaying modes. However, in our case, its magnitude stays bounded and fluctuates for a long time, before the transient growth starts and US2 emerges. It suggests there should be a slow growing mode. The observation of the time evolution of all azimuthal modes, over all z, r and k θ , reveals no steady increasing of energy. Instead, in addition to the regular unsteady pattern, a perturbation appears in the steady zone, z < 20d, and evolves into turbulence. A closer look into the time evolution of the perturbation amplitude in the vicinity of transition, for 1900 s < t < 2200 s, is shown in Fig. 7 for an arbitrary value of r, here r = d/4, and three different streamwise positions. The other values of r exhibit similar behaviours. The fluctuations of energy in streamwise direction, e z (r, k θ , z, t), or spanwise direction, e xy = e x + e y , are recorded as a signal in time for given values of z, r and k θ . The magnitude of the energy of the azimuthal modes is several orders smaller than the streamwise energy. From these signals, a local energy peak in the time-evolution can be observed and is located close to the transition. This peak can be recorded as a time t p corresponding to the maximum of the energy signal. The values of collected t p will be a function of spatial positions r, z and the mode k θ . By comparing the fluctuation patterns, it is found that fluctuations in streamwise direction are approximately two orders of magnitude stronger than the sum of energy fluctuations in the two other spanwise directions in the laminar region, z < 20d. However, the gap in the energy levels gets closer at the transition point in the second zone: z ≈ 50d and the non-linear region, in the third zone: z > 60d. Another observation is that the position of the energy peaks mainly depends on the streamwise position, z, and seems to be independent of spanwise variables, r and k θ . For z ≈ 17d, the peaks start to be observed at t = 2080 s.

Fig. 8 presents the contour plots at z = 17d of different streamwise velocity in the cross-sectional plane. It is found that the most intense perturbation in the energy peak appears close to the strongest shear rate position in the streamwise velocity mean profile u z , see Fig. 8(b). This suggests the energy peak is related to a shear instability and is independent of the regular unsteady pattern u z . The regular unsteady pattern, see Fig. 8(c), is the fluctuation collected before the peak emerges. The latter structure is different from both the peak velocity and the mean flow profile.

Considering the first four modes (k θ ∈ {0, 1, 2, 3}), as shown in Fig. 7, one can notice the position of local energy peak evolves smoothly and linearly in the streamwise direction in the range of 10d < z < 50d. This suggests that the energy peak is a perturbation that appears, gets carried downstream by the main flow and amplifies. In Fig. 9(a), the tracking of the local peaks velocity in space, z, and time, t, is shown for all the value of r, where a linear fit of z/V 0 + c 1 with c 1 is a constant and V 0 ≈ 0.667. The error of the tracking process mainly comes from the output frequency of the data and the fluctuation in z direction. All the peaks are found in the range 10d < z < 40d. It is interesting to mention that the peak evolution tracking in space, z, and in time, t, collapse for all the values of spanwise variables (r and k θ ). As a consequence the magnitude of the energy peaks over time grows following a power law as shown in Fig. 9(b). The growth for the first four azimuthal modes takes place until z ≈ 50d where the energy saturates.

It is noted that the phenomenon of the convective instability is also observed in the numerical simulations of flow over a backward-facing step (see i.e. [START_REF] Blackburn | Convective instability and transient growth in flow over a backward-facing step[END_REF]).

Hysteresis

In general, a hysteresis appears when two flow solutions can exist for the same Re. Therefore, the initial conditions and the parameters of the disturbance control the appearance of the two solutions. To test the hys- teresic behaviour, two branches of the simulations with Re increasing and Re decreasing are investigated. The Re increasing branch starts from laminar flow, whereas the Re decreasing branch begins from an unsteady state, here US2. This approach was also considered by Sanmiguel-Rojas and Mullin (2012) using a transverse velocity disturbance of amplitude δ. Depending on δ and Re, a domain of hysteresis was observed. Specifically, for δ = 0.001, the coexistence region was reported for 1475 < Re < 1850. Moreover, they found that the hysteresis region grows as δ decreases. It is not possible to directly compare the effect of the transverse velocity disturbance with our vortex disturbance because these perturbations are of different nature: the vortex disturbance perturbation introduces rotation, whereas the transverse velocity disturbance introduces a translation to the flow. However, the following results discuss the universality of the hysteresis behaviour. A series of eight simulations were performed, with a fixed amplitude of vortex perturbation, A = 0.2. The initial condition is case 3a with Re = 1300. When the simulation time reached 1500 s, the final state is analysed and used as the initial condition for the next run at a higher Re. Then, Re, is increased with steps of 25, up to Re = 1400 (case 4a, 5a, 7a, 8a), as depicted in the space-time diagrams of Fig. 10(a). In Fig. 10(b), the decreasing Re branch is initialised with a laminar state (case 3a), then Re is directly increased to 1400 (case 8b) and the decreasing path down to 1300 with a step of 25 (case 7b, 5b, 4b). The simulations from the increasing and decreasing paths are compared for the same Re. The results are presented in Fig. 10(a) and 10(b) and show minor changes. Looking at the drag coefficient also represented in Fig. 10(a) and (b), every change in Re initiates a peak or a transient increase of C z corresponding to a disordered patch that propagates and decays downstream. A small hysteresis is found while following this procedure. Additional simulations were performed using different Re steps and a larger range of Re. The decreasing branch started with Re = 2000. Then Re, is decreased to 1700, 1350 and 1000 (simulation L1d, L2d, L3d). The results is compared with the increasing branch which start from Re = 1000 then increase to 1350, then 1700 (simulation L3d, L2i, L1i). Again, the results are represented in the form of space-time diagrams, in Fig. 11(a) and 11(b). The data show minor differences in the space-time behaviour. However, the drag, C z , now suggests a small loop of hysteresis.

When comparing two series of simulations, one could notice: in the first loop, with small steps of Re (∆Re = 25), the two cases with Re = 1350 are laminar in both increasing and decreasing branches. Whereas, in the second series, with larger steps of Re (∆Re = 350), the two cases with same Re = 1350 show US2. The space-time diagram of 4 cases with Re = 1350 are extracted in the two series and shown in Fig. 12. The different behaviour of the same Re suggests that larger steps of Re could eventually trigger the unsteady behaviour and potentially the hysteresis loop sooner than smaller steps.

To quantify the hysteresis behaviour, S, the integral of the curve C z over Re for both increasing and decreasing branches is defined:

S = C z (Re)dRe. ( 12
)
The relative difference, H, is defined to quantify the hysteresis behaviour:

H = ∆S S . ( 13
)
where ∆S is the difference between two branches and S is the mean value. It is noted that in order to observe the hysteresis phenomenon clearly, a specific procedure needs to be implemented. In the decreasing branch, the variation of Re should be large enough to avoid the transformation from US1 to US2. On the other hand, in the increasing branch, the variation of Re should be small enough to keep the flow laminar. The loops are presented in Fig. 13, with the decreasing branch initiated at Re = 2000, then decreasing consecutively to 1700, 1350 and 1325 (simulation L1d, L2d, L3dbis). The increasing branch is initiated at Re = 1300 and is then increased consecutively to 1325, 1350 and 1700 (simulation 3a, 4a, 10a). Based on the criteria defined in Eq. ( 13), this last procedure leads to a hysteresis of H = 27.87% compared to the two previous loops with the measure of hysteresis are H = 1.15% and H = 0.75% respectively. In Fig. 13, C z is obtained from the last 100 seconds of each case and is plotted against Re. The laminar states, Re = 1325 and Re = 1350, lead to C z with almost the same value within 0.04%. The unsteady states at Re = 1375 and Re = 1400 have slightly different final value of C z within 3% because of the unsteady nature of the flow. The systematic study of all possible steps with extremely long time scales would be a tedious investigation thus beyond the scope of this study. The fact that the hysteresis measure, H, increases with the Re steps height suggests the hysteresis is related to the re-establishment of the recirculation region.

Conclusion

The results of numerical simulations of the flow through a circular pipe with a sudden expansion have been reported. The expansion ratio is 1:2 and the inlet velocity profile is parabolic together with a finite amplitude vortex perturbation. The spatio-temporal velocity fluctuation have been presented when the perturbation amplitude is larger than the threshold. In the subcritical range of Reynolds number, the critical threshold scales with Re -3 , which is consistent with recent experiments (Lebon et al., 2018a,b). The present vortex perturbation, added at the inlet, distorts the flow and the recirculation region before the appearance of disordered motion close to the reattachment point of the laminar recirculation region.

Using the present DNS, it was possible to monitor spatially and temporally the velocity fluctuations, the drag coefficient and the energy of the flow. Furthermore, it was found that the velocity fluctuations in the streamwise direction are dominant, i.e. about two orders of magnitude larger than the transverse flow components. Additionally, a peak of energy fluctuations was observed around to the region of high shear rate close to point of reattachment of the recirculation region and its amplitude grows exponentially along the axial position of the expansion. Finally, a hysteresis quantification procedure, increasing and decreasing the Reynolds number, indicate a small hysteresis region. Its measure depends on the step and computation time, suggesting a transient effect due to the re-establishment of the recirculation region. In the future, it is planned to extend the present simulations to other types of disturbances and to expansions with various expansion ratio and diverging angles [START_REF] Lanzerstorfer | Global stability of the twodimensional flow over a backward-facing step[END_REF][START_REF] Jotkar | Instability mechanisms in straight-diverging-straight channels[END_REF][START_REF] Jotkar | Two-dimensional modal and non-modal instabilities un straight-diverging-straight channel flow[END_REF][START_REF] Kfuri | Friction coefficients for Bingham and power-law fluids in abrupt contractions and expansions[END_REF].

The computations were conducted using HPC resources of the Centre Régional Informatique et d'Applications Numériques de Normandie (CRI-ANN). The authors acknowledge financial support of the Agence National de la Recherche (ANR) thought the programme 'Investissement d'Avenir' from the laboratoire d'excellence Energy Materials and Clean Combustion Center (LabEx EMC3). Our work has also benefited from helpful discussions with Ashley P. Willis (University of Sheffield, UK), who suggested this form for the vortex perturbation.

Figure 1 :

 1 Figure 1: (a) Axial vorticity of the vortex perturbation (R Ω = d/4, (x Ω , y Ω ) = (0, -d/4) and A = 1) and (b) Sketch of fluid domain with 2 cross-sections inlet and outlet mesh.

Figure 2 :

 2 Figure 2: Space-time diagrams of the centreline perturbed streamwise velocity, u z , and drag coefficient, C z , as a (red) line, indicated on the top, for Re = 2000 at (a) A = 0.0782, (b) A = 0.09, (c) A = 0.2 and (d) A = 0.5. All the diagrams use the same colour code shown on the top right corner, such that light colour corresponds to laminar flow. The whole pipe section downstream the expansion, up to z = 150d, is also sketched in the top right corner below the colour bar.

Figure 3 :

 3 Figure 3: Contour plot of instantaneous streamwise velocity at Re = 2000, y = 0 and zoomed into the range of 0 < z < 100d. (a) LS at A = 0.09, t = 750, (b) US2 at A = 0.09, t = 1250 and (c) US1 at A = 0.2, t = 1000.

Figure 4 :

 4 Figure 4: Critical amplitude, A, of the vortex perturbation as a function of Re. The (green) triangles down indicate US and the (blue) triangles up LS. The (red) dashed line is power-law fit: A ∝ Re -3 .

Figure 5 :

 5 Figure 5: (a) Space-time diagram for Re = 1360 and A = 0.2, corresponding to case 5a. (b) u z (0, 0, z, t) profiles along the centreline at t = 1000 (LS ) and t = 2500 (US2 ).

Figure 6 :

 6 Figure 6: u z (0, 0, z, t) signal over time for Re = 1360 and A = 0.2 recorded at z = 40d (a ) and z = 80d (b ). The time range of (a) and (b) are in LS, whereas the time rage in (c) and (d) are in US2, recorded at z = 40d (c ) and z = 80d (d ). The time range in this figure are marked with the corresponding colour in Fig. 5(a).

Figure 7 :Figure 8 :

 78 Figure 7: Time-evolution of the energy of the most relevant azimuthal modes at r = d/4 for Re = 1360 and A = 0.2. The colour code is: blue line for the mode 0, red (thick) line for the mode 1, black line for mode 2, and magenta line for mode 3. The three plots in upper row (a,b,c) show the energy of fluctuation pattern in streamwise direction, e z , and the three plots in the lower row (d,e,f) show the energy of fluctuation pattern in spanwise direction. The three columns correspond to (a,d) z = 17d, (b,e) z = 47d and (c,f) z = 67d. The (green) circles indicate the peaks of energy.

Figure 9 :

 9 Figure 9: Tracking position of the energy peak for Re = 1360 and A = 0.2. (a) Time and axial position with the scaling represented as a (green) line: t p = z/V 0 + c 1 , c 1 is a constant and the different colours represented different radial position. (b) Energy peak amplitude, e p z , versus axial position, z, for the first four modes: mode 0 (×), mode 1 (o), mode 2 ( ), mode 3 ( * ), with the scaling in (green) line e p z = c 2 z 20 where c 2 is a constant.

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure 10: Hysteresis loop of transitional flows for A = 0.2. (a) Space-time diagrams of the increasing Re, (cases 3a, 4a, 6a and 7a) and (b) the decreasing Re, (cases 7b, 6b, 4b and 3b) branches. The diagrams also indicate the value of C z as a (red) continuous line labeled on the top and the corresponding values of Re on the right.

Table 1 :

 1 Summary of the simulations. Abbreviations: LS: laminar state and US: unsteady state. The initial condition could be the Hagen-Poiseuille profile, noted as H-P, or the final times of a past simulation.

	Case	Re	A	Initial condition Remark
	1	1100	0	H-P	LS
	1a	1100 0.458	1	LS
	1b	1100 0.480	1	LS, US2
	1c	1100 0.494	1	US1
	2	1300	0	H-P	LS
	2a	1300	0.2	2	LS
	2b	1300 0.239	2	LS, US2
	2c	1300 0.385	2	US1
	3a	1325	0.2	2a	LS
	4a	1350	0.2	3a	LS
	5a	1360	0.2	4a	LS, US2
	6a	1375	0.2	4a	LS, US2
	7a	1400	0.2	6a	LS, US2
	7b	1400	0.2	2a	LS, US2
	6b	1375	0.2	7b	LS, US2
	5b	1360	0.2	6b	LS, US2
	4b	1350	0.2	6b	LS, US2
	3b	1325	0.2	4b	LS, US2
	2b	1300	0.2	3b	LS, US2
	8	1600	0	H-P	LS
	8a	1600 0.123	8	LS
	8b	1600 0.128	8	LS, US2
	8c	1600 0.16	8	US1
	9a	1700	0.2	9e	US1
	9a	2000 0.0773	9	LS
	9b	2000 0.0782	9	LS, US2
	9c	2000 0.09	9	LS, US2
	9d	2000	0.1	9	US2
	9e	2000	0.2	9	US1
	9f	2000	0.5	9	US1
	10a	1700	0.2	3a	US1
	L1d	1700	0.2	9e	US1
	L2d	1350	0.2	L1d	US2
	L3d	1000	0.2	L2d	LS
	L2i	1350	0.2	L3d	US2
	L1i	1700	0.2	L2i	US1
	L3dbis 1325	0.2	L2d	LS