Seule la force de traînée, introduite à travers le temps de relaxation de la particule, et la pesanteur agissent sur la particule. En régime de Stokes, le temps de relaxation de la particule est donné par l'équation (B.11). Pour traiter des nombres de Reynolds particulaires supérieurs à l'unité, le temps de relaxation est défini par :

$$ \tau_p = \frac{\rho_p}{\rho} \frac{4d_p}{3C_D |\mathbf{v} - \mathbf{v}_p|} \quad \text{avec} \quad C_D = \frac{24}{Re_p} (1 + 0.15Re_p^{0.687}) \quad \text{pour } Re_p < 1000 $$

$$ C_D = 0.44 \quad \text{pour } Re_p \ge 1000 $$

<figure><img src="image_25.png" /><figcaption>Figure 2.5 : Éléments fluides et parcours de la particules solide S [Pozorski 98]</figcaption></figure>

Détaillons l'équation du mouvement (2.15). $\mathbf{v}_p$ est le vecteur vitesse de la particule solide. $\mathbf{v}[\mathbf{x}_p(t), t]$ est la vitesse du fluide porteur, c'est-à-dire la vitesse instantanée de différentes particules fluides, schématisées par F sur la figure 2.5, vue par la particule solide le long de sa trajectoire. Elle correspond à la vitesse du fluide non perturbée par la particule, prise au centre de la particule solide $\mathbf{x}_p(t)$ au temps $t$. La position de la particule n'est fonction que de $t$, dont la valeur est reliée à sa vitesse par :

$$ \frac{d\mathbf{x}_p}{dt} = \mathbf{v}_p \qquad (2.16) $$

$\mathbf{v}_p$ est donc une fonction de la vitesse instantanée du fluide, vue par la particule solide. Il est à noter que la particule doit être suffisamment petite, pour considérer que le champ de vitesse est homogène autour de l'inclusion.

Pour les écoulements turbulents, l'utilisation de l'équation du mouvement de la particule exige en plus une description du champ de vitesse instantanée du fluide. Or, à l'exception de la simulation numérique directe ou de la simulation directe des grandes échelles, les modèles utilisés pour les écoulements turbulents ne donnent accès qu'aux vitesse moyennes. En particulier, la moyenne de Reynolds ou la moyenne spatio-temporelle contribue à une perte d'information sur les grandeurs turbulentes instantanées. Il faut
donc reconstruire le champ turbulent instantané, et cela en générant les fluctuations de vitesse du fluide (selon un processus euléro-lagrangien), vues par la particule. Ce dernier point est la difficulté majeure des modèles lagrangiens pour traiter le transport turbulent de scalaires. L'approche déterministe complète pour la résolution de l'équation (2.15) étant hors de portée à l'heure actuelle, le mouvement turbulent du fluide ainsi que la trajectoire des particules solides sont décrites comme des fonctions régies par des lois probabilistes.

L'approche lagrangienne est utilisée pour traiter la dispersion des particules biotiques dans la couche de surface et dans les couverts végétaux. Il existe dans la littérature de nombreux modèles probabilistes plus ou moins complexes. Tous se basent sur l'équation différentielle stochastique de la trajectoire qui fournit les coordonnées de la particule solide à un instant donné, en fonction de l'instant précédent, par intégration de sa vitesse. Leurs différences proviennent du degré de prise en compte des phénomènes qui affectent la particule dans son mouvement, et en particulier les caractéristiques de la turbulence. De manière succincte, nous pouvons présenter les deux principales approches stochastiques utilisées.

## Les modèles d'interaction tourbillon-particule

Les modèles d'interaction tourbillon-particule considèrent que les particules sont soumises à une succession de tourbillons « aléatoires », caractérisés par leur temps de vie et leur taille (eddy interaction model) [Pétrissans 01]. À chaque pas de temps, une probabilité de mort du tourbillon est calculée, et permet de maintenir ou non la particule dans ce tourbillon. Les fluctuations de la vitesse du fluide sont générées en fonction de cette probabilité (un nouveau tourbillon entraîne de nouvelles fluctuations). Cela revient donc à considérer la durée et la distance en deçà de laquelle la fluctuation est corrélée ou totalement décorrélée d'elle-même. Cette information peut être évaluée (de manière grossière toutefois [Chassaing 00]) par l'échelle intégrale spatiale eulérienne $L_E$ (ou macro-échelle spatiale de Taylor) et l'échelle de temps intégrale lagrangienne $T_L$ qui est une mesure de temps de la persistance des structures tourbillonnaires. Pour les modèles $k - \varepsilon$, ces grandeurs peuvent se déduire par [Rodi 84, Pozorski 98]:

$$L_E = C_E \frac{k^{3/2}}{\varepsilon} \qquad (2.17)$$

$$T_L = C_L \frac{k}{\varepsilon} \qquad (2.18)$$
$C_E$ et $C_L$ sont des constantes à déterminer expérimentalement. Cette mise en oeuvre passe par le calcul de la fluctuation de la vitesse, vue par la particule :

$$u'_{i} = \sqrt{\overline{u'_{i}u'_{i}}} \xi = \sqrt{\frac{2}{3}k} \xi \qquad (2.19)$$

où $\xi$ est un nombre aléatoire à distribution gaussienne de moyenne nulle, et de variance égale à 1. Les autres grandeurs turbulentes sont issues de la résolution des équations du mouvement du fluide, de l'énergie cinétique et du taux de dissipation. Enfin, la trajectoire de la particule est calculée par son équation dynamique.

### Les modèles basés sur l'équation de Langevin

La première formulation de Langevin a été introduite pour modéliser le mouvement brownien d'une particule lourde, immergée dans un fluide, subissant des collisions avec les molécules fluides. L'effet du fluide sur la particule brownienne a été modélisé suivant deux forces distinctes : une force de frottement visqueux et une force fluctuante due aux impacts des particules fluides. Cette force fluctuante, appelée force de Langevin, est indépendante de la vitesse de la particule et n'est fonction que du temps. Par suite, l'équation de Langevin est la première équation différentielle stochastique puisqu'elle contient une force aléatoire.

Cette approche est la plus utilisée pour la dispersion atmosphérique de particules suffisamment petites, pour suivre les fluctuations du fluide, mais qui peuvent toutefois sédimenter à leur vitesse de chute terminale (même restrictions que pour les modèles eu- lériens, que nous verrons à la prochaine section). Le mouvement du fluide ayant une cer- taine mémoire ou persistance, la trajectoire (fluide ou solide) est simulée par une marche pseudo-aléatoire. Le mouvement de la particule est décrit comme une succession de dépla- cements discrets, assimilés à une marche aléatoire possédant une mémoire déduite de la connaissance des propriétés statistiques de la turbulence [Wilson 81, Aylor 90]. La vitesse (ou la fluctuation) d'un élément fluide porteur, prise au pas de temps $t + \Delta t$, est la somme d'un terme de mémoire qui s'exprime en fonction de sa vitesse à l'instant $t$, et d'un terme d'accélération aléatoire, pour rendre compte de la fluctuation turbulente. En turbulence homogène, nous pouvons écrire :

$$u_{i}^{n+1} = \underbrace{u_{i}^{n} R_{L}(\Delta t)}_{\text{mémoire}} + \underbrace{\sigma_{u_{i}} \sqrt{(1 - R_{L}^{2}(\Delta t))} \xi^{n}}_{\text{partie aléatoire}} \qquad (2.20)$$
$\xi^n$ est un nombre aléatoire issu d'un tirage gaussien, et $\sigma_{u_i}$ est l'écart type de $u_i^n$. Le terme de mémoire est représenté par la fonction d'autocorrélation lagrangienne $R_L(\Delta t)$ de $u_i^n$ qui doit s'exprimer comme une fonction de l'échelle de temps intégrale lagrangienne $T_L$ [McCartney 85]:

$$R_L(\Delta t) = \exp\left(-\frac{\Delta t}{T_L}\right) \qquad (2.21)$$

La suite {$u_i^n(t_n)$} est appelée séquence Markovienne (par opposition au processus Markovien) à des instants $t_n$ [Legg 82]. Cette séquence impose au pas de temps $\Delta t$ d'être supérieur au temps pour lequel l'accélération aléatoire reste corrélée à elle-même. Cette échelle de temps n'est autre que l'échelle de temps de Kolmogorov $\tau$, imposée par les structures dissipatives:

$$\tau = (\nu/\epsilon)^{1/2} \ll \Delta t \qquad (2.22)$$

Cette hypothèse fondamentale permet d'appliquer cette séquence à la dispersion turbulente.

Dans le cas de couverts végétaux, la turbulence étant non homogène, la dispersion des particules est fortement influencée par la mémoire ou persistance du mouvement –les particules doivent voyager sur des distances de l'ordre de la hauteur du couvert, pour que le temps de voyage excède $T_L$ [Raupach 89]. Par conséquent, le pas de temps $\Delta t$ doit vérifier [Loubet 00]:

$$\tau \ll \Delta t \ll T_L \qquad (2.23)$$

Nous avons essayé d'entrevoir les paramètres nécessaires pour les modèles lagrangiens. Il apparaît indispensable de disposer d'informations sur les corrélations statistiques du mouvement turbulent (moments d'ordre 3 et plus élevés pour des écoulements inhomogènes et/ou avec de forts cisaillements) et sur le taux de dissipation. En particulier, l'échelle de temps lagrangienne joue un grand rôle dans ces modèles. Malheureusement, cette grandeur est difficilement mesurable, surtout dans les couverts végétaux. De plus, les processus de dépôts au sol ou sur les éléments des végétaux exigent des lois de rebond ou de piégeages spécifiques (fond probabiliste [Aylor 01]). Cependant, contrairement aux modèles gaussiens ou eulériens qui ne prennent pas en compte les phénomènes de remise en suspension par rafales ou qui supposent implicitement que la source émettrice est continue, les modèles lagrangiens peuvent intégrer ces informations, et donc apporter des réponses plus pertinentes [Braaten 90]. Enfin, les modèles lagrangiens n'utilisent pas la modélisation phénoménologique, qui relie le flux turbulent de masse au gradient, relation
qui est souvent contestée (notamment dans le couvert végétal). La dispersion des particules est calculée de manière plus physique, à partir de considérations statistiques sur les champs turbulents, ce qui rend ces modèles de plus en plus appréciés pour la simulation de la dispersion de particules biotiques.

## 2.4 Approche Eulérienne

Contrairement à l'approche lagrangienne, ici l'observateur est fixé en un point et examine les grandeurs physiques du fluide qui passent au point d'observation. La première formulation eulérienne pour les écoulements diphasiques, dont la phase dispersée ne constitue pas un continuum connecté (c'est le cas en particulier pour des particules solides dans l'air) mais dont le fluide porteur est continu, a été proposée par Spalding [Crowe 96]. Bien que l'approche eulérienne semble moins naturelle que la description lagrangienne pour traiter le transport des particules solides, elle correspond souvent à une réalité expérimentale, puisque les grandeurs acquises sont pour la plupart de type eulérien, étant mesurées en un point. L'approche eulérienne est couramment utilisée pour la simulation d'écoulements gaz-solides dilués<fnref n="7" /> et peut être mise en oeuvre à l'aide de plusieurs considérations.

Le modèle diphasique eulérien-eulérien (ou modèle à deux fluides) est un modèle qui distingue les deux constituants, en introduisant une fraction volumique $R_α$. $R_p$ pour la phase dispersée correspond au volume $V_p$ (volume occupé par les particules) sur le volume total $V$ (dans lequel coexistent les particules solides et fluides). De manière similaire, la fraction volumique du fluide porteur est le volume $V_f = V - V_p$ contenu dans $V$. Chaque phase $R_α$ est définie par ses équations de conservation de masse et de quantité de mouvement, les particules ayant une vitesse propre, différente de celle de l'écoulement du fluide porteur. $R_α$ étant une variable turbulente, la résolution des équations du mouvement des particules implique l'introduction d'un modèle de turbulence propre aux particules [Laslandes 98]. Cette approche permet de prendre en compte les interactions telles que les échanges de quantité de mouvement entre les particules et le fluide, l'influence des particules sur la turbulence de l'écoulement du fluide porteur ainsi que les effets diphasiques (effets de ségrégation spatiale des particules et effets de croisement de trajectoire dus à l'inertie des particules).

Le modèle eulérien monophasique s'appuie sur le principe que les particules se com-

<footnote n="7">7. Un écoulement est dilué si $R_p = \frac{N_p V_p}{V} \le 0.001$ où $N_p$ est le nombre de particules, de volume $V_p$, contenu dans le volume $V$</footnote>
portent comme un scalaire passif, au sein de l'écoulement. Cette hypothèse est admise lorsque :

- Les particules considérées ont une densité beaucoup plus grande que celle du fluide et ont un diamètre inférieur ou du même ordre que les plus petites échelles de la turbulence (échelle de Kolmogorov);

- Le nombre de Stokes est bien inférieur à l'unité, c'est à dire que le temps de relaxation des particules solides est petit devant l'échelle de temps de la turbulence.

Ces conditions impliquent que la présence des particules n'a pas d'effet sur la turbulence de l'écoulement du fluide porteur. D'autre part, que les effets d'inertie des particules (et de façon plus générale, les effets diphasiques) sont négligeables : la réponse des particules aux accélérations du mouvement turbulent du fluide porteur est parfaite, et la vitesse relative entre la phase fluide et la phase dispersée est égale à la vitesse de sédimentation. La phase particulaire et la phase fluide sont donc découplées. La démarche à suivre consiste alors à résoudre les champs aérodynamiques (vitesse, turbulence...) du fluide porteur dont la solution permet de traiter les particules.

Pour traiter la dispersion des particules biotiques, nous avons employé la modalité eu- lérienne monophasique, dont les hypothèses sont relativement vérifiées dans l'atmosphère, à l'exception des particules lourdes ($d_p > 500 \mu m$ et $v_s > 1.83 \text{ms}^{-1}$) pour lesquelles les ef- fets d'inertie ne peuvent être négligés [Csanady 63]. Nous avons également restreint notre étude à une couche limite atmosphérique neutre.

## 2.5 Équation instantanée de convection-diffusion

Le champ particulaire, considéré comme continu, est caractérisé par la concentration $C$ définie comme étant le nombre de particules présentes par unité de volume de fluide. L'évolution de la concentration dans l'écoulement est décrite par l'équation de convection-diffusion au sein du fluide porteur :

$$ \underbrace{\frac{\partial C}{\partial t}}_{I} + \underbrace{u_j \frac{\partial C}{\partial x_j}}_{II} - \underbrace{v_s \frac{\partial C}{\partial x_j}}_{III} \delta_{j3} = \underbrace{\frac{\partial}{\partial x_j} \left[ k_c \frac{\partial C}{\partial x_j} \right]}_{IV} + \underbrace{S_p}_{V} \quad (2.24) $$

Le transport des particules solides dans l'écoulement traduit en fait la conservation de la concentration $C$ qui est fonction du temps et des coordonnées spatiales. Les différents
termes sont identifiés par :

-I- La variation temporelle de la concentration locale;

-II- Le transport par convection. La concentration est transportée par le fluide en mouvement, à la même vitesse;

-III- Les particules étant pesantes, la vitesse de sédimentation $v_s$ est ajoutée, pour prendre en compte la vitesse de glissement des particules solides par rapport aux particules fluides. Ainsi, dans la direction de l'écoulement, les particules ont une vitesse qui correspond à celle du fluide, mais peuvent sédimenter sous l'action de la gravité;

-IV- Le transport par diffusion provoqué par l'agitation moléculaire du fluide. Ce transfert de masse est décrit par la loi de Fick, et $k_c$ est le coefficient de diffusion moléculaire de C;

-V- $S_p$ [particule $m^{-3} s^{-1}$] correspond à un terme source, dans le cas d'un apport de particules (le signe est alors positif) ou un terme puits, s'il y a une perte de particules (le signe est alors négatif).

La concentration est considérée comme une variable turbulente et par suite, la prise en compte du transport turbulent se fait à l'aide d'un traitement statistique de Reynolds, dans la couche limite.

### 2.5.1 Dispersion dans la couche limite atmosphérique

À l'exclusion des particules lourdes qui ont une masse qui influence leur mouvement du fait de leur inertie, l'équation de convection-diffusion est valide pour traiter le transport turbulent des particules biotiques telles que le pollen ou les spores [McCartney 90]. Dans ce cas, la variable C admet une décomposition de Reynolds dans la couche limite atmosphérique :

$$C = \bar{C} + C' \qquad (2.25)$$

En n'incluant que des particules dont le diamètre $d_s$ est supérieur à 0.1µm (pas de mouvement brownien) [Chamberlain 67], les termes de diffusion moléculaire sont négligeables, et l'équation moyenne de transport turbulent de la phase particulaire obéit à :

$$\frac{\partial \bar{C}}{\partial t} + \bar{u}_j \frac{\partial \bar{C}}{\partial x_j} - v_s \frac{\partial \bar{C}}{\partial x_j} \delta_{j3} = - \frac{\partial \bar{C}' u_j'}{\partial x_j} + \bar{S}_p \qquad (2.26)$$
En régime turbulent, les différents termes constituant cette équation sont directement transposables de l'équation de convection-diffusion dont nous avons identifié tous les termes. Ici, le terme convectif est un transport par le mouvement moyen. Le terme de corrélation turbulente $\overline{C'u'_j}$ entre la vitesse fluctuante et la concentration fluctuante, correspond à un transfert diffusif par agitation turbulente. C'est donc un flux turbulent de masse qui requiert une modélisation.

L'évaluation des flux turbulents de masse peut s'appuyer sur une modélisation phénoménologique [Legg 79], à condition que les particules considérées aient une vitesse de sédimentation suffisamment petite pour que les effets d'inertie soient négligeables (notamment l'effet de croisement de trajectoire [Csanady 63] qui illustre le fait qu'une particule lourde se déplacera dans plusieurs structures turbulentes au cours de sa trajectoire, alors qu'une particule fluide restera dans la même structure, jusqu'à la dissipation de celle-ci: ceci a pour conséquence une réduction de la diffusivité turbulente de la phase particulaire). De plus, si la vitesse de sédimentation est du même ordre ou plus petite que la fluctuation de la vitesse $w'$ (typiquement de l'ordre de 1.25 $u_*$ [McCartney 90]), la diffusivité de la phase particulière peut être prise identique à celle d'un gaz [Csanady 63]. Compte tenu de l'analogie entre les transferts de type diffusif par agitation turbulente et par agitation moléculaire, le flux turbulent de masse est proportionnel au coefficient de diffusivité turbulente de la concentration $k_{ct}$ et au gradient de $\overline{C}$ dans la direction $j$:

$$ \overline{C'u'_j} = -k_{ct} \frac{\partial \overline{C}}{\partial x_j} \qquad (2.27) $$

La diffusion turbulente de masse est supposée isotrope, contrairement aux modèles de dispersion atmosphérique [Belot 76, Di-Giovanni 90] par exemple, qui considèrent des diffusivités turbulentes différentes suivant les trois directions. Ce choix s'impose de lui-même car le modèle de turbulence utilisé est isotrope (que ce soit le modèle $k-l$ ou le modèle $k-\epsilon$), la viscosité turbulente calculée étant un scalaire.

Par suite, en introduisant le nombre de Prandtl-Schmidt turbulent $\sigma_c$ qui correspond au rapport de la viscosité turbulente $\nu_t$ sur le coefficient de diffusion turbulent de la phase dispersée $k_{ct}$, le flux turbulent de masse devient :

$$ \overline{C'u'_j} = - \frac{\nu_t}{\sigma_c} \frac{\partial \overline{C}}{\partial x_j} \qquad (2.28) $$

Le nombre de Prandtl-Schmidt a la même signification que le nombre de Prandtl turbulent, appliqué à la concentration. La plupart des modèles de dispersion atmosphérique fixent la valeur de $\sigma_c$ à 1 [Aylor 90] pour décrire la diffusivité turbulente verticale
des particules. Par conséquent, les échelles caractéristiques de vitesse et de longueur du mouvement turbulent sont supposées identiques pour les deux phases. Cette hypothèse explique alors toutes les restrictions (que nous venons de citer) imposées pour l'utilisation de l'analogie phénoménologique. Le nombre de Prandtl-Schmidt est pourtant largement utilisé pour traiter le transport de la neige dans l'atmosphère ou le transport hydraulique des sédiments [Celik 88], qui notamment a testé la sensibilité de $\sigma_c$ sur le modèle de transport et conclu que la valeur optimale de $\sigma_c$ est 0.5. Cette valeur laisse entendre que l'efficacité de diffusion turbulente des sédiments est 2 fois supérieure à la diffusion de la quantité de mouvement.

Nous conservons donc le nombre de Prandtl-Schmidt pour la description du flux turbulent de masse, dont la présence permet d'une part, d'étudier son influence sur la dispersion des particules qui nous intéressent, d'autre part, de modifier l'hypothèse d'isotropie de la diffusivité turbulente des particules, en lui affectant différentes valeurs suivant les directions de l'espace. L'équation de convection-dispersion modélisée devient :

$$ \frac{\partial \bar{C}}{\partial t} + \bar{u}_j \frac{\partial \bar{C}}{\partial x_j} - v_s \frac{\partial \bar{C}}{\partial x_j} \delta_{j3} = \frac{\partial}{\partial x_j} \left[ \frac{\nu_t}{\sigma_c} \frac{\partial \bar{C}}{\partial x_j} \right] + \bar{S}_p \quad (2.29) $$

En l'absence de source de particules, $\bar{S}_p$ est nulle dans la couche de surface [McCartney 90]. L'évolution du champ de concentration des particules, près du sol, doit être précisée par des conditions aux limites adéquates et sera abordée à la section 2.6.1 concernant les dépôts.

### 2.5.2 Dispersion dans le couvert végétal

La décomposition de Reynolds ne suffit plus à décrire le vent moyen, ainsi que les champs scalaires moyens (température, concentration...) du fait de leur variabilité spatiale dans le couvert végétal. De plus, le processus de transport turbulent est étroitement lié à la distribution des sources (et puits) des scalaires à la surface des éléments du couvert [Finnigan 87]. Pour introduire ces caractéristiques dans les équations (de conservation de masse, de température...), il est nécessaire d'appliquer une moyenne spatiale, en plus de la moyenne temporelle. La concentration admet une décomposition spatio-temporelle et nous rappelons que la moyenne spatiale est effectuée dans un volume d'intégration contenant au moins un élément du couvert (mais qui est exclu du volume) :

$$ C = \bar{C} + C' = \langle \bar{C} \rangle + \bar{C}'' + C' \quad (2.30) $$
Après un double traitement statistique de l'équation instantanée d'advection-diffusion, l'équation moyenne de transport turbulent d'un champ scalaire dans le couvert végétal, devient [Finnigan 85]:

$$ \frac{\partial \langle \overline{C} \rangle}{\partial t} + \langle \overline{u}_j \rangle \frac{\partial \langle \overline{C} \rangle}{\partial x_j} - v_s \frac{\partial \langle \overline{C} \rangle}{\partial x_j} \delta_{j3} = \frac{\partial}{\partial x_j} \left\langle k_c \frac{\partial \overline{C}}{\partial x_j} \right\rangle - \frac{\partial}{\partial x_j} \left[ \langle \overline{C' u'_j} \rangle + \underbrace{\langle \overline{C'' u''_j} \rangle}_{\text{I}} \right] \\ - \underbrace{\frac{k_c}{V} \iint_{S_I} \frac{\partial \overline{C}}{\partial n} dS}_{\text{II}} \qquad (2.31) $$

Dans cette équation, nous identifions de gauche à droite le terme temporel, le terme de transport par convection et par sédimentation, le transport par diffusion moléculaire et par diffusion turbulente qui s'exprime de manière similaire à l'équation (2.29). En revanche, cette équation diffère de son équation homologue, par la présence des termes apparus lors de la moyenne spatio-temporelle :

-I- Ces termes sont assimilés aux flux dispersifs de masse. Ces flux ont été mesurés dans le cas d'un couvert homogène (pour la température) et sont faibles devant les flux turbulents de masse [Raupach 86]. Bien qu'ils soient susceptibles d'être plus importants lorsque la structure du couvert est fortement hétérogène [Finnigan 87], nous les négligerons devant les flux turbulents de masse;

-II- Ce terme est fondamental. Il rend compte de l'activité des éléments végétaux sur le champ scalaire moyen, dans le volume d'intégration. Les éléments végétaux peuvent ainsi modifier le bilan de conservation du champ moyen, en l'augmentant (source) ou en le diminuant (puits):

$$ \langle \overline{S_p} \rangle = - \frac{k_c}{V} \iint_{S_I} \frac{\partial \overline{C}}{\partial n} dS \qquad (2.32) $$

$\langle \overline{S_p} \rangle$ dépend évidemment du champ scalaire étudié. Pour les particules solides, il va englober les processus de dépôt et de remise en suspension qui peuvent intervenir entre les éléments du couvert et le microclimat. En particulier, en l'absence de source, $\langle \overline{S_p} \rangle$ n'est pas nul contrairement à $\overline{S_p}$, car les particules peuvent chuter en impactant sur les différents organes de la plante et sédimenter sur des éléments végétaux, tels que les feuilles par gravité [Legg 79, Aylor 01]. Ce terme représente donc les échanges nets de particules entre le couvert et l'air. Il va dépendre de la cadence de remise en suspension, ainsi que de
la cadence de dépôt des particules. La modélisation de ce terme fera l'objet de la section suivante.

En supposant que $\langle \overline{S_p} \rangle$ soit connu, le problème de fermeture consiste maintenant à modéliser les flux turbulents de masse $\langle \overline{C'u'_j} \rangle$. Ce flux peut s'exprimer suivant une relation flux-gradient, dont le coefficient de diffusivité turbulent est :

$$k_{ct} = - \langle \overline{C'u'_j} \rangle \left( \frac{\partial \langle \overline{C} \rangle}{\partial x_j} \right)^{-1} \qquad (2.33)$$

$k_{ct}$ est par définition positif et s'exprime à l'aide de $u'$ et $l'$ qui sont les échelles caractéristiques du mouvement turbulent, soit $k_{ct} = u'l'$. Dans le cas d'un couvert végétal, cette relation de flux-gradient (2.33) est régulièrement non vérifiée. En effet, des études expérimentales ont montré que le comportement de $k_{ct}$, à proximité du couvert, diffère suivant trois régions :

- Pour $z > 2h$, $k_{ct}$ et $\nu_t$ sont donnés par les relations de similitude de Monin-Obukhov où les échelles de longueur et de vitesse sont respectivement $z-d$ et $u_*$ dans le cas d'un écoulement en stratification neutre. De plus, $k_{ct}=\nu_t$, donc dans la sous-couche inertielle, les diffusivités turbulentes s'écrivent :

$$k_{ct}^* = \nu_t^* = \kappa u_* (z-d) \qquad (2.34)$$

- Pour $h \le z \le 2h$, il a été observé que les diffusivités turbulentes sont supérieures à leur valeur dans la sous-couche inertielle, au fur et à mesure que l'on se rapproche du couvert [Finnigan 00]. En particulier, la diffusivité turbulente des scalaires est 2 à 3 fois plus importante que celle de la quantité de mouvement :

$$\nu_t = 1.1 \nu_t^*$$

$$k_{ct} = 2 \nu_t^* \text{ à } 3 \nu_t^*$$

$$\sigma_{ct} = 0.5 \text{ à } 0.3$$

- Dans le couvert végétal $z < h$, les observations de $\nu_t$ et de $k_{ct}$ sont très erratiques, notamment dans les couverts forestiers, avec de nombreuses singularités, pouvant donner lieu à des diffusivités nulles ou négatives [Raupach 89, Raupach 96, Finnigan 00]. En effet, dans un couvert de pin, Denmead & Bradley (1985) ont constaté des flux turbulents de chaleur et de diooxide de carbone dirigés dans le sens des gradients [Kaimal 94]. Ces flux turbulents à «contre-gradient» mettent en défaut la relation
flux-gradient pour décrire le transport turbulent par diffusion, qui requièrent l'introduction d'une autre relation de fermeture, par exemple, en utilisant des moments d'ordre supérieur.

<table><caption>Tableau 2.3 : Comparaison des propriétés moyennes de la couche de surface, de la couche de mélange et de l'écoulement dans un couvert végétal proche de h d'après [Raupach 96]. T = transport par le mouvement turbulent</caption><thead><tr><th>Propriétés</th><th>Couche de surface</th><th>Couche de mélange</th><th>Couvert en z = h</th></tr></thead><tbody><tr><td>Point d'inflexion</td><td>non</td><td>oui</td><td>oui</td></tr><tr><td>P<sub>rt</sub><sup>-1</sup> = k<sub>t</sub>/ν<sub>t</sub></td><td>1.1</td><td>2</td><td>2</td></tr><tr><td>Bilan de k</td><td>0 ≈ P - ε</td><td>0 = P + T - ε</td><td>0 = P + T - ε</td></tr></tbody></table>

L'origine de ces anomalies s'explique par les échelles mises en jeu pour les transferts turbulents dans le couvert, par leur caractère très intermittent et par la présence de multitudes d'échelles de longueurs et de vitesses qui gouvernent l'écoulement [Finnigan 87]. Nous avons évoqué lors de la modélisation du mouvement dans le couvert végétal, la présence de structures cohérentes propres à la turbulence de couvert. L'analogie avec la couche de mélange (tableau 2.3) permet de mettre en évidence le rôle des différentes structures tourbillonnaires (eddies) présentes dans le couvert, et cela en introduisant l'échelle de cisaillement :

$$L_s = \bar{u}(h) \left(\frac{\partial \bar{u}}{\partial z}\right)_{z=h}^{-1} \qquad (2.35)$$

Cette échelle est identifiée comme l'échelle de longueur verticale à laquelle s'associent de forts transferts turbulents. $L_s$ est de taille comparable à $h$, et plus précisément, de l'ordre de $0.1h$, $0.5h$ et $h$ pour des couverts denses, modérés et clairsemés. Les tourbillons de tailles supérieures à $L_s$ constituent une turbulence inactive et ne peuvent contribuer que faiblement aux transferts verticaux. Les tourbillons de tailles comparables à $L_s$ vont, en revanche, constituer la turbulence active dans le couvert et ainsi contribuer fortement aux échanges turbulents. Enfin, les fines échelles inférieures à $L_s$, dont le principal rôle est de dissiper l'énergie du mouvement turbulent, n'ont pas un poids significatif sur les échanges verticaux [Raupach 96].
La turbulence active à l'origine des transferts verticaux est donc caractérisée par des structures cohérentes de l'ordre de la taille de $h$. Ces structures sont de tailles supérieures aux échelles de variation verticale des scalaires et ne peuvent, de ce fait, décrire ces transferts localement [Raupach 81]. Par ailleurs, [Chassaing 00] ne parle plus de transfert, mais de transport turbulent et définit les structures cohérentes, en reprenant la citation de Roland Stull: *Les structures cohérentes sont des entités qui assurent un transport (de quantité de mouvement, de scalaire...) de type non diffusif sur des distances finies.*

En ce qui concerne l'intermittence des transferts, c'est essentiellement la turbulence inactive (grandes échelles) qui va être à l'origine des structures cohérentes au sein du couvert, et cela, par le biais de courtes, intermittentes mais extrêmes rafales (gust ou sweep). Par conséquent, chaque rafale à grande échelle va initier un «paquet» de structures cohérentes, qui proviennent de l'instabilité du point d'inflexion du vent au sommet du couvert, séparées d'une distance moyenne de $\Lambda_x \approx 8L_s$ [Raupach 96, Kaimal 94]. C'est essentiellement lors de ces rafales, que les transferts turbulents (de masse, de quantité de mouvement) efficaces ont lieu.

Des techniques d'analyse par quadrant seuil, ont permis d'affiner les connaissances des transferts turbulents dans les couverts. Cette méthode consiste à représenter le signal turbulent dans l'espace des phases ($u'w'$) (mais reste applicable dans le cas du flux de masse ou de chaleur), chaque quadrant étant attribué à un évènement spécifique [Finnigan 00]:

Quadrant 1 $u' > 0 \ w' > 0$ interaction externe

Quadrant 2 $u' > 0 \ w' < 0$ balayage

Quadrant 3 $u' < 0 \ w' < 0$ interaction interne

Quadrant 4 $u' < 0 \ w' > 0$ éjection

Ces études ont mis en évidence la présence des structures cohérentes, ainsi que leurs rôles importants dans l'expansion et le mélange à travers deux mécanismes: la majeure contribution des phénomènes de balayage ou injection, qui induit des transferts turbulents vers le bas (quadrant 2). Viennent en suite les phénomènes d'éjection qui sont à l'origine de transferts ascendants (quadrant 4). Enfin, les interactions internes et externes interviennent peu dans les transferts turbulents.

Bien que la relation flux-gradient semble inadaptée pour la description des transferts turbulents dans le couvert, elle reste néanmoins très utilisée [Belot 76, Aylor 89, Aylor 99]. Dans de nombreux modèles eulériens de dispersion, la diffusivité turbulente est spécifiée
<figure><img src="image_26.png" /><figcaption>Figure 2.6 : Vitesse moyenne et viscosité turbulente dans et au-dessus d’un couvert de blé en stratification neutre ($u_*$ = 0.43 $ms^{-1}$, $h$ = 1 m et le coefficient d’at-ténuation $\alpha$ = 4). La diffusivité turbulente des spores est égale à 3 fois la viscosité turbulente, qui suit une loi exponentielle $\nu_t = \nu_t(h)\exp[\alpha(z/h - 1)]$ pour $0 < z < h/2$ et $\nu_t = \nu_t(h)$ pour $h/2 < z < h$ d’après [Aylor 89]</figcaption></figure>

(c'est-à-dire qu'elle n'est pas explicitement calculée), suivant une loi exponentielle ana-logue à la vitesse dans le couvert (figure 2.6). Legg & Powell (1979) ont notamment obtenu des résultats satisfaisants pour la modélisation de la dispersion de spores *Lycopodium* dans et au-dessus d'un champ d'orge [Legg 79].

En considérant que notre modélisation de l'écoulement dans la canopée végétale permet de retranscrire correctement les flux turbulents, nous avons utilisé la relation flux-gradient pour décrire les flux turbulents de masse. L'équation de dispersion dans le couvert végétal s'exprime alors de façon similaire à son équation (2.29) homologue, valable dans la couche de surface :

$$ \frac{\partial \langle \overline{C} \rangle}{\partial t} + \langle \overline{u_j} \rangle \frac{\partial \langle \overline{C} \rangle}{\partial x_j} - v_s \frac{\partial \langle \overline{C} \rangle}{\partial x_j} \delta_{j3} = \frac{\partial}{\partial x_j} \left[ \frac{\nu_t}{\sigma_c} \frac{\partial \langle \overline{C} \rangle}{\partial x_j} \right] + \langle \overline{S_p} \rangle \quad (2.36) $$

Le terme source ou puits $\langle \overline{S_p} \rangle$ nécessite une modélisation, car il n'est pas nul en l'absence d'émission de particules dans le couvert. Par ailleurs, le nombre de Prandtl-Schmidt est utilisé afin de tenir compte de l'augmentation de la diffusivité, par rapport à la viscosité turbulente vers le sommet du couvert, dont la valeur est proche de $\sigma_c = 0.5$ (tableau 2.3).
## 2.6 Les processus de dépôt

Nous avons évoqué à plusieurs reprises l'importance des processus de dépôts dans la couche de surface et dans le couvert végétal. Ces processus de dépôts sont essentiellement dus à la sédimentation et à l'impact. Que ce soit dans la couche de surface ou dans le couvert végétal, l'équation de dispersion s'exprime de manière similaire, à ces différences près :

- Pour la dispersion à l'échelle de la couche de surface, les conditions de dépôt sont prises en compte par le biais d'une condition aux limites adéquate. Cette condition doit décrire globalement les processus de dépôt qui sont accentués par la présence des éléments rugueux dans la sous-couche rugueuse (sans résoudre à cette échelle);

- Dans le cas du couvert végétal, les éléments constituant le couvert interviennent de façon explicite dans les échanges particules-air.

Bien que ces processus soient relativement bien appréhendés d'un point de vue physique, ils restent néanmoins complexes, et leur description repose sur des hypothèses simplificatives.

### 2.6.1 Les dépôts dans la couche limite atmosphérique

Les difficultés majeures rencontrées par les modèles de dispersion (et cela reste vrai pour toutes les approches) sont les conditions limites, notamment près du sol. En effet, l'équation de conservation de la concentration ne tient pas compte des processus de dépôt au voisinage du sol, et des conditions aux limites doivent être introduites, afin de préciser l'évolution de la phase particulaire près du sol.

Afin d'écrire la condition limite au sol, considérons l'équation de convection-diffusion, dans le cas d'un écoulement développé stationnaire, à une distance $z_p$, près du fond du domaine. L'équation se réduit alors à faire un bilan de flux, qui se met sous la forme générale suivante :

$$v_s \overline{C} - \overline{w'C'} = F_d - F_e = F_n \quad (2.37)$$

$v_s \overline{C}$ est associé au flux de particules dû à la gravité et représente donc un dépôt par sédimentation et $\overline{w'C'}$ est le flux turbulent de masse. $F_d$ et $F_e$ sont respectivement le flux de dépôt et le flux d'érosion [particule $s^{-1}m^{-2}$]. Cette équation représente donc le flux
net de particules à travers la surface. Nous pouvons distinguer deux situations :

- Lorsque le flux net est nul ($F_n = 0$), le régime d'équilibre est atteint. Les échanges de masse dans le plan horizontal sont nuls, car ils se compensent exactement. Le profil de concentration ne varie plus le long de l'écoulement;

- En revanche, si le flux net est différent de zéro, c'est-à-dire en situation de non équilibre, il y a bien un échange. Il est alors nécessaire de spécifier les flux de dépôt, ainsi que les flux d'érosion. Il apparaît de façon claire que c'est la situation de non équilibre qui nous intéresse.

Nous négligerons les flux d'érosion ($F_e = 0$) puisque nous ne traitons pas la remise en suspension dans l'air des particules déposées, mais nous signalons qu'ils sont souvent décrits en terme de contrainte seuil de mise en mouvement (dont nous avons vu que ce n'est pas approprié pour les particules biotiques) ou en terme de capacité de transport de l'écoulement [Celik 88]. Sous ces conditions, l'équation (2.37) se simplifie, et le flux de dépôt doit être précisé.

Chamberlain a étudié le dépôt de *Lycopodium* pour différentes configurations d'écoulement. Il a constaté que le flux de dépôt s'exprime en fonction de la concentration moyenne, prise à une distance $z_p$ de l'ordre de $h$ (ou $10 z_0$) et d'une vitesse $v_d$ appelée vitesse de dépôt. La condition limite au sol, dans la couche de surface, devient donc :

$$ \frac{\nu_t}{\sigma_c} \frac{\partial \bar{C}}{\partial z} + v_s \bar{C} = F_d = v_d \bar{C} \quad (2.38) $$

$v_d$ est la vitesse intégrale de dépôt dans la sous-couche rugueuse.

La vitesse de dépôt caractérise les dépôts qui s'effectuent par sédimentation mais également par impact inertiel sur les éléments de rugosité. Par conséquent, $v_d$ dépend de nombreux facteurs (tableau 2.4) tels que :

- La taille de la particule. Ceci implique que $v_d$ est au moins de l'ordre de la vitesse de sédimentation et McCartney utilise une vitesse de dépôt [McCartney 91]:

$$ v_d = 2v_s \quad (2.39) $$

- La structure et la densité du couvert. Cela caractérise donc l'influence de la surface sur les processus de dépôt [Yao 97];

$$ v_d = v_s (1 + LAI) \quad (2.40) $$
<table><caption>Tableau 2.4 : Vitesse de dépôt et vitesse de sédimentation pour différentes particules sur des surfaces naturelles, d'après [McCartney 85]</caption><thead><tr><th>Particules</th><th>Taille [µm]</th><th>Surface</th><th>v<sub>d</sub> [cms<sup>-1</sup>]</th><th>v<sub>s</sub> [cms<sup>-1</sup>]</th><th>u/u<sub>*</sub> [ms<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Spore Lycopodium</td><td>32</td><td>Herbe courte</td><td>0.7</td><td>1.9</td><td>0.12 u<sub>*</sub></td></tr><tr><td>Spore Lycopodium</td><td>32</td><td>Herbe courte</td><td>3.5</td><td>1.75</td><td>1.1-3.6 u</td></tr><tr><td>Spore Lycopodium</td><td>32</td><td>Orge</td><td>1.4-2.1</td><td>1.8</td><td>-</td></tr><tr><td>Pollen Zea maïs</td><td>90-100</td><td>Sol rugueux</td><td>63 ± 25</td><td>20</td><td>-</td></tr><tr><td>Gouttelettes</td><td>30-60</td><td>Blé (0.8 m)</td><td>6.8 ± 0.9</td><td>4.9</td><td>0.37 u<sub>*</sub></td></tr><tr><td>Gouttelettes</td><td>30-60</td><td>Blé (1 m)</td><td>21 ± 5</td><td>10.5</td><td>0.39 u<sub>*</sub></td></tr></tbody></table>

LAI est l'indice de surface foliaire qui représente la surface foliaire totale par unité de surface du sol [$m^2 m^{-2}$].

- La vitesse de l'écoulement et la turbulence, qui influencent les échanges avec la surface rugueuse. En particulier $v_d/v_s$ augmente avec la vitesse de frottement $u_*$ [Chamberlain 67]. Ces résultats sont confirmés par Callender, qui définit la vitesse de dépôt de gouttelettes de 30 à 60 µm de diamètre dans un champ de blé, comme une fonction de $u_*$:

$$v_d = v_s + 0.1 u_* \qquad (2.41)$$

Cette vitesse peut également traiter le dépôt des spores *Lycopodium* dans un champ de blé humide [McCartney 85].

Généralement, $v_d$ est supérieur à la vitesse de sédimentation, de l'ordre de $2v_s < v_d < 5v_s$, sauf si l'écoulement est peu turbulent et que les processus de dépôt sont principalement gouvernés par la gravité (dans ce cas $v_d = v_s$). Pour traiter le dépôt des particules, il est nécessaire d'utiliser une vitesse de dépôt. Chaque auteur met en avant un paramètre pour expliciter la vitesse de dépôt, dont l'expression dépend de nombreux paramètres. À l'heure actuelle, aucune n'est pleinement satisfaisante et de nombreux modèles tentent de calculer explicitement cette vitesse [Slinn 82, Ferrandino 85] ou de la mesurer pour des couverts et des expériences variés.
## 2.6.2 Les dépôts dans le couvert végétal

En l'absence de source, $\langle \overline{S_p} \rangle$ représente les processus de dépôt qui ont lieu dans le couvert. Ces processus de dépôt ont été identifiés comme :

- La sédimentation, qui découle principalement de la vitesse de sédimentation $v_s$ des particules et de la distribution spatiale du feuillage. Ce processus est donc gouverné par la gravité;

- Le dépôt par impact inertiel sur les éléments végétaux, qui dépend de la structure du couvert, de la vitesse de l'écoulement et de la turbulence, ainsi que du temps de relaxation des particules. Ce processus de dépôt est gouverné par l'inertie de la particule.

$\langle \overline{S_p} \rangle$ peut s'exprimer suivant une loi linéaire de la concentration locale moyenne, dont la constante de proportionnalité est :

$$ \langle \overline{S_p} \rangle = r \langle \overline{C} \rangle \qquad (2.42) $$

$r$ est la cadence de dépôt [$s^{-1}$] qui dépend de la structure du couvert, de la taille des particules et des conditions météorologiques (vent, turbulence...).

Sous cette forme, l'équation (2.42) englobe les deux mécanismes de dépôt, qui ne sont pas dissociés : le dépôt dépend fortement de la valeur de $r$. Lorsque le vent est relativement faible, la prise en compte du dépôt par $r$ est acceptable, puisque la sédimentation est prépondérante devant l'impact. Cependant, pour des vents plus forts, et en particulier dans le cas de rafales qui génèrent une turbulence active pour les transferts de masse, les particules sont violemment éjectées contre les éléments du couvert, et l'impact inertiel ne peut plus être ignoré.

Afin de dissocier et de discerner chaque processus de dépôt, Legg & Powell ont supposé que le dépôt par sédimentation concerne principalement les surfaces horizontales des végétaux (les feuilles) et que les particules chutent, à la suite de leur impact, sur les éléments végétaux perpendiculaires à l'écoulement (troncs, par exemple) [Legg 79]. Cette hypothèse implique implicitement des processus de dépôt indépendants. En réalité, au contact des feuilles, en particulier, la vitesse du vent s'annule et une couche limite se développe. Au sein de cette couche limite, il se crée une composante verticale de la vitesse de l'écoulement, qui affecte le dépôt par sédimentation : ce processus n'est donc pas découplé des conditions micro-climatiques. Cependant, très peu de théories permettent de traiter
<figure><img src="image_27.png" /><figcaption>Figure 2.7 : Processus de dépôt des particules sur un élément d’un couvert végétal assimilé à un cylindre</figcaption></figure>

ce problème, et le taux de dépôt net par unité de volume s’exprime simplement comme la somme d’un terme qui résulte du dépôt par sédimentation $\langle \overline{S_s} \rangle$ et d’un terme de dépôt par impact $\langle \overline{S_i} \rangle$ [Legg 79, Aylor 01]:

$$ \langle \overline{S_p} \rangle = \langle \overline{S_s} \rangle + \langle \overline{S_i} \rangle \qquad (2.43) $$

Avec :

$$ \langle \overline{S_s} \rangle = E_s f_z A v_s \langle \overline{C} \rangle \qquad (2.44) $$

$$ \langle \overline{S_i} \rangle = E_i f_x A \langle \overline{u} \rangle \langle \overline{C} \rangle \qquad (2.45) $$

$E_s$ et $E_i$ (sans dimension) désignent l’efficacité de sédimentation et l’efficacité d’impact. $f_z$ et $f_x$ sont respectivement la fraction de la densité de surface foliaire $A$ dans le plan de normale $\vec{n}_z$ et de normale $\vec{n}_x$ (figure 2.7). Ils dépendent donc de l’orientation des éléments du couvert et sont variables suivant la hauteur, dans le cas d’un couvert hétérogène. Enfin, la vitesse de sédimentation $v_s$ permet de décrire le dépôt par sédimentation, alors que le dépôt par impact s’exprime en fonction de la vitesse locale moyenne $\langle \overline{u} \rangle$.

En supposant que la structure du couvert soit connue, l’évaluation de $\langle \overline{S_p} \rangle$ passe par la détermination des efficacités de collecte par sédimentation et par impact des éléments du couvert.
## Efficacité de sédimentation

Le processus de dépôt gravitationnel est généralement considéré comme optimal. Ceci implique donc une efficacité de collection par sédimentation très proche de l'unité :

$$E_s = 1$$

Toutes les particules qui arrivent sur une surface foliaire par sédimentation y adhèrent et ne sont pas remises en suspension. Cette approximation s'appuie sur l'hypothèse que la sédimentation est un processus suffisamment lent pour que la surface humide retienne la particule [Aylor 78, Legg 79, Aylor 01].

## Efficacité d'impact

Le processus de dépôt par impact dépend de l'efficacité d'impact de l'obstacle étudié. Cette efficacité est définie par [Chamberlain 67, Bache 79] :

$$E_i = \frac{Q_I}{Q}$$

- $Q_I$ est la quantité de particules impactées par seconde sur la surface de l'obstacle perpendiculaire à l'écoulement;

- $Q$ est la quantité de particules qui traversent la surface perpendiculaire à l'écoulement en l'abscence de l'obstacle;

- Lorsque $E_i$ est égale à 1, toutes les particules impactent sur la surface.

De nombreuses études ont été consacrées à la détermination de l'efficacité de collecte par impact pour des spores, sur des plantes dont les organes sont assimilés à des cylindres (pour les tiges), des plaques ou des rubans (pour les feuilles), en soufflerie [Chamberlain 67] ou dans des champs [Bainbridge 79]. Elles révèlent que $E_i$ est une fonction de la dimension et de l'orientation de l'obstacle, mais également de la vitesse de l'écoulement et de la taille de la particule. Pour des nombres de Reynolds (basés sur une dimension caractéristique de l'obstacle) supérieurs à 10, l'efficacité de collection par impact inertiel dépend principalement du nombre de Stokes [Bainbridge 79]. $E_i$ est une fonction non-linéaire croissante (figure 2.8) du nombre de Stokes calculé, en régime de Stokes, par l'expression :

$$S_t = \frac{\tau_p}{\tau_f} = \frac{v_s \langle \bar{u} \rangle}{g l} \quad \text{avec} \quad \tau_p = \frac{v_s}{g} \qquad (2.46)$$
<figure><img src="image_28.png" /><figcaption>Figure 2.8 : Efficacité d'impact sur un cylindre placé perpendiculairement à l'écoulement en fonction du nombre de Stokes, pour des nombres de Reynolds supérieurs à 100. Expérience de May & Clifford (1967) [Aylor 81]</figcaption></figure>

$l$ est une longueur caractéristique de l'obstacle (un diamètre pour une tige ou la largeur pour une feuille) et $\langle \bar{u} \rangle$ est la vitesse moyenne locale de l'écoulement, à proximité de l'obstacle.

En reprenant les études expérimentales de May & Clifford (1967), Aylor propose une expression du coefficient d'impact pour des cylindres et des rubans [Aylor 82, Aylor 01]:

$$E_i = \frac{0.86}{1 + 0.442 S_t^{-1.967}} \qquad (2.47)$$

L'équation (2.47) suppose que toutes les particules qui impactent la plante sont retenues, ce qui semble effectivement vérifié pour la majorité des pollens et des spores. Cependant, certaines spores (E. Graminis) se déplacent par bonds successifs sur des surfaces suffisamment sèches. Par suite, ces mouvements réduisent l'efficacité de collection par impact de la plante, que l'équation (2.47) tend à surestimer. Enfin, chaque élément végétal de la plante n'a pas la même efficacité de collection. Le dépôt par impact s'explicite alors, comme la somme des contributions des différents éléments $n$ (tige, feuille, épi...) constituant le végétal:

$$\langle \overline{S_i} \rangle = \sum_n E_{in} f_{xn} A_n \langle \overline{u} \rangle \langle \overline{C} \rangle \qquad (2.48)$$

$E_{in}$ est donc l'efficacité de collection de chaque organe $n$ de densité de surface $A_n$.

Par ailleurs, [Aylor 81] a étudié le dépôt par impact inertiel dans le cas de rafales de vent. De cette étude, il a déduit que l'efficacité d'impact $E_i$, basée sur la vitesse moyenne
de l'écoulement $\langle \bar{u} \rangle$, tend à sous-estimer les dépôts. Par conséquent, il a défini une efficacité moyenne d'impact pour traiter les phénomènes de rafales :

$$ \bar{E}_i = \frac{\int_0^{u_{max}} E_i(u) f(u) du}{\int_0^{u_{max}} f(u) du} \qquad (2.49) $$

$u_{max}$ est la vitesse maximale observée durant l'expérience, et $f(u)$ est la densité de probabilité de la distribution de la vitesse instantanée du vent, qui suit une loi de Gumbel dans un couvert de blé [Shaw 79].

Cet aspect très intéressant ne sera pas traité. Cependant il peut facilement s'intégrer par la suite dans notre modèle. Pour expliciter le dépôt des particules dans le couvert végétal, nous le décrirons comme la somme des dépôts par sédimentation et par impact (relations (2.43), (2.44), (2.45)) en admettant une efficacité de sédimentation optimale et une efficacité d'impact fonction du nombre de Stokes, donnée par l'équation (2.47).
## Conclusion

Nous venons de présenter notre modèle de dispersion dans la couche de surface en stratification neutre. Ce modèle se base sur l'équation moyenne turbulente de conservation du champ de concentration des particules. Il permet de traiter la dispersion de particules pesantes, dont les effets d'inerties sont négligeables. Près du sol, les processus de dépôt sont spécifiés à cette échelle, par des conditions aux limites utilisant la vitesse de dépôt.

L'originalité de notre travail est d'utiliser la même équation de transport pour traiter la dispersion à l'échelle d'un couvert végétal, et cela en décrivant l'activité des éléments végétaux sur le champ scalaire. En effet, dans le couvert, le dépôt des particules est gouverné par deux processus : la sédimentation qui dépend principalement des caractéristiques de la particule, et de manière cruciale, le dépôt par impact sur les éléments végétaux. Ces processus sont explicitement calculés dans le couvert, et s'appuient sur la définition de l'efficacité de sédimentation et de l'efficacité d'impact de chaque élément du couvert.

Néanmoins, nous avons utilisé une relation de flux-gradient pour exprimer les flux turbulents de masse. Or, cette relation est contestée dans le couvert, car la turbulence est caractérisé par de grosses structures cohérentes qui assurent un transfert turbulent, relevant plus d'un transport convectif que diffusif. Ce point délicat qui nécessitera une attention particulière, sera abordé dans le chapitre consacré aux validations.

Enfin, le modèle de dispersion, que ce soit à l'échelle du couvert ou à plus grande échelle, est couplé au modèle dynamique. De par sa construction, il permet de modéliser des paysages hétérogènes (transitions couvert-sol nu, couvert en rangs) et cela en spécifiant ou non les paramètres du couvert. Il permet donc de répondre à nos objectifs.

Le système final d'équations ne possédant pas de solutions analytiques connues, l'objet du chapitre suivant est de présenter le modèle numérique que nous utilisons.
<figure><img src="image_29.png" /></figure>
# Chapitre 3

## Méthodologie numérique du code de calcul

### Introduction

Nous avons présenté les phénomènes à prendre en compte pour la dispersion des particules dans le couvert et dans la couche limite atmosphérique, à l'aide de formalismes mathématiques. Cependant, ces équations n'ont pas de solutions analytiques connues. Pour résoudre notre problème physique, nous utiliserons l'outil de calcul numérique Aquilon. Ce code de calcul est apte à traiter les équations générales de la mécanique des fluides.

L'objet de ce chapitre est d'exposer les concepts et les méthodes numériques disponibles dans le code de calcul. En particulier, nous pouvons citer les points suivants :

- Le maillage est dissocié de la géométrie étudiée;

- Les équations sont discrétisées en utilisant la méthode des volumes finis;

- La méthode du Lagrangien augmenté permet de découpler les champs de vi- tesse et de pression.

Les équations de la mécanique des fluides disponibles dans le code de calcul, ont été modifiées et adaptées aux écoulements atmosphériques. Le modèle de canopée a été développé, en accord avec les procédures propres au code de calcul. Enfin, nous présenterons les conditions limites qui ont été rajoutées pour modéliser un écoulement atmosphérique en configuration hétérogène, ainsi que les conditions limites pour les particules.
<figure><img src="image_30.png" /></figure>
### 3.1 D'un problème continu à un problème discret

Les équations étant fortement non-linéaires et couplées, il est impossible d'évaluer exactement les solutions aux problèmes posés. Dans ces conditions, nous utilisons une grille de discrétisation afin d'approcher la solution associée à cette grille finie qui constitue notre domaine d'étude.

Pour cela, la méthodologie numérique utilisée, est celle développée pour le code de calcul Aquilon du laboratoire Master. Ce code de calcul est largement appliqué à la résolution de problèmes en mécanique des fluides. De ce fait, il dispose de nombreuses méthodes numériques développées sur des critères de stabilité et de solidité. Avant de présenter les techniques numériques employées pour la résolution des systèmes d'équations physiques de notre problème, nous allons préciser les concepts fondamentaux du code de calcul, ainsi que les modèles de turbulence retenus.

#### 3.1.1 Une grille unique

Une particularité du code de calcul est l'utilisation d'une grille de discrétisation dissociée de la géométrie étudiée. La définition d'une configuration particulière revient à imposer à chaque noeud de la grille pré-définie, une perméabilité, une viscosité ou une caractéristique quelconque du milieu. Ces paramètres physiques sont variables dans l'espace et dans le temps, ce qui permet une étude sur le domaine hétérogène entier. Pour cela, un terme de Darcy de pénalisation a été rajouté aux équations de Naviers-Stokes. Le code de calcul permet de traiter des milieux fluide/solide par le biais d'un coefficient de perméabilité $K$: le milieu est fluide lorsque la perméabilité $K$ tend vers l'infini, et solide lorsque $K$ tend vers 0.

C'est en s'appuyant sur ce principe, qui à l'origine a été développé pour traiter des hétérogénéités fluide/solide, que nous avons développé un système d'équations, dit à l'échelle du paysage et à l'échelle locale. En fait, les équations sont semblables à la différence près :

- À l'échelle locale, le domaine d'étude est un milieu de porosité variable. L'utilisation du maillage déconnecté de la géométrie étudiée permet de traiter, dans sa globalité, le milieu hétérogène fluide/couvert. Chaque maille de calcul est dotée d'une densité de surface végétale et d'un coefficient de traînée;

- À l'échelle du paysage, le couvert n'existe plus, l'écoulement n'est pas explicitement résolu à cette échelle. Par conséquent, les paramètres physiques du couvert tendent
vers 0.

<figure><img src="image_31.png" /><figcaption>Figure 3.1 : Concept de grille unique appliqué aux couverts végétaux</figcaption></figure>

Le modèle numérique, ainsi défini, permet de simuler les phénomènes physiques à différentes échelles, sans aucune modification des équations, et cela par l'ajout ou par la pénalisation des termes spécifiques à l'échelle du phénomène observé.

De manière similaire, la dispersion atmosphérique des particules est traitée dans le milieu fluide/couvert, en affectant à chaque maille de calcul, un terme puits qui modélise les processus de dépôt par impact et par sédimentation et qui ont lieu sur les éléments du couvert.

### 3.1.2 Choix des modèles de turbulence et récapitulatifs du système d'équations final

Nous avons adapté le code de calcul pour la simulation des écoulements atmosphériques aux différentes échelles, en 2 et 3 dimensions. Nous résumons ci-dessous les équations retenues et les modèles de turbulence choisis. Il est à remarquer que de nombreux modèles de turbulence sont disponibles dans le code de calcul : $k - \varepsilon$ étendu, RNG, $v^2 - f$ [Durbin 91]. Bien que ces modèles soient plus élaborés, notre choix s'est porté sur les modèles $k - l$ et $k - \varepsilon$ standard.

En effet, les processus turbulents sont complexes dans le couvert végétal. Il nous a semblé plus intéressant, dans un premier temps, de juger les capacités de tels modèles pour la simulation de la turbulence de couvert. Afin d'être cohérent dans notre démarche, le choix du modèle de turbulence à l'échelle du paysage, est dicté par le modèle de turbulence choisi à l'échelle du couvert.
À l'échelle du paysage

- Pour la modélisation du vent, nous utilisons l'équation de continuité (1.24), les équations de Reynolds (1.31) et l'équation de conservation de la température (1.26) présentées au premier chapitre aux sections 1.3.2 et 1.3.3;

- Pour le mouvement turbulent, nous avons utilisé le modèle $k-l$ (équations (1.37), (1.38) et constantes présentées dans le tableau 1.1 à la section 1.3.3) que nous avons développé, ainsi que le modèle $k-\epsilon$ standard (équations (1.32), (1.40), (1.41) et constantes présentées dans le tableau 1.2 à la section 1.3.4) qui a été validé en 2 et 3 dimensions sur le cas-test complet de l'écoulement turbulent dans un canal derrière une marche descendante [Glockner 00].

- Le transport turbulent des particules est traité par une équation de convection-diffusion (équation (2.29) section 2.5.1) dont l'évolution aux bords du domaine doit être précisée.

À l'échelle du couvert

Nous avons développé le modèle de dispersion atmosphérique à l'échelle du couvert en stratification adiabatique. Nous avons donc discrétisé les termes induits par la présence du couvert, dont les caractéristiques sont introduites dans chaque volume élémentaire :

- Le vent est modélisé par le biais de l'équation de conservation de masse et de quantité de mouvement :

$$ \frac{\partial \langle \overline{u_j} \rangle}{\partial x_j} = 0 $$

$$ \rho \frac{\partial \langle \overline{u_i} \rangle}{\partial t} + \rho \langle \overline{u_j} \rangle \frac{\partial \langle \overline{u_i} \rangle}{\partial x_j} = -\frac{\partial}{\partial x_i} \left[ \langle \overline{P} \rangle + \frac{2}{3} \rho k \right] + \frac{\partial}{\partial x_j} \left[ (\mu + \mu_t) \left( \frac{\partial \langle \overline{u_i} \rangle}{\partial x_j} + \frac{\partial \langle \overline{u_j} \rangle}{\partial x_i} \right) \right] \\ - \rho C_d a \sqrt{\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle \langle \overline{u_i} \rangle} $$

- Le modèle de turbulence de canopée est un modèle de type $k-\epsilon$. Les termes de production de sillage et de dissipation de sillage pour $k$ et $\epsilon$ ont été rajoutés. Une étude de sensibilité sur la constante $C_{D\epsilon w}$, présentée au chapitre consacré aux valuations, a permis de fixer sa valeur à $C_{D\epsilon w} = 4.45$; les autres constantes empiriques étant conservées et exposées dans le tableau 1.5 à la section 1.6.4.
Les équations d'évolution de $k$ et de $\varepsilon$ s'écrivent :

$$ \rho \frac{\partial k}{\partial t} + \rho \langle \overline{u_j} \rangle \frac{\partial k}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \mu + \frac{\mu_t}{\sigma_k} \right) \frac{\partial k}{\partial x_j} \right] + \rho P_k - \rho \varepsilon \\ + C_{Pw} \rho C_d a (\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle)^{3/2} - C_{\varepsilon w} \rho C_d a \sqrt{\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle} k $$

$$ \rho \frac{\partial \varepsilon}{\partial t} + \rho \langle \overline{u_j} \rangle \frac{\partial \varepsilon}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \mu + \frac{\mu_t}{\sigma_\varepsilon} \right) \frac{\partial \varepsilon}{\partial x_j} \right] + C_{\varepsilon 1} \rho \frac{\varepsilon}{k} P_k - C_{\varepsilon 2} \rho \frac{\varepsilon^2}{k} \\ + C_{P\varepsilon w} \rho \frac{\varepsilon}{k} C_d a (\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle)^{3/2} - C_{D\varepsilon w} \rho C_d a \sqrt{\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle} \varepsilon $$

- L'équation de transport du champ de concentration a été modifiée pour prendre en compte les processus de dépôt dans le couvert à travers $\langle \overline{S_p} \rangle$ (sans remise en suspension ni source, ce terme est négatif) :

$$ \frac{\partial \langle \overline{C} \rangle}{\partial t} + \langle \overline{u_j} \rangle \frac{\partial \langle \overline{C} \rangle}{\partial x_j} - v_s \frac{\partial \langle \overline{C} \rangle}{\partial x_j} \delta_{j3} = \frac{\partial}{\partial x_j} \left[ \frac{\nu_t}{\sigma_c} \frac{\partial \langle \overline{C} \rangle}{\partial x_j} \right] + \langle \overline{S_p} \rangle $$

$\langle \overline{S_p} \rangle$ est le dépôt par unité de volume et de temps, par les mécanismes de dépôts par sédimentation et par impact :

$$ \langle \overline{S_p} \rangle = E_s f_z A v_s \langle \overline{C} \rangle + E_i f_x A \langle \overline{u} \rangle \langle \overline{C} \rangle $$

L'efficacité de sédimentation est considérée comme optimale, égale à l'unité, et l'efficacité d'impact s'écrit suivant le formalisme de [Aylor 82] (équation (2.47) à la section 2.6.2).

- Enfin, la viscosité turbulente s'exprime par la relation suivante :

$$ \mu_t = \rho C_\mu \frac{k^2}{\varepsilon} \quad \text{et} \quad \nu_t = \frac{\mu_t}{\rho} $$

Le système d'équation obtenu est donc un système fermé. Rappelons que toutes ces équations sont aptes à simuler la dispersion atmosphérique des particules à deux échelles d'observation, sans aucune modification. Dans le cas d'une modélisation à plus grande échelle, les paramètres physiques du couvert tendent vers 0. Enfin, il est nécessaire d'expliciter nos conditions aux limites, pour l'ensemble des équations.
### 3.1.3 Traitement des conditions aux limites

Le concept de maillage unique requiert la possibilité de pouvoir imposer à chaque itération, et en tout point du domaine de calcul, des conditions de contrôle de la maille. Les équations aux limites sont intégrées aux équations grâce à des termes supplémentaires. Concrètement, cela se traduit par l'apport d'une condition aux limites de flux surfacique de type Fourier, qui s'écrit en notant par $\phi$ une des variables du problème, et $\Gamma$ le coefficient de diffusion associé :

$$ -\Gamma \left. \frac{\partial \phi}{\partial \mathbf{n}} \right|_{paroi} = B i_{\phi} (\phi_{paroi} - \phi_{\infty}) + \varphi \quad (3.1) $$

Où $B i_{\phi}$ représente un coefficient d'échange surfacique propre à la grandeur $\phi$ considérée, $\phi_{\infty}$ une valeur de référence adaptée au problème, $\varphi$ un flux surfacique et $\mathbf{n}$ le vecteur normal unitaire sortant à la paroi.

Suivant les valeurs du nombre $B i_{\phi}$ sur les noeuds frontières du domaine, les conséquences sur la condition à la limite imposée, seront les suivantes :

$$
\begin{array}{lll}
B i_{\phi} = 0 & \text{Condition de Neumann homogène} & \frac{\partial \phi}{\partial \mathbf{n}} = 0 \\
\\
B i_{\phi} \to \infty & \text{Condition de Dirichlet} & \phi_{paroi} = \phi_{\infty} \\
\\
0 < B i_{\phi} < \infty & \text{Condition de Fourier généralisée} & \\
\\
B i_{\phi} = 0 \text{ et } \varphi \neq 0 & \text{Flux surfacique} &
\end{array}
$$

Ce traitement des conditions aux limites est fondé sur la définition physique des coefficients d'échanges thermiques<fnref n="1" />. Il est adapté à toutes les équations de notre système, et inclus directement dans le système à résoudre. Pour l'équation du mouvement, il convient de remarquer qu'il s'agit d'une équation vectorielle qui porte sur toutes les composantes du vecteur vitesse. Chaque composante peut être affectée d'une condition limite différente.

De plus, ces conditions aux limites sont étendues à l'ensemble du domaine de calcul; c'est-à-dire qu'il est possible d'imposer une condition de Dirichlet ou de Neumann à n'importe quel endroit de la grille. En revanche, la condition de Fourier généralisée ne peut s'appliquer que sur le bord du domaine de calcul.

L'écriture des conditions aux limites de notre problème sera abordée à la section consacrée à cet effet, car ces conditions sont particulièrement importantes pour la modélisation d'un paysage hétérogène.

<footnote n="1">1. $\Phi_{paroi} = h(T_{paroi} - T_{\infty})$. $\Phi_{paroi}$ est la densité de flux de chaleur et $h$ le coefficient d'échange convectif [$W m^{-2} K^{-1}$]</footnote>
## 3.2 Procédure numérique des équations

Ces équations qui régissent notre problème vont faire l'objet de techniques numériques. Afin de rechercher la solution sur le domaine de calcul, nous utiliserons une série de méthodes itératives emboîtées.

Pour expliciter la méthodologie employée, nous illustrerons à chaque stade par un exemple en se limitant aux équations de Reynolds (écrites sous forme classique et sans les opérateurs moyennes pour plus de clarté) et de l'énergie. Cependant, cette méthodologie s'applique à l'ensemble des équations retenues.

Les équations de conservation de masse, de quantité de mouvement et de l'énergie auxquelles nous avons ajouté le traitement propre aux conditions limites, s'écrivent :

$$ \nabla \cdot \mathbf{u} = 0 \qquad (3.2) $$

$$ \begin{aligned} \rho \left[ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right] = & - \nabla \left[ P + \frac{2}{3} \rho k \right] + \nabla \cdot \left[ (\mu + \mu_t) (\nabla \mathbf{u} + {}^t \nabla \mathbf{u}) \right] \\ & + \rho \mathbf{g} - \rho \mathbf{L} \mathbf{u} - B i_{\mathbf{u}} [\mathbf{u} - \mathbf{u}_{\infty}] \end{aligned} \qquad (3.3) $$

où $\mathbf{F} = \rho \mathbf{L} \mathbf{u}$ est une force volumique qui s'exprime en fonction de la vitesse. Dans le cadre de notre travail, elle correspond à la force totale due à l'absorption de la quantité de mouvement par le couvert (où $\mathbf{L}$ est en fait $\mathbf{L}(\mathbf{u})$).

$$ \rho C_P \left[ \frac{\partial T}{\partial t} + (\mathbf{u} \cdot \nabla) T \right] = \nabla \cdot \left[ \left( \lambda + C_P \frac{\mu_t}{P_{rt}} \right) \nabla T \right] - B i_o [T - T_{\infty}] \qquad (3.4) $$

### 3.2.1 Discrétisation temporelle

Les équations précédemment définies sont tout d'abord discrétisées en temps, avant d'en effectuer une résolution spatiale. Un traitement implicite est effectué pour les termes formant le système d'équations à résoudre. Ces termes sont tous évalués à l'instant $n+1$, à l'exception de ceux issus de la discrétisation temporelle et des termes faisant apparaître une non linéarité.

L'intervalle de temps, sur lequel nous observons les phénomènes, est discrétisé en un nombre fini d'intervalleles $[t_n, t_{n+1}]$. L'approximation de la variable $\phi$ à l'instant $n$ est notée $\phi^n$.
Pour atteindre la solution stationnaire du problème, il est nécessaire d'effectuer plusieurs itérations. À chacune de ces itérations, la solution réelle est avancée temporellement de la quantité $\Delta t = t_{n+1} - t_n$ appelée pas de temps.

Lors de la phase de discrétisation en temps, la dérivée temporelle $\frac{\partial \phi}{\partial t}$ à l'instant $t_{n+1}$ est alors obtenue à partir d'un développement en série de Taylor tronqué à un ordre $N$. Le calcul de cette dérivée à l'itération $t_{n+1}$ à partir des valeurs aux temps précédents peut se faire de différentes façons :

- En utilisant les résultats de l'itération $n$. C'est le schéma d'Euler, basé sur un développement décentré de Taylor à l'ordre $N = 1$:

$$ \left(\frac{\partial \phi}{\partial t}\right)^{n+1} = \frac{\phi^{n+1} - \phi^n}{\Delta t} \qquad (3.5) $$

- En utilisant les résultats des itérations $n$ et $n-1$, suivant le schéma de Gear, basé sur un développement décentré de Taylor à l'ordre $N = 2$:

$$ \left(\frac{\partial \phi}{\partial t}\right)^{n+1} = \frac{3\phi^{n+1} - 4\phi^n + \phi^{n-1}}{2\Delta t} \qquad (3.6) $$

Le schéma d'Euler est principalement utilisé pour les écoulements stationnaires et le schéma de Gear pour les écoulements instationnaires.

Ces traitements implicites ou semi-implicites permettent de s'affranchir de critères de stabilité du type CFL liant le pas de temps au pas d'espace. De ce fait, la solution d'un problème stationnaire s'obtient rapidement en utilisant un grand pas de temps. Enfin, contrairement aux schémas explicites, qui de plus sont moins robustes, le choix du pas de temps pour les schémas implicites est lié et adapté à la physique des phénomènes, sans tenir compte de considérations d'ordre numérique.

### 3.2.2 Découplage et linéarisation des équations

Lorsque les équations d'un problème physique sont fortement couplées, la résolution globale d'un système unique, issu de la discrétisation en bloc de toutes les équations, est souvent difficilement réalisable.

Pour pallier cette difficulté, il est nécessaire de découpler les équations entre elles. Pour cela, il est possible de réaliser des itérations de couplage, qui consistent à évaluer par exemple le couple $(\mathbf{u}, P)$ à l'itération $n+c/C$ ($c$ étant la $c^{ème}$ itération de couplage et $C$ le nombre d'itérations total de couplage à effectuer) en fonction de la troisième variable
$T$ prise à l'itération précédente $n+(c-1)/C$. À l'aide du nouveau champ calculé ($\mathbf{u}, P$), il reste à évaluer la variable $T$ à l'instant $n+c/C$. Le processus est mené jusqu'à convergence et stabilisation des champs entre deux itérations successives de couplage.

Toutefois, dans certaines conditions, en particulier lorsque les effets de compressibilité sont négligeables, et que les caractéristiques physiques (viscosité) du fluide sont indépendantes (ou faiblement dépendantes) de la température, une seule itération de couplage est réalisée.

La linéarisation des équations consiste à écrire les termes non-linéaires à deux temps différents. Cette procédure est utilisée pour les termes de transport par convection, pour l'ensemble des grandeurs, mais s'applique également aux termes sources et puits non-linéaires de chaque équation.

À l'issue de ces transformations, le système discrétisé devient :

- Résolution à l'instant $n + 1$ du couple $(\mathbf{u}, P)$:

$$ 
\begin{aligned}
\rho^n \left[ \frac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\Delta t} + (\mathbf{u}^n \cdot \nabla) \mathbf{u}^{n+1} \right] &= -\nabla \left[ P^{n+1} + \frac{2}{3} \rho^n k^n \right] \\
&\quad + \nabla \cdot \left[ (\mu + \mu_t^n) (\nabla \mathbf{u}^{n+1} + {}^t \nabla \mathbf{u}^{n+1}) \right] \\
&\quad + \rho^n \mathbf{g} - \rho^n \mathbf{L}^n \mathbf{u}^{n+1} - B i_u^n [\mathbf{u}^{n+1} - \mathbf{u}_\infty] \quad (3.7)
\end{aligned}
$$

- La solution est ré-injectée dans l'ensemble des équations du système, afin d'obtenir la résolution à l'instant $n + 1$:

$$ 
\begin{aligned}
\rho^n C_P \left[ \frac{T^{n+1} - T^n}{\Delta t} + (\mathbf{u}^{n+1} \cdot \nabla) T^{n+1} \right] &= \nabla \cdot \left[ \left( \lambda + C_P \frac{\mu_t^n}{P_{rt}^n} \right) \nabla T^{n+1} \right] \\
&\quad - B i_o^n [T^{n+1} - T_\infty] \quad (3.8)
\end{aligned}
$$

Les équations sont ainsi découplées et linéarisées. Il est à noter que les équations de l'énergie cinétique moyenne du mouvement turbulent et du taux de dissipation sont résolues de manière couplée.
### 3.2.3 Couplage vitesse-pression : méthode du Lagrangien Augmenté

Un des problèmes rencontrés lors de la résolution des équations de Naviers-Stokes (ou de Reynolds) concerne le traitement du couplage entre la vitesse et la pression. En effet, il n'existe pas d'équations explicites permettant le calcul de la pression. Toute la difficulté réside dans la détermination de la pression sous la contrainte de la relation de conservation de la masse.

Nous utiliserons la méthode du Lagrangien Augmenté pour découpler les champs de vitesse et de pression [Khadra 99]. Cette méthode itérative consiste, à partir d'un champ de pression initial, à résoudre les équations de Naviers-Stokes modifiées par l'introduction de la contrainte de divergence nulle. Un nouveau champ de pression est calculé à l'aide d'une équation de continuité modifiée. Le calcul du couple vitesse-pression est effectué par l'algorithme itératif d'Uzawa [Fortin 82].

Les itérations du Lagrangien Augmenté, notées $k$, sont répétées jusqu'à ce que la valeur moyenne de la divergence de la vitesse dans l'ensemble du domaine soit suffisamment petite (dont la valeur est fixée par l'utilisateur). L'algorithme est présenté ci-dessous :

- Initialisations : $\mathbf{u}^{k=0} = \mathbf{u}^n$ et $P^{k=0} = P^n$

- Pour $k = 0$ à $K - 1$ faire :

* Calcul de $\mathbf{u}^{k+1}$ solution de l'équation :

$$ 
\begin{aligned}
\rho^n \left[ \frac{\mathbf{u}^{k+1} - \mathbf{u}^n}{\Delta t} + (\mathbf{u}^k \cdot \nabla) \mathbf{u}^{k+1} \right] &= -\nabla \left[ P^k + \frac{2}{3} \rho^n k^n \right] \\
&\quad + \nabla \cdot \left[ (\mu + \mu_t^n) (\nabla \mathbf{u}^{k+1} + {}^t \nabla \mathbf{u}^{k+1}) \right] \\
&\quad + \rho^n \mathbf{g} - \rho^n \mathbf{L}^k \mathbf{u}^{k+1} - B_i^n \mathbf{u}^{k+1} [\mathbf{u}^{k+1} - \mathbf{u}_\infty] \\
&\quad + dr \nabla (\nabla \cdot \mathbf{u}^{k+1})
\end{aligned}
$$

* Calcul de $P^{k+1}$ à partir de l'équation :

$$ P^{k+1} = P^k - dp \nabla \cdot \mathbf{u}^{k+1} \quad (3.9) $$

- Solutions : $\mathbf{u}^{n+1} = \mathbf{u}^{k=K}$ et $P^{n+1} = P^{k=K}$
Les paramètres *dp* et *dr*, positifs, sont évalués selon des critères de convergence. Cependant, le choix de ces valeurs est parfois délicat et en pratique, plusieurs tests sont nécessaires pour trouver un couple adéquat au problème physique étudié. Enfin, cet algorithme permet de s'affranchir de conditions aux limites sur la pression.

### 3.2.4 Discrétisation spatiale : méthode des volumes finis

À la suite de cette série de transformations des équations différentielles, nous obtenons un système d'équations découplées qu'il reste à discrétiser sur une grille pré-définie.

La structure du maillage que nous utilisons est un maillage cartésien à pas variables, décalé en vitesse (maillage type MAC) et composé d'une superposition de trois grilles en deux dimensions.

<figure><img src="image_32.png" /><figcaption>Figure 3.2 : Maillage 2D du domaine de calcul. Grille principale, grille décalée</figcaption></figure>

Les champs scalaires sont disposés sur le maillage principal, tandis que les deux composantes de la vitesses sont calculées sur les grilles décalées (figure 3.2). L'utilisation d'un
maillage décalé permet non seulement un calcul facilité du gradient de pression et de la divergence de la vitesse, mais évite en plus des oscillations en pression observées pour des maillages non décalés, ce qui tend à rendre le schéma numérique moins stable.

La discrétisation spatiale des équations sur le maillage se fait par la méthode des volumes finis [Patankar 80]. La méthode des volumes finis est basée sur la recherche des valeurs des inconnues du problème sur un nombre fini de points du domaine de calcul. Les équations, écrites sous forme conservative<fnref n="2" />, sont discrétisées de façon à ce que le flux soit conservé lors du passage d'un volume de contrôle à l'autre. Cette manière de procéder fait des volumes finis une méthode naturelle pour résoudre les équations de la mécanique des fluides.

Le domaine de calcul est découpé en volumes de contrôle, associés à chaque noeud du maillage. Les équations sont ensuite intégrées sur le volume de contrôle, centré en chacune des mailles. À ce stade, aucune approximation n'est faite, et l'écriture des équations dans leur forme globale permet de s'affranchir d'un ordre de dérivation par rapport à l'écriture locale.

Les équations sont ensuite discrétisées en terme de flux physiques à travers les interfaces des volumes de contrôle en utilisant le théorème de Green-Ostrogradski. Il est à noter que cette formulation permet d'assurer les propriétés de conservation. Au niveau des interfaces des volumes de contrôle, les variables ne sont pas définies. Le problème consiste à reconstruire leur valeur aux interfaces. Pour cela, il est nécessaire de disposer d'un schéma d'interpolation, qui peut être différent selon les termes considérés (termes diffusifs ou convectifs). L'approximation de la solution se fait par des fonctions affines par morceaux entre les noeuds du maillage.

Les termes de diffusion sont discrétisés avec un schéma centré d'ordre 2. En ce qui concerne la discrétisation spatiale des termes de convection, la précision de l'ordre des schémas et leur stabilité jouent un rôle essentiel à grand nombre de Reynolds. Nous disposons de plusieurs schémas d'interpolation dans notre code de calcul :

- Le schéma Centré, qui évalue la variable en $i + \frac{1}{2}$ à partir d'une formulation de Taylor tronquée à l'ordre 1. Cela revient à faire une moyenne pondérée des deux valeurs adjacentes au point considéré (figure 3.3a). Bien que ce schéma soit le plus naturel et le plus simple, il est rapidement instable et inadapté aux écoulements à fort nombre de Reynolds ;

<footnote n="2">2. Une équation sous forme conservative s'écrit $\frac{\partial \rho \phi}{\partial t} + \frac{\partial}{\partial x_j} (\rho u_j \phi) = \frac{\partial}{\partial x_j} \left( \Gamma \frac{\partial \phi}{\partial x_j} \right) + S$</footnote>
<figure><img src="image_33.png" /><figcaption>Figure 3.3 : Schémas numériques pour l’évaluation de la variable à l’interface $i + \frac{1}{2}$</figcaption></figure>

- Le schéma Upwind ou simple-amont, l'interpolation consiste à remplacer la valeur inconnue en $i + \frac{1}{2}$ par celle située au noeud en amont de l'écoulement (figure 3.3b). Ce schéma d'ordre 1 implique une diffusion numérique importante, mais reste adapté aux écoulements pour lesquels les termes convectifs sont prépondérants;

- Le schéma Hybride s'appuie sur plusieurs schémas d'interpolation. Son principe est de choisir un schéma afin de s'adapter au mieux aux caractéristiques locales de l'écoulement. Pour cela, un nombre sans dimension est introduit et permet d'évaluer localement le rapport qui existe entre les termes convectifs et les termes diffusifs. Ce nombre appelé nombre de Peclet de maille est défini sur l'intervalle d'espace $[X_{i+1}, X_i]$ par :

$$Pe = \frac{\rho u (X_{i+1} - X_i)}{\Gamma} \qquad (3.10)$$

Lorsque le nombre de Peclet tend vers 0, les termes diffusifs sont prépondérants. Dans ce cas, le schéma hybride utilise un schéma centré. En revanche, lorsque les termes convectifs sont prépondérants, le schéma est identique au schéma Upwind. Finalement, le schéma hybride est défini par :

-2 ≤ Pe ≤ 2 schéma centré

$Pe > 2$ schéma upwind

$Pe < -2$ schéma upwind
- Le schéma Quick est décentré d'ordre 2. L'évaluation de l'inconnue sur la face $i + \frac{1}{2}$ est basée sur une interpolation quadratique à partir des deux noeuds situés en amont et d'un noeud situé en aval de l'interface [Parneix 95] (figure 3.3c). Ce schéma diffuse faiblement à l'approche de forts gradients, mais provoque des oscillations ;

- Enfin, un schéma de type Lax-Wendroff TVD d'ordre 2 qui s'avère particulièrement adapté aux équations d'évolution des fronts et qui est peu sensible au problème de diffusion numérique et au problème d'oscillation [Vincent 99].

Chacun de ces schémas est une approximation et le fait d'interpoler des valeurs dé- nature les équations originelles. Il a été montré par des développements de Taylor que l'utilisation de tels schémas équivaut à résoudre une équation à laquelle est ajoutée un terme de diffusion plus ou moins important, selon l'ordre de développement du schéma.

Au vu des caractéristiques numériques des schémas que nous venons de présenter, nous avons choisi un schéma hybride pour la discrétisation spatiale des équations du mouvement et des équations de $k$ et $\varepsilon$. Ce schéma présente un bon compromis entre la précision du calcul et sa rapidité de convergence.

Pour les équations d'advection-diffusion des scalaires, nous avons testé différents sché- mas. Bien que le schéma TVD semble le plus approprié pour le transport des particules, puisque très peu sensible au problème de diffusion numérique qui apparaît à grand nombre de Reynolds, nous avons préféré un schéma Hybride, qui présente une rapidité de conver- gence appréciable, surtout pour les simulations tridimensionnelles.

### 3.2.5 Résolution numérique des systèmes linéaires

La discrétisation implicite des équations engendre des systèmes linéaires (dont la ma- trice est creuse) de taille importante, proportionnelle au nombre de points de discrétisa- tion. Les méthodes itératives classiques (Jacobi, Gauss-Seidel, relaxation) possèdent un taux de convergence lent du fait de la détérioration de leur conditionnement avec l'aug- mentation de la dimension de la matrice.

Nous utiliserons un solveur itératif de type bi-gradient conjugué stabilisé (Bi-CGSstb) [VanderVorst 92] utilisé dans le code de calcul du laboratoire. Afin d'améliorer l'effi- cacité de la résolution de grands systèmes linéaires, les matrices sont préconditionnées [Khadra 99] avec un préconditionnement de type MILU ou ILU. Ce dernier permet d'ac- célérer la convergence des calculs et s'avère fort utile pour des écoulements complexes.
En plus du test sur la valeur moyenne de la divergence, il existe un critère de convergence qui consiste à calculer la norme $L^2$ de la différence entre deux pas de temps d'une variable considérée. Ce test d'arrêt est donc adapté au problème stationnaire. Le choix des variables sur lesquelles il s'effectue, ainsi que la précision désirée est introduit par l'utilisateur. La fréquence (en itération) à laquelle nous effectuons ce test est également commandée par l'utilisateur. Enfin, lorsque la précision souhaitée est atteinte, le calcul s'arrête.

## 3.3 Conditions limites au sol

La résolution du système d'équations nécessite la connaissance des conditions aux limites. Dans cette partie, nous présenterons les conditions aux limites qui ont été développées pour permettre d'intégrer des paysages hétérogènes et les perturbations associées dans la couche de surface en stratification neutre. Puis, nous généraliserons ces relations pour un écoulement en stratification thermique.

Nos conditions aux limites, à proximité d'une paroi, utilisent le concept de loi de paroi [Launder 74] avec une implémentation comparable à [Schiestel 93, Costes 96]. L'écoulement n'est pas résolu explicitement en région proche de la paroi, les effets de surface étant reproduits globalement par le biais de lois empiriques. Il est toutefois possible d'adapter le modèle de turbulence pour résoudre explicitement l'écoulement dans la sous-couche rugueuse -en excluant les couverts végétaux. Bien que cette démarche permette d'accéder à de plus amples informations, elle nécessite des maillages particulièrement fins. Compte tenu de nos objectifs de modéliser des paysages hétérogènes de caractéristiques très variables, l'utilisation des lois de parois à cette échelle est la plus adaptée. Enfin, précisons que l'ensemble des lois de parois définies s'appuie sur le concept des $Bi_φ^3$.

### 3.3.1 Loi de paroi dynamique

La condition limite à la paroi est construite à partir des relations d'équilibre pour une couche à flux constant, sur un sol homogène infini, que nous avons établies précédemment. Chaque cellule proche de la paroi est donc affectée de paramètres physiques tels que la rugosité dynamique et la hauteur de déplacement.

Afin d'écrire les conditions aux limites pour les composantes de la vitesse à la paroi,
nous avons considéré un écoulement en équilibre local avec la surface. Le premier noeud de calcul est pris à une hauteur suffisamment éloignée du sol pour que les flux turbulents soient dominants. Sous ces conditions :

- La vitesse moyenne horizontale suit une loi logarithmique :

$$ \bar{u}(z) = \frac{u_*}{\kappa} \ln \left( \frac{z-d}{z_0} \right) \qquad (3.11) $$

La vitesse de frottement peut donc se déduire :

$$ u_* = \frac{\kappa \bar{u}(z)}{\ln \left( \frac{z-d}{z_0} \right)} \qquad (3.12) $$

- L'hypothèse de couche à flux constant considère que le flux turbulent de quantité de mouvement est constant et égal à sa valeur près du sol :

$$ -\rho \overline{u'w'} = \rho u_*^2 \qquad (3.13) $$

- Les flux turbulents de quantité de mouvement suivent une modélisation phénoméno- logique. Par conséquent, nous pouvons relier les tensions de Reynolds au cisaillement local :

$$ -\rho \overline{u'w'} = \mu_t \frac{\partial \bar{u}}{\partial z} = \rho u_*^2 \qquad (3.14) $$

- L'écoulement est en équilibre local avec la paroi (production=dissipation), ce qui se traduit par :

$$ k = C_{\mu}^{-1/2} u_*^2 \qquad (3.15) $$

Nous pouvons donc exprimer à nouveau la vitesse de frottement en fonction de l'énergie cinétique moyenne du mouvement turbulent :

$$ u_* = C_{\mu}^{1/4} k^{1/2} \qquad (3.16) $$

En combinant les équations (3.12), (3.14) et (3.16), la condition limite à la paroi, considérée à l'altitude $z = z_p$, s'écrit pour un écoulement bidimensionnel :

$$ \mu_t \left. \frac{\partial \bar{u}}{\partial z} \right|_{z=z_p} = \frac{\rho \kappa C_{\mu}^{1/4} k_p^{1/2}}{\ln \left( \frac{z_p - d}{z_0} \right)} [\bar{u}_p - \bar{u}_{\infty}] \qquad (3.17) $$
avec $\bar{u}_{\infty}$, la valeur de la vitesse horizontale à la paroi ($\bar{u}_{\infty}=0$). Cette condition est donc une condition implicite de type Fourier généralisée. Elle introduit un couplage entre la dynamique et la turbulence car la contrainte de cisaillement est estimée par une combinaison de deux échelles de vitesses différentes. La première vitesse de frottement s'appuie sur la loi logarithmique, et la seconde vitesse de cisaillement est spécifiée en fonction de l'énergie cinétique moyenne du mouvement turbulent 4.

<figure><img src="image_34.png" /><figcaption>Figure 3.4 : Volume de contrôle en 2D proche pariétal. Définition de la distance $z_p$</figcaption></figure>

Pour la composante verticale de la vitesse moyenne, nous considérons qu'il y a adhérence. La condition limite est une condition de Dirichlet :

$$ \bar{w}|_{z=z_p} = 0 \qquad (3.18) $$

En ce qui concerne la distance à la paroi, la première couche de cellules est située dans la couche à flux constant, c'est-à-dire à une distance $z_p = 10 z_0$ (figure 3.4). L'écoulement est donc résolu sur ces premiers points de calcul.

Par ailleurs, nous soulignons que la condition limite en $Bi_{\phi}$ est discrétisée directement dans les équations. La condition de Fourier généralisée ne peut s'appliquer que pour des mailles prises en bord du domaine physique.

Enfin, nous avons généralisé la condition limite (équation (3.17)) pour des parois aérodynamiquement lisses. Pour cela, la rugosité est calculée par le biais de la relation suivante :

$$ z_0 = \frac{z_*}{E_0} = \frac{\nu}{E_0 C_{\mu}^{1/4} k^{1/2}} \quad \text{avec} \quad E_0 = 8.8 \qquad (3.19) $$

<footnote n="4">4. Ce jeu de condition limite est nommé à « deux échelles de vitesse »</footnote>
À l'intérieur des couverts végétaux, les relations logarithmiques ne sont pas valables, sauf à une distance suffisamment proche du sol terrestre. Les conditions limites sont alors définies de manière identique (équation (3.17)) en considérant toutefois que le sol est de faible rugosité dynamique $z_0$.

### 3.3.2 Loi de paroi thermique

Pour modéliser les échanges de chaleur à la paroi, il est possible de procéder de différentes façons. Ce choix se fait en fonction des données du problème. En effet, si la densité de flux de chaleur sensible est connue, il suffit d'imposer une condition de flux surfacique, soit :

$$ - \lambda_t \left. \frac{\partial \bar{T}}{\partial z} \right|_{z=z_p} = Q_0 \qquad (3.20) $$

avec $Q_0$ constante ou éventuellement variable au cours du temps.

En revanche, lorsque seule la température pariétale est connue, il est nécessaire d'expliciter les échanges en région proche de la paroi. Pour cela, chaque maille de calcul est affectée d'une température de surface, d'une rugosité thermique, ainsi que d'une hauteur de déplacement.

La loi de paroi thermique se construit en s'appuyant sur les mêmes considérations énumérées précédemment. En stratification neutre, les relations suivantes sont vérifiées avec l'hypothèse d'équilibre local :

- La température moyenne suit une loi logarithmique dans la couche de surface :

$$ \bar{T}(z) - \bar{T}_s = \frac{P_{rt} T_*}{\kappa} \ln \left( \frac{z-d}{z_{0h}} \right) \qquad (3.21) $$

Cette relation permet de déduire la température de frottement :

$$ T_* = \frac{\kappa}{P_{rt} \ln \left( \frac{z-d}{z_{0h}} \right)} \left[ \bar{T}(z) - \bar{T}_s \right] \qquad (3.22) $$

- Le flux de chaleur turbulent est constant et s'exprime en fonction des échelles caractéristiques $u_*$ et $T_*$ :

$$ \rho C_P \overline{T'w'} = -\rho C_P u_* T_* \qquad (3.23) $$

- Ce flux s'explicite par une modélisation phénoménologique :

$$ \rho C_P \overline{T'w'} = -\lambda_t \frac{\partial \bar{T}}{\partial z} \qquad (3.24) $$
En combinant les équations (3.16), (3.22) et (3.24), le flux de chaleur turbulent normal à la paroi s'écrit :

$$ - \lambda_t \left. \frac{\partial \bar{T}}{\partial z} \right|_{z=z_p} = - \frac{\rho C_p \kappa C_\mu^{1/4} k_p^{1/2}}{P_{rt} \ln \left( \frac{z_p - d}{z_{0h}} \right)} \left[ \bar{T}_p - \bar{T}_s \right] \quad (3.25) $$

C'est une condition limite implicite de type Fourier généralisée.

Enfin, nous rappelons que la rugosité thermique $z_{0h}$ est définie par :

$$ z_{0h} = z_0 \exp(-\kappa P(H^+)) \quad (3.26) $$

$P(H^+)$ est la résistance au transfert de chaleur dans le cas d'une paroi rugueuse. Son expression est donnée par l'équation (1.93) à la section 1.4.3.

### 3.3.3 Traitement des variables turbulentes

Il est également nécessaire de spécifier les conditions aux limites pour les variables turbulentes $k$ et $\epsilon$. La liste n'est pas exhaustive, chaque code de calcul trouve un jeu de conditions limites adéquat. Nous avons retenu des conditions limites de type « anglo-saxon » [Laurence 97].

Dans la zone logarithmique, l'énergie cinétique moyenne du mouvement turbulent est constante. La condition limite qui en découle est donc une condition de Neumann homogène :

$$ \left. \frac{\partial k}{\partial z} \right|_{z=z_p} = 0 \quad (3.27) $$

La valeur à la paroi est alors le résultat de l'équation d'équilibre (transport, production, diffusion et dissipation) en spécifiant que l'apport diffusif à travers la paroi est nul.

Pour le taux de dissipation, elle est définie par les hypothèses de compatibilité avec la longueur de mélange et d'équilibre :

- La production par cisaillement, pour un écoulement développé suivant $x$, s'écrit :

$$ P_k = - \overline{u'w'} \frac{\partial \bar{u}}{\partial z} = \frac{u_*^3}{\kappa (z-d)} \quad (3.28) $$

- L'hypothèse d'équilibre local conduit à :

$$ P_k = \epsilon = \frac{u_*^3}{\kappa (z-d)} \quad (3.29) $$
En exprimant dans l'équation (3.29) la vitesse de frottement en fonction de $k$ (relation 3.15), la condition limite pour $\varepsilon$ est donnée par :

$$ \varepsilon]_{z=z_p} = \frac{C_{\mu}^{3/4} k_p^{3/2}}{\kappa (z_p - d)} \qquad (3.30) $$

C'est une condition limite de Dirichlet. La valeur à la paroi du taux de dissipation est directement couplée à la valeur calculée de $k$, au même noeud.

Ces conditions limites peuvent être déconnectées de la loi de paroi dynamique. Par conséquent, elles peuvent s'appliquer dans tout le domaine de calcul et sont notamment utilisées pour la prise en compte d'inclusion solide dans le domaine. Enfin, soulignons que cette formulation « anglo-saxonne » s'applique aussi près des points d'arrêt, alors qu'imposer des conditions limites en fonction de $u_*$ (relation 3.29) aboutit à un niveau de turbulence nul, lorsque la vitesse s'annule.

L'utilisation de ces conditions limites (loi de paroi dynamique et jeu « anglo-saxon » pour $k$ et $\varepsilon$) ne reste toutefois justifiée que pour des écoulements pleinement développés. Cependant, ces formulations donnent des résultats assez satisfaisants en situation de non-équilibre, du fait du couplage entre la dynamique et le mouvement turbulent à travers $k$.

## 3.4 La couche de surface et les conditions limites qui en découlent

Les conditions aux limites de type Neumann homogène, appliquées à l'ensemble des variables dans la couche de surface, ne sont pas satisfaisantes aussi bien d'un point de vue physique que numérique [Costes 96]. En effet, en appliquant la condition limite de Neumann pour $\bar{u}$ au haut du domaine physique en $z = z_H$, il en découle :

$$ \left. \frac{\partial \bar{u}}{\partial z} \right|_{z=z_H} = 0 \qquad (3.31) $$

Cette relation est en totale contradiction avec une couche à flux constant, dont le gradient de $\bar{u}$ vérifie la relation suivante :

$$ \frac{\partial \bar{u}}{\partial z} = \frac{u_*}{\kappa (z - d)} \qquad (3.32) $$

Bien qu'il soit possible d'intercaler une couche de vitesse constante entre le bord du domaine et la couche à flux constant, il est toutefois plus satisfaisant d'adapter nos conditions limites aux phénomènes physiques modélisés. Nous avons donc redéfini des conditions au
bord supérieur du domaine physique, en s'appuyant sur les relations et les propriétés de la couche de surface.

Les conditions limites au plafond pour modéliser une couche à flux constant, s'écrivent pour les différentes variables :

- La vitesse est couplée à l'énergie cinétique moyenne du mouvement turbulent. La vitesse de frottement $u_*$ est spécifiée en fonction de $k$, suivant la relation (3.16) :

$$ \left. \frac{\partial \bar{u}}{\partial z} \right|_{z=z_H} = \frac{C_{\mu}^{1/4} k_H^{1/2}}{\kappa (z_H - d)} \quad (3.33) $$

Pour la composante verticale du vent, une condition de Neumann est appliquée.

- La condition limite de $k$ est une condition de Neumann homogène car elle est compatible avec un écoulement en équilibre (production=dissipation) :

$$ \left. \frac{\partial k}{\partial z} \right|_{z=z_H} = 0 \quad (3.34) $$

- Enfin, la dissipation est couplée à l'énergie cinétique $k$ :

$$ \varepsilon|_{z=z_H} = \frac{C_{\mu}^{3/4} k_H^{3/2}}{\kappa (z_H - d)} \quad (3.35) $$

Ces conditions limites doivent s'appliquer à une hauteur $z_H$ suffisante, compatible avec une couche à flux constant, non perturbée. Nous notons également que la condition limite en entrée du domaine pour $k$ est une donnée importante, puisque par son biais, il s'applique un forçage du cisaillement en haut du domaine.

## 3.5 Généralisations des conditions limites en stratification thermique

Lorsque l'atmosphère est en stratification stable ou instable, la turbulence est inhibée ou au contraire accentuée par ces effets thermiques. Les relations flux-gradients ne sont plus valables et doivent être modifiées. Afin de tenir compte de la stabilité thermique de l'atmosphère, nous avons corrigé les profils et les flux turbulents selon la théorie de similitude de Monin et Obukhov. De ces relations, il est possible de généraliser les lois de paroi et les conditions au plafond.
Les procédures pour expliciter les conditions limites sont identiques au cas adiabatique :

- À l'aide des profils moyens corrigés par les fonctions de similitude $\psi$ (équations (1.113) et (1.114) définies à la section 1.5.3) nous obtenons des conditions de type Fourier généralisées pour la vitesse et la température en région proche paroi :

$$ \mu_t \left. \frac{\partial \bar{u}}{\partial z} \right|_{z=z_p} = \frac{\rho \kappa C_\mu^{1/4} k_p^{1/2}}{\left[ \ln \left( \frac{z_p - d}{z_0} \right) - \psi_m \left( \frac{z_p - d}{L} \right) + \psi_m \left( \frac{z_0}{L} \right) \right]} \left[ \bar{u}_p - u_\infty \right] \quad (3.36) $$

$$ - \lambda_t \left. \frac{\partial \bar{T}}{\partial z} \right|_{z=z_p} = - \frac{\rho C_p \kappa C_\mu^{1/4} k_p^{1/2}}{P_{rt} \left[ \ln \left( \frac{z_p - d}{z_{0h}} \right) - \psi_h \left( \frac{z_p - d}{L} \right) + \psi_h \left( \frac{z_{0h}}{L} \right) \right]} \left[ \bar{T}_p - \bar{T}_s \right] \quad (3.37) $$

Pour l'énergie cinétique moyenne du mouvement turbulent, la condition limite est une condition de flux nul, puisque l'écoulement est en équilibre local. La production totale de l'énergie cinétique (production d'origine dynamique et d'origine thermique) est égale à son taux de dissipation. La condition limite pour le taux de dissipation est rectifiée, puisque la dissipation contrebalance la production d'origine mécanique et thermique :

$$ \varepsilon|_{z=z_p} = \frac{C_\mu^{3/4} k_p^{3/2}}{\kappa (z_p - d)} \phi_\varepsilon \left( \frac{z_p - d}{L} \right) \quad (3.38) $$

- Les conditions aux limites à la frontière supérieure du domaine s'établissent en corrigeant les flux turbulents à l'aide des fonctions de similitude $\phi$ (équations (1.111) et (1.112)) définies à la section 1.5.3. Les conditions aux limites généralisées sont :

$$ \left. \frac{\partial \bar{u}}{\partial z} \right|_{z=z_H} = \frac{C_\mu^{1/4} k_H^{1/2}}{\kappa (z_H - d)} \phi_m \left( \frac{z_H - d}{L} \right) \quad (3.39) $$

$$ \left. \frac{\partial \bar{T}}{\partial z} \right|_{z=z_H} = \frac{T_*}{\kappa (z_H - d)} \phi_h \left( \frac{z_H - d}{L} \right) \quad (3.40) $$

$$ \varepsilon|_{z=z_H} = \frac{C_\mu^{3/4} k_H^{3/2}}{\kappa (z_H - d)} \phi_\varepsilon \left( \frac{z_H - d}{L} \right) \quad (3.41) $$

Tous les processus turbulents ne dépendent que de l'indice de Monin et Obukhov. Il est donc nécessaire de déterminer la longueur et l'indice de Monin et Obukhov pour l'élabo-ration des conditions limites. Nous rappelons que la longueur $L$ et l'indice de stabilité $\xi$
sont donnés par :

$$L = \frac{u_*^2}{\kappa\beta T_*} \qquad \xi = \frac{z-d}{L}$$

Il existe plusieurs méthodes pour la détermination de ces grandeurs. Il est possible d'établir des relations analytiques entre $\xi$ et le nombre de Richardson volumique [Byun 90, Guilloteau 99]. Bien que coûteuse en temps de calcul, nous avons choisi une méthode itérative.

Cela consiste à évaluer la vitesse de frottement $u_*$, ainsi que la température de frottement $T_*$, en imposant une valeur initiale infinie pour $L$ ($L$ étant infinie pour un écoulement en stratification thermique neutre). Cette étape permet alors de calculer $u_*$ et $T_*$, à une hauteur proche du sol. Des valeurs obtenues, nous pouvons déduire l'échelle de longueur $L$ ainsi que $\xi$ à cette même hauteur, et élaborer ainsi les conditions aux limites. Ces étapes sont réitérées jusqu'à convergence.

## 3.6 Les conditions limites pour l'équation de transport des particules

Pour la dispersion des particules, nous avons défini deux échelles d'étude qui utilisent une équation de transport de type convection-diffusion. Cette équation de transport nécessite des conditions aux limites en haut et en bas du domaine physique. Les conditions limites appliquées au sol diffèrent selon l'échelle d'observation.

### Conditions limites au sol

À l'échelle du paysage, nous avons vu que les processus de dépôt dépendent de nombreux paramètres, tels que la turbulence, le couvert et les caractéristiques de la particule. Près du sol, le flux de particules déposées $F_d$ est proportionnel à la concentration, la constante de proportionnalité définissant la vitesse de dépôt. La condition limite s'écrit à la paroi en $z = z_p$ :

$$\left. \frac{v_t}{\sigma_c} \frac{\partial \bar{C}}{\partial z} + v_s \bar{C} \right|_{z=z_p} = F_d = v_d [\bar{C}_p - C_\infty] \qquad (3.42)$$

Pour la définition de la vitesse de dépôt, nous avons utilisé l'expression suggérée par [McCartney 91] :

$$v_d = 2v_s \qquad (3.43)$$
$v_d$ étant la vitesse de dépôt globale dans la sous-couche rugueuse et par conséquent, $z_p$ est définie comme une hauteur supérieure à 10 fois la rugosité (contrainte identique aux conditions limites pour la dynamique).

Dans le couvert, les processus de dépôt par impact et par sédimentation sur les éléments du couvert sont calculés explicitement. Dans ces conditions, seul le processus de dépôt par sédimentation est pris en compte près du sol. La condition limite qui en découle devient, puisque $v_d = v_s$ :

$$ \left. \frac{\nu_t}{\sigma_c} \frac{\partial \langle \overline{C} \rangle}{\partial z} + v_s \langle \overline{C} \rangle \right|_{z=z_p} = v_s \left[ \langle \overline{C} \rangle_p - C_\infty \right] \quad (3.44) $$

avec $C_\infty = 0$. Les conditions limites s'écrivent de manière identique, et cela en imposant une vitesse de dépôt différente suivant que nous résolvons explicitement l'écoulement dans le couvert ou non.

## Conditions limites supérieures

Dans l'atmosphère, la diffusion turbulente est sévèrement limitée au sommet de la couche limite atmosphérique car elle est coiffée par une couche d'entraînement, caractérisée par une stabilité thermique qui bloque les ascensions d'air. Il en résulte que peu de particules franchissent ce « mur ». En revanche, les particules peuvent franchir la couche de surface, dont la limite ne constitue pas une barrière physique nettement identifiable (et encore moins un « mur » pour les particules).

Pour prendre en compte ce flux de particules en haut du domaine physique, certains auteurs utilisent la notion de vitesse de fuite introduite par Maul (1975) :

$$ \frac{\nu_t}{\sigma_c} \frac{\partial \overline{C}}{\partial z} + v_s \overline{C} = -v_f \overline{C} \quad (3.45) $$

Toute la difficulté réside dans la détermination expérimentale de la vitesse de fuite $v_f$, ce qui explique que l'utilisation de cette condition limite s'est révélée peu satisfaisante.

Des études expérimentales ont montré que la concentration des particules tend toujours à décroître, à mesure que l'on s'éloigne du sol. La condition limite supérieure communément admise est alors de considérer que le flux net de particules, à une hauteur suffisamment éloignée du sol ou du couvert, est nul :

$$ \left. \frac{\nu_t}{\sigma_c} \frac{\partial \overline{C}}{\partial z} + v_s \overline{C} \right|_{z=z_H} = 0 \quad (3.46) $$
Étant donné les hauteurs auxquelles nous modélisons l'écoulement (couche de surface), cette condition limite est plus ou moins satisfaisante. Enfin, nous insistons sur le fait que lorsque les particules franchissent la couche de surface, les distances de dispersion peuvent s'accroître fortement, ce qui ne peut être pris en compte par notre modèle.

Bilan global de conservation de masse

<figure><img src="image_35.png" /><figcaption>Figure 3.5 : Flux convectifs de particules, intégrés entrant et sortant du domaine physique. Dépôt cumulé au sol</figcaption></figure>

Nous avons effectué un bilan global sur l'ensemble du domaine de calcul pour vérifier que les méthodes numériques sont bien conservatives. Pour cela, nous avons calculé les flux intégrés de particules convectées, entrant et sortant du domaine de calcul (figure 3.5).

Le flux de particules convectées intégré à travers le plan vertical, correspond au dé-
bit solide en suspension dans le fluide $[particule s^{-1}m^{-1}]$ (pour une largeur unitaire en
simulation bidimensionnelle) :

$$
\Phi_x(x) = \int_{z_p}^{z_H} \bar{u}(x, z) \bar{C}(x, z) \, dz \qquad (3.47)
$$

Près du sol, le flux déposé intégré suivant la direction horizontale sur toute la longueur du domaine rend compte du dépôt cumulé :

$$
\Phi_{cums} = \int_{0}^{L_x} v_d \bar{C}(x, z) \, dx \qquad (3.48)
$$

Dans le couvert, le flux cumulé de matière déposée sur la totalité de la végétation est évalué à l'aide de :

$$
\Phi_{cumh} = \int_{0}^{L_x} \int_{z_p}^{h} \langle \overline{S_p} \rangle (x, z) \, dx \, dz \quad (3.49)
$$
où $\langle \overline{S_p} \rangle$ [particule $m^{-3}s^{-1}$] représente le terme puits qui correspond aux dépôts des particules par unité de volume et de temps, à chaque tranche du couvert.

Le bilan global est respecté, lorsque le flux convectif intégré entrant, est égal au flux sortant et aux dépôts cumulés près du sol et dans la végétation :

$$ \Phi_x(0) = \Phi_x(L_x) + \Phi_{cums} + \Phi_{cumh} \quad (3.50) $$

Finalement, nous pouvons évaluer en pourcentage, l'erreur globale sur la conservation de la masse :

$$ E = 100 \left| 1 - \frac{\Phi_x(L_x) + \Phi_{cums} + \Phi_{cumh}}{\Phi_x(0)} \right| \quad (3.51) $$

Pour les simulations en 3 dimensions, il est nécessaire de considérer la troisième dimension introduite. L'expression (3.50) devient si $L_y$ est la largeur du domaine de calcul :

$$ \int_{0}^{L_y} \Phi_x(0) dy = \int_{0}^{L_y} \Phi_x(L) dy + \int_{0}^{L_y} \Phi_{cums} dy + \int_{0}^{L_y} \Phi_{cumh} dy \quad (3.52) $$

L'erreur globale est alors évaluée de manière similaire à l'équation (3.51).

Bien que la méthode des volumes finis soit par définition une méthode conservative, nous avons voulu vérifier que les modifications implantées restent cohérentes. Des bilans de masse ont été effectués pour différentes configurations de calcul et confirment que le modèle numérique est bien conservatif. Une perte de masse inférieure à 0.06% a été obtenue avec un schéma Hybride. Comme nous le verrons, ces flux sont très utiles, car ils permettent d'illustrer l'appauvrissement moyen d'un panache, ainsi que la répartition des dépôts.
<figure><img src="image_36.png" /></figure>
Conclusion

Ce chapitre présente la façon dont nous résolvons numériquement les équations qui régissent les phénomènes physiques. Le code de calcul numérique utilise la méthode des volumes finis et dispose de nombreuses techniques numériques.

Un concept original de maillage déconnecté permet simplement de résoudre
l'écoulement atmosphérique à différentes échelles d'observation. En se basant sur ce
concept, nous avons développé un modèle de canopée. Ce modèle a été complètement
intégré au code de calcul.

Des conditions limites propres aux écoulements atmosphériques, ont été discré-
tisées et implantées. Ces conditions limites permettent de traiter des configurations
de paysages hétérogènes, en stratification thermique.

Dans la partie qui suit, nous testerons la qualité et le comportement des modèles
physiques et des méthodes numériques.
<figure><img src="image_37.png" /></figure>
# Chapitre 4

## Validations des modèles physiques et numériques

### Introduction

L'objet de ce chapitre est de vérifier le bon comportement des modèles physiques et numériques. Pour cela, nous reproduirons une série de simulations numériques basées sur des expériences menées en soufflerie ou *in situ*, issues de la littérature. Ce chapitre présente successivement :

- Des simulations à l'échelle d'un paysage hétérogène. Cette première section regroupe les effets dynamiques et thermiques qui prennent naissance, au sein de l'écoulement, après une hétérogénéité simple de la surface;

- La seconde section est consacrée à la validation du modèle dynamique et du modèle de turbulence à l'échelle du couvert. Cette section est découpée en deux parties distinctes. Dans la première partie, nous nous sommes attachés à vérifier en détail, que les nombreuses échelles qui gouvernent la turbulence de couvert sont bien reproduites par notre modèle. La seconde partie porte sur l'étude de transitions complexes végétation-sol nu-végétation;

- Enfin, la dernière section traite de la dispersion à courte distance de spores *Lycopodium* à l'échelle d'un couvert végétal continu et hétérogène suivant la direction verticale.

Ces cas de validation du code de calcul ont été choisis, de telle sorte qu'ils permettent au mieux de s'assurer, que chaque étape prise séparément est correctement simulée.
<figure><img src="image_38.png" /></figure>
### 4.1 Les transitions dans la couche de surface

Les échanges convectifs près du sol sont influencés par les caractéristiques de la surface. Ils peuvent être modifiés lorsque l'écoulement rencontre un changement brutal de rugosité de surface ou de température pariétale. Afin d'illustrer nos propos et valider notre modèle physique et numérique, nous avons regroupé dans cette section deux configurations différentes d'écoulement sur une surface qui présente une hétérogénéité simple :

- Le premier cas test consiste à simuler le développement d'une couche limite interne dans la couche de surface, après un saut de rugosité dynamique, en stratification thermique neutre;

- Le deuxième cas test est l'étude d'une couche limite turbulente sur une plaque lisse qui présente un saut de température pariétale, en régime de convection forcée.

Nous avons choisi ces deux cas tests, car ils permettent de dissocier les effets de surface dynamiques des effets thermiques. Ces cas tests ont pour but de vérifier le comportement du modèle numérique de manière générale, et en particulier, de vérifier l'aptitude des lois de paroi dynamique et thermique à faire « ressentir » les effets de surface sur la couche limite atmosphérique.

#### 4.1.1 Transition lisse-rugueux

Dans la couche de surface, un saut de rugosité dynamique de la surface entraîne la formation de couches limites internes. Pour mettre en évidence ce phénomène, nous nous sommes appuyés sur l'expérience de [Bradley 68]. Cette simulation a pour but de vérifier que la loi de paroi rugueuse définie permet de traiter des configurations simples de saut de rugosité dynamique, et d'autre part, que les conditions aux limites au plafond du domaine reproduisent une couche de surface.

##### Description de l'expérience

L'expérience de [Bradley 68] a été menée en site naturel, sur une piste d'un aérodrome aménagée de pointes verticales pour créer un saut de rugosité. Les paramètres de surface relevés sont les suivants :

- À l'amont du saut de rugosité, la surface est lisse :

$$z_{01} = 2 \times 10^{-5} \text{ m et } d_1 = 0 \text{ m pour } x < 0 \text{ m}$$
– À l'aval du saut de rugosité, la surface est rugueuse :

$$z_{02} = 2.5 \times 10^{-3} m \text{ et } d_2 = 0 m \text{ pour } x \ge 0 m$$

Le paramètre de rugosité de cette expérience est $m = \frac{z_{02}}{z_{01}} = 125$.

## Les relations dans la couche de surface

Afin d'étudier l'influence d'un changement de rugosité sur l'écoulement, nous utilisons des grandeurs pertinentes dont l'évolution permet d'identifier le comportement de l'écoulement vis-à-vis du changement de rugosité.

Près du sol, une surface homogène provoque un forçage (cisaillement, flux turbulents, énergie cinétique, dissipation) sur l'écoulement. Cette zone est identifiée comme la couche de surface. Les flux turbulents sont constants au sein de cette couche et l'écoulement est en équilibre local avec la surface. En notant par $u_{*1}$ la vitesse de frottement au-dessus de cette surface, la situation d'équilibre permet de dégager les relations suivantes, dans le cas d'un écoulement en stratification thermique neutre :

$$\text{– } \phi_m \text{ est l'échelle de cisaillement sans dimension : } \phi_m = \frac{\kappa (z-d)}{u_{*1}} \frac{\partial \bar{u}}{\partial z} = 1$$

$$\text{– } \phi_k \text{ est l'énergie cinétique turbulente sans dimension : } \phi_k = \frac{C_{\mu}^{0.5} k}{u_{*1}^2} = 1$$

$$\text{– } \phi_{\varepsilon} \text{ est le taux de dissipation sans dimension : } \phi_{\varepsilon} = \frac{\kappa (z-d) \varepsilon}{u_{*1}^3} = 1$$

Donc, dans la couche de surface en équilibre, ces profils sans dimension ne varient pas, et en particulier, leur valeur est égale à 1.

## Paramètres de simulation

Les simulations, réalisées avec les modèles de turbulence $k-l$ et $k-\varepsilon$, sont menées en imposant en entrée du domaine de calcul un écoulement en équilibre avec la surface amont. La vitesse moyenne horizontale suit une loi logarithmique de paramètre de rugosité $z_{01}$ et $u_{*1} = 0.104 \text{ ms}^{-1}$. Les valeurs pour l'énergie cinétique moyenne du mouvement turbulent et du taux de dissipation, pour le modèle $k-\varepsilon$, sont déduites des relations présentées plus haut. Les conditions limites sont définies au sol, par la loi de paroi rugueuse et en haut du domaine, par les conditions de couche de surface.
Enfin, précisons qu'aucun calage numérique n'a été effectué au niveau des constantes des modèles de turbulence pour ces simulations. Par conséquent, les modèles de turbulence s'appuient sur les jeux des constantes classiques du modèle $k-l^1$ et du modèle standard $k-\varepsilon^2$.

### Résultats qualitatifs : les perturbations dans la couche de surface

<figure><img src="image_39.png" /><figcaption>Figure 4.1 : *Modifications de la structure de la couche de surface induites par le saut de rugosité défini par m = 125.* Résultats obtenus avec le modèle $k-\varepsilon$</figcaption></figure>

Les figures 4.1a, 4.1b, 4.2a et 4.2b représentent respectivement les profils sans dimension, de la vitesse moyenne horizontale, de l'échelle de cisaillement, de l'énergie cinétique moyenne du mouvement turbulent et du taux de dissipation depuis le saut de rugosité.

Dès que la couche limite atmosphérique atteint le saut de rugosité, la vitesse moyenne horizontale diminue brusquement, puis tend à se stabiliser à l'aval de la discontinuité (figure 4.1a) de manière progressive. Les profils de $\phi_m$, $\phi_k$ et $\phi_\varepsilon$ (figures 4.1b, 4.2a et 4.2b) permettent d'affirmer que l'écoulement est perturbé et complètement déséquilibré puisque leur valeur s'écarte de l'unité. Le développement de cette perturbation est associé à celui d'une couche limite interne $H(x)$ (figure 4.3) qui résulte du changement des

<footnote n="1">1. Tableau 1.1 à la section 1.3.3</footnote>

<footnote n="2">2. Tableau 1.2 à la section 1.3.4</footnote>
<figure><img src="image_40.png" /><figcaption>Figure 4.2 : Modifications de l'énergie cinétique moyenne du mouvement turbulent et de sa dissipation à l'aval du saut de rugosité m = 125. Modèle k-ε</figcaption></figure>

caractéristiques de la surface. Au-dessus de la couche interne, les profils $\phi_i$ retrouvent leur valeur amont : à cette hauteur, la couche limite atmosphérique n'est pas affectée par l'hétérogénéité de la surface et conserve sa structure amont.

Les courbes des paramètres $\phi_i$ dont la valeur est égale à 1 au-delà de la couche limite interne, confirment l'existence d'une couche à flux constant. La condition limite définie au plafond du domaine, simule correctement une couche de surface.

<figure><img src="image_41.png" /><figcaption>Figure 4.3 : Représentation schématique du développement d'une couche interne et d'une sous-couche interne lors d'un changement brutal de rugosité $z_{02} > z_{01}$</figcaption></figure>

De ces courbes, il est possible de déterminer différentes caractéristiques de l'écoulement : le cisaillement au sol, l'évolution de la couche interne et de la sous-couche interne.
## Le cisaillement au sol

Le calcul du cisaillement (sans dimension) au sol est obtenu à partir de la vitesse de cisaillement adimensionnelle, à la limite inférieure du domaine de calcul en $z_p$ :

$$ \phi_{sol}(x) = \phi_m(x, z_p) \qquad (4.1) $$

La courbe 4.4 représente l'échelle de cisaillement adimensionnée issue de notre simulation. Ces résultats numériques regroupent les résultats obtenus par le modèle de turbulence $k - \epsilon$ et le modèle $k - l$. Nous avons également reporté les mesures de Bradley pour valider la loi de paroi rugueuse et tester le comportement des deux modèles de turbulence vis-à-vis d'une hétérogénéité de surface.

<figure><img src="image_42.png" /><figcaption>Figure 4.4 : Évolution du cisaillement adimensionné au sol pour une transition de rugosité de surface lisse-rugueux [Bradley 68]. Le paramètre de rugosité vaut $m = 125$</figcaption></figure>

À l'amont du changement de rugosité, le cisaillement adimensionné est constant et sa valeur est égale à 1, car la couche limite est en équilibre. À la discontinuité, la contrainte de cisaillement adimensionnée augmente brusquement pour atteindre son maximum. Ce pic est bien reproduit par les deux modèles de turbulence. L'écoulement réagit donc au saut de rugosité, dont la présence entraîne un flux turbulent de quantité de mouvement important suivant la direction verticale au contact de la paroi la plus rugueuse, où l'écoulement est brusquement ralenti.
Pour de faibles fetchs (définis comme la distance de la discontinuité) le cisaillement diminue très rapidement vers une valeur asymptotique. Nous notons que la décroissance est mieux reproduite par le modèle $k - \varepsilon$. Nous pouvons donc conclure que l’écoulement se stabilise à une certaine distance du saut de rugosité, et ainsi s’ajuste aux nouvelles caractéristiques de la surface.

Le cisaillement adimensionné calculé est en bon accord avec les résultats de Bradley, notamment pour le modèle de turbulence $k-\varepsilon$. L’écoulement simulé est sensible aux effets dynamiques induits par l’hétérogénéité de la surface.

## Évolution de la couche limite interne et de la sous-couche interne en équilibre

En aval de la discontinuité, l’écoulement se réajuste à la nouvelle surface et le cisaillement tend vers une valeur limite constante. Cette mince couche définie près du sol, qui prend naissance au sein de la couche limite interne $H(x)$, est appelée sous-couche en équilibre $h(x)$. Elle tient son nom du fait que les flux turbulents sont conditionnés uniquement par les caractéristiques géométriques de la surface à l’aval du saut de rugosité.

Dans la couche interne, nous pouvons donc distinguer deux zones de propriétés différentes. La sous-couche interne en équilibre $h(x)$, qui est une couche à flux constant, suivie d'une zone de transition où l'écoulement est fortement imprégné des caractéristiques amont et aval de la surface. Cette zone assure la continuité des vitesses et des flux turbulents au-dessus de $H(x)$. Il s'agit du phénomène d'advection locale.

La couche limite atmosphérique résulte donc d'une succession de couches internes qui s'emboîtent. Lorsque la surface est homogène sur une échelle d'espace suffisante, et que l'écoulement n'est pas perturbé, la sous-couche en équilibre, qui se développe au sein de la couche interne, diffuse et forme alors une nouvelle couche à flux constant. Par conséquent, la structure d'un paysage (la taille d'une parcelle, la rugosité) influence l'écoulement non seulement à l'échelle locale, mais également à l'échelle régionale : on parle alors d'une rugosité régionale. En terme d'aménagement du territoire, la compréhension de l'évolution des couches internes suscite donc un grand intérêt.

Notre calcul de $H(x)$ est basé sur une estimation de la hauteur à laquelle la vitesse s'éloigne de 1 % de sa valeur amont non perturbée. $H(x)$ peut également se déduire du profil de l'échelle de cisaillement, de l'énergie cinétique pour des écarts de 1 % à 10 %. Il existe de nombreuses lois empiriques qui prédissent l'évolution asymptotique de la couche
limite interne. De façon générale, ces lois s’expriment comme une puissance du fetch :

$$ \frac{H}{z_{02}} = a \left( \frac{x}{z_{02}} \right)^b \qquad (4.2) $$

Le facteur de proportionnalité $a$ dépend de la définition utilisée pour le calcul de la couche limite interne. Sa valeur diffère fortement suivant les auteurs, pouvant varier de 0.25 à 1. En revanche, la puissance $b$, qui dépend de la stabilité thermique de l'atmosphère, est universelle. En condition adiabatique, elle vaut 0.8 [Elliott 92].

<figure><img src="image_43.png" /><figcaption>Figure 4.5 : Évolutions de la couche limite interne et de la sous-couche en équilibre –sans dimension– pour $m = 125$ [Bradley 68]. $H(x)$ est définie pour un écart de 1 % de la vitesse moyenne horizontale non perturbée et $h(x)$ est déduite pour un écart de 10 % du cisaillement au sol</figcaption></figure>

La définition de la sous-couche interne $h(x)$ est basée sur le cisaillement au sol. C'est la hauteur à laquelle le cisaillement s'écarte de 10 % de sa valeur au sol. Les deux couches étant homothétiques, l'évolution de la sous-couche en équilibre suit donc une loi puissance 0.8 du fetch, en condition adiabatique.
Sur la figure 4.5, nous avons tracé l'évolution de la couche limite interne et de la sous-couche interne en équilibre, adimensionnées par la longueur de rugosité à l'aval du saut, issue de notre simulation numérique. Nous y avons reporté les résultats expérimentaux du suivi de la couche limite interne de Bradley.

La couche limite interne fournie par le modèle de $k - \varepsilon$ est légèrement surestimée mais reste fidèle aux mesures de Bradley. Ce résultat est appréciable, puisque ce modèle de turbulence est connu pour ses bonnes prédicitions des écoulements cisaillés simples en équilibre ($P_k = \varepsilon$) mais reste peu adapté aux écoulements complexes marqués par une forte anisotropie et en situation de non-équilibre. Or, au saut de rugosité, l'écoulement a bien une direction privilégiée puisque l'hétérogénéité de la surface induit des transferts turbulents verticaux importants. De plus, l'écoulement est loin d'être en équilibre. Nous sommes donc satisfaits du comportement du modèle $k - \varepsilon$ et des résultats fournis. En revanche, le modèle $k - l$ suretime fortement la couche interne, mais nous notons toutefois, que la croissance en puissance 0.8 de la couche limite interne reste conforme aux observations expérimentales.

La sous-couche interne en équilibre (figure 4.5) n'est obtenue que pour un fetch supérieur à 1000, pour le modèle $k - \varepsilon$. En-dessous de cette distance, l'écoulement est trop marqué par la discontinuité de la surface et le retour à l'équilibre s'effectue de manière progressive. À l'aval de cette distance, l'écoulement s'adapte à la surface et la croissance de la sous-couche en équilibre obtenue suit bien une loi en puissance 0.8.

Nous n'avons pas obtenu de résultats cohérents sur l'évolution de la sous-couche interne en équilibre avec le modèle $k - l$. L'écoulement n'arrive donc pas à retrouver un état d'équilibre, et laisse supposer que la longueur de mélange $l_m = \kappa z$ n'est sûrement pas adaptée pour des sauts de rugosité. Toutes ces considérations confirment que le choix de la longueur de mélange est essentiel et que le modèle reste trop marqué par l'empirisme qui réside dans la formulation de $l_m$. Il est possible d'améliorer ce modèle en spécifiant une longueur de mélange propre à cette configuration d'écoulement, mais le manque d'universalité de $l_m$ est un handicap certain.

## Simulation tridimensionnelle

Nous avons étendu la modélisation de la couche limite de surface au cas tridimensionnel. Le calcul est entrepris avec le modèle $k - \varepsilon$, car il a donné de meilleurs résultats en 2 dimensions. Ce passage est nécessaire pour simuler des paysages constitués de mosaïques irrégulières, dont la position ne peut être traitée en 2 dimensions.
Afin de tester le comportement du code de calcul en simulation tridimensionnelle, nous avons repris l'expérience de [Bradley 68]. Il s'agit toutefois d'une simulation faux 3D car l'écoulement est bidimensionnel. Toutes les conditions aux limites sont identiques au cas bidimensionnel, à l'exception des conditions limites sur les faces latérales du domaine de calcul, qui sont des conditions de symétries.

<figure><img src="image_44.png" /><figcaption>Figure 4.6 : Évolution du cisaillement au sol pour un saut de rugosité dont le rapport est $m = 125$. Simulation 3D obtenue à l'aide du modèle $k - \varepsilon$.</figcaption></figure>

Sur la figure 4.6, nous avons représenté le cisaillement au sol (sans dimension) obtenu avec les paramètres de Bradley. La couche limite est sensible au saut de rugosité. L'écoulement dans sa totalité réagit de manière similaire au cas bidimensionnel. En particulier, nous pouvons apprécier la parfaite symétrie du code, ce qui reste important en 3D.

## Conclusion

Nous avons simulé une perturbation dynamique due à un saut de rugosité dans la couche de surface en 2 et en 3 dimensions. Bien que le modèle $k-l$ soit facile à mettre en
oeuvre en 2 dimensions et peu onéreux en temps de calcul, il ne donne pas de résultats
assez précis. Par conséquent, dans la suite de notre étude, nous utiliserons uniquement le
modèle $k - \varepsilon$ standard.

En vue des résultats numériques obtenus, le modèle $k - \varepsilon$ standard sans calage, ainsi que les conditions aux limites permettent :

- De modéliser la couche à flux constant ;

- De représenter le cisaillement au sol après un saut de rugosité ;

- De simuler le phénomène d'advection locale et de prédire l'évolution de la couche interne et de la sous-couche en équilibre.

Le modèle numérique défini est donc apte à reproduire les effets de surface sur les écoulements atmosphériques, en situation hétérogène simple, en utilisant le paramètre de surface $z_0$.

### 4.1.2 Saut de température pariétale

Nous avons considéré jusqu'ici les effets dynamiques induits par des surfaces hétérogènes. Les paysages présentent souvent une grande variabilité spatiale des températures de surface. En effet, le flux de chaleur sensible est souvent très différent d'une culture à l'autre, ou d'un couvert à l'autre, car l'activité biologique y est différente. Afin d'étudier l'influence d'un saut de température pariétale, nous nous sommes inspirés de l'expérience de [Charnay 79]. Ce cas test permet de valider la loi de paroi thermique, et le modèle physique en régime de convection forcée.

#### Description de l'expérience

Il s'agit de l'étude d'une couche limite turbulente qui se développe sur une plaque plane, sans gradient de pression. La plaque est lisse et présente un saut de température pariétale. L'écoulement d'air incident a une vitesse horizontale de $\bar{u}_0 = 10 \text{ ms}^{-1}$ et sa température, notée $T_0$, est constante. L'intensité turbulente résiduelle est de 0.3%. La plaque est d'une longueur de $2 \text{ m}$ de paramètres thermiques imposés :

- À l'amont du saut de température pariétale, la plaque est chauffée et sa température s'écarte de $23^{\circ}\text{C}$ par rapport à la température du fluide $T_0$ :

$$T_c = T_0 + 23^{\circ}\text{C} \quad \text{pour} \quad 0 \text{ m} \le x \le 0.7 \text{ m}$$
- La plaque est brusquement refroidie, et sa température devient inférieure de $3^\circ C$ par rapport à $T_0$:

$$T_f = T_0 - 3^\circ C \quad \text{pour} \quad 0.7 \, m < x \le 1.3 \, m$$

Cette configuration représente une transition chaud-froid. Les effets thermiques n'influencent pas l'écoulement dynamique, le régime est donc en convection forcée. Nous introduisons les grandeurs calculées pour ce cas test:

- La température adimensionnée :

$$\frac{T(z) - T_c}{T_0 - T_c}$$

- Le flux de chaleur turbulent sans dimension :

$$\frac{\overline{T'w'}}{\overline{u_0(T_0 - T_c)}}$$

### Paramètres de simulation

En entrée du domaine de calcul, la température $T_0$, la vitesse $\overline{u}_0 = 10 \, ms^{-1}$ et l'énergie cinétique du mouvement turbulent $k_0$ sont pris constants, le taux de dissipation est une fonction de $z$ et de $k_0$. Les conditions limites utilisées au sol sont la loi de paroi dynamique, généralisée au cas lisse et la loi de paroi thermique. Les conditions imposées en haut du domaine sont de type Dirichlet pour $\overline{u}$ (vitesse moyenne horizontale égale à la vitesse infinie non perturbée) et Neumann pour les autres variables.

### Développement de la couche limite dynamique

La première étape de cette simulation consiste à reproduire la couche limite dynamique turbulente obtenue expérimentalement par [Charnay 79]. L'évolution de l'épaisseur de la couche limite n'a été relevée expérimentalement qu'à partir du saut de température pariétale. Pour des conditions numériques et expérimentales identiques (plaque de 2 mètres) nous n'avons pu obtenir l'épaisseur de la couche limite expérimentale, représentée sur la figure 4.7. Cette différence, qui est notable, peut provenir d'origines diverses.

Tout d'abord, nous ne disposions pas d'informations complètes sur les grandeurs turbulentes en dehors de l'intensité résiduelle de turbulence. Les conditions d'entrée pour l'énergie cinétique et le taux de dissipation sont donc fixées de manière peu précise. De
plus, le modèle $k-\varepsilon$ reste peu adapté pour décrire l'établissement d'une couche limite turbulente sur une paroi lisse, dont les mécanismes de transition sont complexes (illustrés par la figure 4.8). Cette différence peut également provenir des conditions aux limites en haut du domaine qui ne sont pas très satisfaisantes car elles tendent à confiner l'écoulement.

<figure><img src="image_45.png" /><figcaption>Figure 4.7 : Évolution de l'épaisseur de la couche limite dynamique turbulente obtenue expérimentalement par [Charnay 79]. Le bord d'attaque est en $x = 0$ m. Couche limite obtenue numériquement en introduisant une plaque adiabatique de longueur 1.5 m à l'amont de la plaque chauffée</figcaption></figure>

Enfin, nous pensons que l'explication la plus plausible pour justifier une telle différence réside dans le dispositif expérimental. En effet, la couche limite turbulente obtenue par [Charnay 79] s'établit sur une distance inférieure à 0.7 mètre. Or, un écoulement d'air à $10 \, m \, s^{-1}$ sur une plaque plane lisse dont la viscosité $\nu$ est de l'ordre de $10^{-5}$ devient turbulent à une distance aval de 1 mètre [Lesieur 94]. Par conséquent, nous pouvons supposer que l'écoulement lors de l'expérience de [Charnay 79] a été nécessairement perturbé pour accélérer le mécanisme naturel de transition –par la présence d'éléments rugueux par exemple. Cette information n'est pas disponible, car la procédure expérimentale n'est pas détaillée.

Nous avons donc rajouté une plaque plane adiabatique de 1.5 mètres de longueur entre l'entrée du domaine de calcul et la plaque chauffée pour retrouver une épaisseur de la couche limite similaire à celle de [Charnay 79], au niveau du saut de température.
Nous obtenons alors une épaisseur de couche limite dynamique comparable à celle de [Charnay 79] au niveau du saut de température, exceptée à la dernière position (figure 4.7).

<figure><img src="image_46.png" /><figcaption>Figure 4.8 : Schéma de la transition en couche limite [Arbez 98]</figcaption></figure>

L'intérêt de cette démarche n'est pas de retrouver en intensité les résultats de [Charnay 79] mais d'étudier qualitativement les effets thermiques lors d'un saut de température. En effet, il n'est pas possible de réaliser la même expérience puisque l'écoulement est en régime de convection forcée. Dans ces conditions, les échanges thermiques sont fixés par le problème dynamique. Or, bien que la couche dynamique soit approximativement la même au saut de température, les échanges thermiques qui s'effectuent entre la position $x = 0 m$ et $x = 0.7 m$ ne peuvent être pris en compte, puisque les caractéristiques dynamiques sont différentes de celles de [Charnay 79] et que nous ne connaissons pas entre ces positions.

### Résultats qualitatifs : les profils de température et les flux de chaleur turbulents

Sur la figure 4.9a, nous avons tracé la température adimensionnée aux positions relevées par [Charnay 79]. Ces positions ont été prises à l'aval du saut de température pariétale. Il apparaît immédiatement que la couche limite thermique, obtenue par notre simulation, est fortement sous-estimée. Seule la température en proche paroi est calculée correctement. Compte tenu de toutes les remarques que nous avons formulées, il semble évident que la couche limite thermique reproduite numériquement ne se développe pas de manière similaire à celle de [Charnay 79], puisque le problème dynamique est différent.

Nous avons représenté les profils verticaux des flux de chaleur turbulents adimensionnés aux mêmes positions, à l'aval du saut de température sur la figure 4.9b. Nous pouvons
<figure><img src="image_47.png" /><figcaption>Figure 4.9 : Profils de la température et des flux de chaleur turbulents sans dimension, relevés à différentes distances (x en mètre) après le saut de température au sol [Charnay 79]</figcaption></figure>

observer une inversion de flux de chaleur pariétal. En effet, à l'amont du saut de température, la température au sol est supérieure à celle du fluide. Le flux de chaleur turbulent est positif : le fluide reçoit de la chaleur. La température au sol est brusquement refroidie et diminue de 26°C par rapport à la température initiale du fluide. Le fluide étant plus chaud que la plaque, à son contact, cède de la chaleur. L'échange est donc inversé et le flux de chaleur turbulent devient négatif. Ce phénomène est reproduit qualitativement par notre simulation numérique.

## Conclusion

Nous ne sommes pas en mesure de retrouver les résultats de [Charnay 79] car notre simulation numérique est différente. En effet, la procédure expérimentale n'étant pas suffisamment détaillée, il n'est pas possible de reproduire la même simulation : ce cas test ne nous permet pas de valider le code de calcul. De surcroît, le choix de l'expérience de [Charnay 79], comme cas de validation, est mauvais puisqu'il s'agit d'une couche limite sur une paroi aérodynamiquement lisse, qui ne met pas en jeu les lois de paroi dynamique et thermique en conditions représentatives de la couche de surface atmosphérique. Il est donc nécessaire de disposer d'un cas test plus adapté pour valider la loi de paroi thermique, et de manière générale, les effets de la stratification thermique de la basse atmosphère.
## 4.2 Simulations à l'échelle du couvert végétal

Nous abordons la partie validation du modèle physique à plus petite échelle, c'est-à-dire à l'échelle du couvert végétal. Les équations du mouvement et le modèle de turbulence calculent de manière explicite l'écoulement dans le couvert végétal. Nous avons choisi de présenter deux études différentes qui ont été réalisées dans la même soufflerie, en stratification thermique neutre :

- La première expérience consiste à étudier l'écoulement turbulent du vent dans un couvert végétal homogène. Toutes les grandeurs turbulentes (vent, flux turbulent de quantité de mouvement, énergie cinétique moyenne du mouvement turbulent, taux de dissipation) dans et au-dessus du couvert sont disponibles;

- La seconde expérience porte sur l'étude de l'écoulement dans le cas de transitions complexes végétation-clairière-végétation. Les effets aérodynamiques induits par la présence de la clairière sur l'écoulement sont explicités.

Ces deux expériences utilisent un couvert végétal constitué de tiges cylindriques en nylon, flexibles. À travers la première simulation, nous allons valider le modèle de canopée loin d'une discontinuité, de façon à mettre en valeur les caractéristiques de l'écoulement dans le couvert végétal, qui sont très différentes de la couche de surface. Enfin, la seconde simulation permet de tester le modèle de canopée en configuration complexe.

### 4.2.1 Écoulement dans un couvert végétal homogène

Pour résoudre l'écoulement turbulent à l'échelle d'un couvert végétal, nous avons développé un modèle de canopée à l'aide de l'opérateur moyenne spatio-temporelle. Ce modèle utilise une approche phénoménologique pour modéliser les flux turbulents de quantité de mouvement dans le couvert, relation qui est contestée. Cependant, dans la mesure où les équations d'évolution de l'énergie cinétique moyenne du mouvement turbulent et du taux de dissipation ont été modifiées pour tenir compte des variations d'échelle dans le couvert, comme il a été fait par [Liu 96], ou par [Wilson 99] pour un modèle de type $k-l$ avec succès, nous avons conservé cette méthode de fermeture pour les corrélations doubles. Afin de valider le modèle de canopée dans sa totalité, nous avons choisi de reproduire l'étude expérimentale menée par [Brunet 94], car elle dispose de données complètes sur l'ensemble des grandeurs turbulentes.
## Description de l'expérience

L'expérience a été réalisée dans une soufflerie de 0.65 m de hauteur, dont la paroi supérieure s'ajuste pour annuler le gradient de pression longitudinal. Les mesures des composantes des vecteurs vitesse ont été réalisées à l'aide d'un anémomètre à trois fils chauds. Cette instrumentation permet de mesurer les variations rapides de la vitesse et son emploi est plus adapté qu'un anémomètre classique, notamment à proximité du sommet du couvert où l'écoulement est caractérisé par une forte turbulence. L'écoulement à l'entrée de la soufflerie est constant, de vitesse moyenne horizontale de 10.2 ms⁻¹. À l'amont du couvert, le sol est recouvert de gravier d'une hauteur moyenne de 14 mm et permet ainsi le développement d'une couche limite turbulente. Le couvert végétal s'étend sur une distance horizontale x = 5.15 m, de caractéristiques physiques suivantes :

- Une hauteur de $h = 47 \text{ mm}$;

- Un coefficient de traînée moyen du couvert $C_d = 1.35$ et non des éléments constituant la végétation;

- Une densité de surface foliaire frontale $a = 10 \text{ m}^{-1}$.

Les profils verticaux de la vitesse moyenne horizontale et du flux turbulents de quantité de mouvement sont pris à 4.08 m du début du couvert. Les écarts types ont également été relevés ce qui nous permet d'accéder à l'énergie cinétique $k$. Enfin, le taux de dissipation a été déduit du spectre de $u$ et de $w$. Tous les résultats expérimentaux sont adimensionnés par le cisaillement au sommet du couvert dont la valeur a été déterminée à $u_* = 0.87 \text{ ms}^{-1}$.

## Paramètres de simulation

Le maillage utilisé, représenté sur la figure 4.10, est à pas variable pour les deux directions de l'espace. Suivant la direction verticale, le pas d'espace est pris constant dans le couvert, raffiné au niveau de sommet du couvert et exponentiel au-delà.

Pour les conditions limites, nous avons imposé en entrée du domaine, une vitesse moyenne horizontale uniforme de 10.2 ms⁻¹, l'énergie cinétique de la turbulence $k$ constante et le taux de dissipation $\varepsilon$ proportionnelle à $1/z$ (relations classiques de la couche de surface pour $k$ et $\varepsilon$ avec $u_* = 0.75 \text{ ms}^{-1}$). La condition limite au sol utilise la loi de paroi rugueuse, avec une faible rugosité dans le couvert, et les conditions au plafond sont des conditions de couche de surface.
<figure><img src="image_48.png" /><figcaption>Figure 4.10 : Exemple de maillage utilisé pour la résolution de l'écoulement à petite échelle. En sombre le couvert végétal</figcaption></figure>

Comme nous l'avons déjà évoqué, l'équation d'évolution du taux de dissipation est empreint d'un fort empirisme, car la dissipation de sillage de $\varepsilon$ est modélisée par le biais de la cascade énergétique de Kolmogorov, qui est accélérée dans le couvert végétal. Du reste, tous les termes clés de l'équation de $\varepsilon$ s'appuient sur des constantes très empiriques. Dans le cadre de ce travail, toutes les constantes du modèle de canopée $k - \varepsilon$ ont été conservées, en dehors de la constante $C_{D\varepsilon w}$ de la dissipation de sillage, car [Liu 96] a montré que le modèle de turbulence est particulièrement sensible à sa valeur et qu'un calage est nécessaire. Les résultats numériques obtenus ont donc été réalisés après une étude de sensibilité, présentée en seconde partie de cette validation, et qui a permis de fixer une valeur de $C_{D\varepsilon w} = 4.45$.

### Profils verticaux de $\langle \bar{u} \rangle$ et du flux turbulent $\langle \overline{u'w'} \rangle$

Nous avons tracé le profil de la vitesse moyenne horizontale obtenu par le modèle numérique sur la figure 4.11 en y reportant les mesures expérimentales de Brunet. Ce profil de la vitesse moyenne a l'allure typique des profils observés dans les couverts végétaux, dont la particularité est de s'atténuer de manière exponentielle, au fur et à mesure qu'on pénètre dans le couvert, et de posséder de très forts gradients à proximité du sommet. Ces gradients de vitesse intenses sont accompagnés par un point d'inflexion. C'est précisément ce point d'inflexion, avec des vitesses très différentes de part et d'autre du couvert, qui permet l'analogie avec la couche de mélange.

Dans la sous-couche rugueuse, qui s'étend sur une distance égale à 2 fois la hauteur du couvert $z_* = 2h$, l'écoulement turbulent est directement influencé par sa structure. Au-
<figure><img src="image_49.png" /><figcaption>Figure 4.11 : Profil de la vitesse moyenne horizontale dans et au-dessus du couvert végétal pris à la position x = 4.08 m [Brunet 94]</figcaption></figure>

delà, la vitesse suit une loi logarithmique classique, caractéristique de l'écoulement dans
la sous-couche inertielle. À cette échelle, le couvert est ressenti par l'écoulement comme
une paroi rugueuse de paramètres de rugosité z₀ et de hauteur de déplacement d.

Le profil de la vitesse moyenne horizontale calculé est en bon accord avec les résultats
expérimentaux de [Brunet 94] dans la canopée végétale et au-dessus. En revanche, nos
résultats s'écartent des mesures expérimentales loin du sommet du couvert. Cet écart est
dû à la paroi ajustable de la soufflerie, que nous ne pouvons reproduire numériquement.

Sur la figure 4.12, nous pouvons visualiser le flux turbulent de quantité de mouvement
$\langle u'w' \rangle$. Le flux turbulent est bien calculé dans les différentes régions de l'écoulement, avec
de forts transferts turbulents au sommet du couvert. Pour $z$ inférieur à $h$, les éléments du
couvert végétal absorbent par pression et viscosité le flux turbulent qui diminue (en valeur
absolue). Juste au-dessus du sommet du couvert, le flux turbulent observé expérimentale-
ment est maximum au sommet, et reste quasi-constant sur une distance qui s'étend jusqu'à
$2h$, et au-delà pour l'expérience de référence ($3h$). C'est la signature d'une sous-couche
en équilibre, ou à flux constant, présente au sein de la sous-couche rugueuse. Le flux
calculé par notre modèle présente bien un maximum à l'interface couvert-air cependant
légèrement surestimé (en valeur absolue). Cependant, nous n'observons pas clairement la
<figure><img src="image_50.png" /><figcaption>Figure 4.12 : *Flux turbulent de quantité de mouvement dans et au-dessus du couvert végétal pris à la position x = 4.08 m* [Brunet 94]</figcaption></figure>

présence d'une sous-couche en équilibre. Enfin, au-dessus de la sous-couche rugueuse, le flux turbulent décroît en s'éloignant du couvert, et reste peu influencé par la présence du couvert. Cette décroissance est souvent rencontrée dans les expériences de soufflerie.

Les résultats dans leur ensemble sont satisfaisants. Bien que légèrement surestimés, le flux turbulent de quantité de mouvement est bien reproduit, malgré l'utilisation d'une relation flux-gradient. Nous avons toutefois voulu vérifier la cohérence du modèle de turbulence en visualisant l'énergie cinétique moyenne du mouvement turbulent et son taux de dissipation.

### Profils de l'énergie cinétique et du taux de dissipation

Sur la figure 4.13, nous avons comparé l'énergie cinétique moyenne du mouvement turbulent simulée aux valeurs issues de [Brunet 94]. La variation de l'énergie cinétique dans les différentes parties de l'écoulement est bien retranscrite par le modèle de canopée. En particulier, nous notons que dans la sous-couche rugueuse, la structure turbulente ne dépend que des caractéristiques du couvert. L'énergie cinétique moyenne du mouvement turbulent est maximale au-dessus du couvert, où la turbulence est produite de manière
<figure><img src="image_51.png" /><figcaption>Figure 4.13 : Énergie cinétique moyenne du mouvement turbulent –sans dimension– dans et au-dessus du couvert végétal [Brunet 94]</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>Figure 4.14 : Taux de dissipation –sans dimension– dans et au-dessus du couvert végétal calculée à l'aide du spectre d'énergie [Brunet 94]</figcaption></figure>
prépondérante. Cette turbulence est induite par des structures cohérentes de dimension de l'ordre de la hauteur du couvert et qui sont générées par le cisaillement au sommet du couvert, ce qui est bien pris en compte par notre modèle de canopée.

Enfin, nous avons comparé le taux de dissipation avec les résultats de [Brunet 94] issus du spectre d'énergie. Bien que le profil exact du taux de dissipation soit l'objet de nombreuses controverses, nous pouvons apprécier nos résultats sur la figure 4.14. La dissipation présente un pic au sommet du couvert, et rend compte de l'activité turbulente importante qui s'effectue à l'échelle de $h$, par les structures cohérentes. Pour $z$ inférieur à $h$, l'extra-dissipation qui résulte de la dissipation de sillage des éléments du couvert, est bien prise en compte.

### Influence de la constante $C_{D\epsilon_w}$ du modèle de turbulence de canopée

La constante $C_{D\epsilon_w}$ qui apparaît dans le terme de dissipation de sillage $D_{\epsilon_w}$ de $\epsilon$ joue un rôle crucial dans le bilan de $\epsilon$ et par conséquent, sur l'énergie cinétique moyenne du mouvement turbulent. Nous avons donc mené une étude de sensibilité concernant l'influence de $C_{D\epsilon_w}$ sur le comportement du modèle turbulent de canopée, les autres constantes étant conservées.

Les profils de $k$ et de $\epsilon$ sont comparés sur la figure 4.15 aux résultats que nous venons de valider, et servent de référence :

- Pour la constante donnée par la relation de Kolmogorov $C_{D\epsilon_w} = 6$, nous pouvons observer sur la figure 4.15a le profil vertical du taux de dissipation. Au sol (premier noeud du maillage) la dissipation est trop importante et est à l'origine du cisaillement au sol, forcé par la loi de paroi. Au-dessus et dans toute la partie basse du couvert ($z < h/2$) l'évolution du taux de dissipation est complètement fausse, puisque la dissipation tend à diminuer, à mesure que l'on s'éloigne du sol et que l'on pénètre dans le couvert. De surcroît, la dissipation est trop faiblement estimée : dans cette région, l'activité visqueuse générée par les éléments du couvert est sous-estimée car mal décrite.

Dans la partie supérieure du couvert ($z > h/2$) le taux de dissipation croît en se rapprochant du sommet pour atteindre des valeurs maximales au-dessus. Bien que ces valeurs soient de trop faible intensité, l'écoulement turbulent à proximité du sommet ressent le couvert, puisque l'évolution du taux de dissipation est correcte. L'extra-dissipation générée par les éléments du couvert est reproduite par le modèle,
mais a lieu essentiellement à hauteur de $h$ et demeure trop insignifiante.

De cette analyse, nous déduisons que la relation de Kolmogorov n’est pas apte à retranscrire la contribution des nombreuses échelles qui gouvernent l’écoulement, puisqu’elle entraîne une extra-dissipation de l’ordre de $h$. L’énergie cinétique du mouvement turbulent représentée sur la figure 4.15b est trop importante et corrobore le fait que l’extra-dissipation est trop faiblement estimée. Toutes ces considérations nous ont incités à réduire le terme puits de l’équation de $\varepsilon$, c’est-à-dire la dissipation de sillage à travers $C_{D\varepsilon w}$, pour tenir compte du second processus de destruction de l’énergie cinétique, imposée par les éléments du couvert végétal;

<figure><img src="image_53.png" /><figcaption>Figure 4.15 : Étude de sensibilité du coefficient $C_{D\varepsilon w}$ du terme de dissipation de sillage $D_{\varepsilon w}$ du taux de dissipation. Les valeurs choisies sont données par la relation de Kolmogorov $C_{D\varepsilon w} = 6$ et par [Liu 96] $C_{D\varepsilon w} = 2.4$</figcaption></figure>

- En se basant sur les travaux numériques de [Liu 96], nous avons diminué la constante à $C_{D\varepsilon w} = 2.4$, valeur qu’il préconise d’utiliser. Les résultats obtenus par notre modèle numérique sont reportés sur la figure 4.15. Dans le couvert, les profils de $k$ et de $\varepsilon$ (figures 4.15b et 4.15a) retrouvent une allure comparable aux profils de référence, mais restent très faibles en intensité. Ce cas est en fait différent de celui que nous venons de voir.
En effet, la production de sillage et l'extra-dissipation sont bien prises en compte dans le couvert, mais demeurent complètement inhibées par le terme de production par cisaillement et sa dissipation associée, dont l'activité est prépondérante au sommet du couvert. La dissipation est concentrée au sommet du couvert où elle atteint un pic important, et diminue en s'éloignant de $z = h$. L'énergie cinétique atteint une valeur maximale en $h$ et varie très faiblement au-delà : la turbulence est donc en équilibre $P_k = \varepsilon$. Nous pouvons déduire que c'est principalement la production et la dissipation par cisaillement qui gouverne l'équation d'évolution de $k$ et de $\varepsilon$. Or, ceci est en totale contradiction avec la turbulence de couvert, pour laquelle la production par cisaillement est faible devant la production par sillage dans le couvert, qui ne retrouve son rôle majeur que pour des distances supérieures à $2h$.

Contrairement aux modèles $k-l$ où il est nécessaire d'expliciter une longueur de mélange appropriée pour traiter la turbulence de couvert, le modèle $k-\varepsilon$ permet de prendre en compte les différentes échelles, ainsi que leurs contributions. Cette prise en compte passe par une diminution de la dissipation de sillage de $\varepsilon$. Toutefois, le choix de cette constante est délicat, car elle accentue ou inhibe les phénomènes physiques, propres à la turbulence de couvert.

## Conclusion

Nous avons développé un modèle de turbulence de type $k-\varepsilon$ pour traiter l'écoulement à l'échelle du couvert végétal. Les résultats obtenus sont satisfaisants, le modèle défini reproductuant la turbulence de couvert. Cette turbulence est caractérisée par de nombreuses échelles très différentes : d'une part, une turbulence active associée à des structures cohérentes de taille de la hauteur du couvert. D'autre part, la présence des éléments du couvert végétal fournit une activité visqueuse à des échelles comparables et supérieures aux éléments. Pour tenir compte de cette activité, nous avons diminué le terme puits du taux de dissipation, en modifiant la constante issue de la cascade énergétique de Kolmogorov, comme il a été recommandé par [Liu 96]. Le modèle numérique a été calé sur l'expérience de [Brunet 94], dont les mesures ont été effectuées loin de la discontinuité. Nous avons retenu une valeur de $C_{D\varepsilon w} = 4.45$, que nous conserverons pour toutes les simulations de ce mémoire. Enfin, il apparaît clairement que cette constante est loin d'être universelle, puisque nous avons utilisé une valeur très différente de celle de [Liu 96] –qui pourtant a obtenu des résultats très satisfaisants. Compte tenu de l'empirisme qui réside dans les
termes de l'équation d'évolution de $k$ et de $\varepsilon$, nous envisageons par la suite de mener une étude de sensibilité plus complète, sur l'ensemble des constantes du modèle de turbulence, en particulier les constantes intervenant dans les termes de production de sillage.

<figure><img src="image_54.png" /><figcaption>Figure 4.16 : Vecteurs vitesse à la transition sol nu-végétation positionnée en $x = 0$ mètre</figcaption></figure>

Pour définir et valider notre modèle de canopée, nous nous sommes placés loin de la discontinuité sol nu-végétation, alors que le calcul a été réalisé sur l'ensemble du dispositif expérimental (figure 4.16). Il est donc nécessaire de confirmer le bon comportement du modèle à proximité de la transition. La section suivante porte précisément sur l'étude d'un écoulement plus complexe, faisant intervenir deux transitions végétation-clairière-végétation.

## 4.2.2 Écoulement dans des transitions végétation-clairière-végétation

Il est maintenant reconnu qu'une discontinuité de type clairière dans une couverture végétale a des conséquences sur le microclimat local. L'écoulement, ainsi que la turbulence sont fortement perturbés par la présence de la discontinuité qui induit des effets aérodynamiques.

De nombreuses études s'intéressent à ces écoulements qui s'avèrent complexes et assez mal connus. L'observation directe sur certains sites naturels a permis de prendre conscience des manifestations occasionnées par ces discontinuités. Nous pouvons citer pour exemple, la délocalisation de nids d'oiseaux constatée, suite à l'aménagement de clairière artificielle dans un site [Raupach 87]. De manière plus spectaculaire, lors de la tempête en décembre 1999, d'importants dégâts sont survenus dans les forêts situées à proximité d'espaces ouverts.
À travers cette étude numérique, qui s'appuie sur l'expérience de [Raupach 87] menée en soufflerie, nous allons donc tenter d'entrevoir les manifestations de ces perturbations et valider le modèle de canopée pour ces transitions complexes.

### Description de l'expérience

L'écoulement dans la soufflerie permet de représenter le vent dans un couvert, les paramètres physiques du modèle de canopée étant :

- Une hauteur de $h = 47 \text{ mm}$;

- Un coefficient de traînée moyen $C_d = 0.9$;

- Une densité de surface foliaire frontale $a = 10 \text{ m}^{-1}$.

Il s'agit de la même maquette que dans le cas précédent. La différence de $C_d$ provient du fait que pour cette expérience, la vitesse était plus élevée et la maquette était en bonne similitude avec une forêt [Raupach 87]. Enfin, une clairière de dimension variable représentée par un sol rugueux $z_0 = 0.15 \text{ mm}$, a été aménagée dans la soufflerie pour caractériser la discontinuité.

Deux configurations expérimentales ont été étudiées : une clairière de 1 mètre de longueur (soit environ $21.3 \text{ h}$) et une clairière de 0.2 mètre de longueur ($4.3 \text{ h}$). Les profils verticaux de la vitesse moyenne horizontale et les écarts types ont été mesurés à l'amont de la discontinuité, située en $x = 0$, dans la clairière et enfin à l'aval de la clairière, c'est-à-dire dans le couvert (figure 4.17). Ce dispositif représente donc une transition forêt-clairière ou rugueux-lisse et clairière-forêt ou lisse-rugueux.

### Paramètres de simulation

Les conditions aux limites sont identiques à la validation que nous avons effectuée précédemment, avec un vent uniforme horizontal de $11 \text{ m s}^{-1}$.

La clairière de longueur notée $L_c$ est positionnée en $x = 0 \text{ m}$. Afin de dissocier les différentes parties de l'écoulement suivant la direction horizontale, les positions relatives à la transition végétation-clairière, sont désignées par l'abscisse $x$ et les positions relatives à la transition clairière-végétation sont désignées par la variable $X = x - L_c$.
## Résultats quantitatifs : les vecteurs vitesse

La figure 4.17 représente les vecteurs vitesse pour ce dispositif discontinu. Nous pouvons observer la réponse de l’écoulement vis-à-vis du changement de nature du milieu. Dans le couvert, la vitesse moyenne suit une variation exponentielle. Dans la clairière, l’écoulement s’accélère et s’ajuste progressivement à la nouvelle surface. Enfin, à l’aval de la clairière, la vitesse moyenne diminue à nouveau lorsque l’écoulement pénètre dans le couvert, et retrouve progressivement sa variation exponentielle.

<figure><img src="image_55.png" /><figcaption>Figure 4.17 : Vecteurs vitesse pour une clairière de 1 mètre de longueur ou $L_c = 21.3 \text{ h}$</figcaption></figure>

Nous pouvons également constater que, contrairement aux écoulements autour d'obstacles solides, où il s'observe des tourbillons à l'amont et à l'aval de l'obstacle, ici il ne s'établit pas de recirculation. En réalité, des visualisations à l'aide de traceur, comme de la fumée, ont permis de déceler la présence de tourbillons très intermittents et non persistants, à l'amont et à l'aval de certains couverts, que des mesures directes ne permettent pas de caractériser. Du reste, [Raupach 87] a confirmé lors de cette expérience, la présence de ces tourbillons.

### Vitesse moyenne horizontale et énergie cinétique pour $L_c = 1 \text{ m} = 21.3 \text{ h}$

Sur les figures 4.18 et 4.19, nous avons représenté respectivement les champs simulations de la vitesse moyenne horizontale et de l'énergie cinétique moyenne du mouvement turbulent. Sur la figure 4.20, nous avons comparé les profils de la vitesse moyenne horizontale aux mesures expérimentales de [Raupach 87] effectuées à diverses positions du dispositif. Les profils de l'énergie cinétique moyenne comparés aux mesures sont tracés sur
la figure 4.21. Les phénomènes physiques étant différents pour la transition végétation-clairière (rugueux-lisse) et clairière-végétation (lisse-rugueux), nous dissocions ces cas.

<figure><img src="image_56.png" /><figcaption>Figure 4.18 : Vitesse moyenne horizontale pour le dispositif $L_c = 1 \text{ m} = 21.3 \text{ h}$</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>Figure 4.19 : Énergie cinétique pour le dispositif $L_c = 1 \text{ m} = 21.3 \text{ h}$</figcaption></figure>

### Transition végétation-clairière

Cette transition végétation-clairière correspond aux profils relevés aux positions $x = -0.1 \text{ m} = -2.1 \text{ h}$, $x = 0.1 \text{ m} = 2.1 \text{ h}$, $x = 0.2 \text{ m} = 4.3 \text{ h}$ et $x = 0.6 \text{ m} = 12.8 \text{ h}$. En sortie du couvert, lorsque l'écoulement rencontre la clairière, il s'ajuste à cette surface de plusieurs façons :

- La clairière perturbe l'écoulement (figure 4.18). Cette perturbation est due à la disparition de la force de traînée imposée par le couvert. Nos résultats numériques, aux
<figure><img src="image_58.png" /><figcaption>Figure 4.20 : Profils de la vitesse moyenne horizontale à différentes positions (x et X = x - Lc en mètres). Clairière de dimension Lc = 21.3 h [Raupach 87]</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>Figure 4.21 : Énergie cinétique moyenne du mouvement turbulent à différentes positions (x et $X = x - L_c$ en mètres). Clairière de dimension $L_c = 21.3$ h
[Raupach 87]</figcaption></figure>
différentes positions à l'aval du couvert, sont comparés au profil initial de la vitesse moyenne horizontale défini à la position $x = -2.1 \, h$. Nous observons l'accroissement progressif de la vitesse moyenne de l'écoulement le long de la clairière, dont le profil tend vers une loi logarithmique. La vitesse moyenne horizontale de l'écoulement est bien calculée par le modèle numérique, à toutes les positions ;

- En-dessous du niveau du couvert, l'écoulement n'est plus amorti par la présence de la canopée. L'énergie cinétique qui était à l'origine faible à cette hauteur, augmente en réponse du cisaillement local, généré sur la surface rugueuse de la clairière. Cette augmentation très importante en début de clairière, s'explique également par une diffusion vers le bas de la turbulence intense qui existait dans la couche de cisaillement, à hauteur du couvert (figure 4.19). Cette diffusion entraîne un phénomène de « overshoot» de $k$, avec de très fortes valeurs aux positions $x = 2.1 \, h$ et $x = 4.3 \, h$ et qui tend à diminuer lentement au-delà. Ce phénomène de sur-dépassement est bien observé sur nos résultats numériques, avec des valeurs légèrement surestimées. Ceci est sûrement dû au caractère trop diffusif du modèle $k-\varepsilon$, mais les résultats restent remarquables ;

- L'énergie cinétique, qui est maximale à hauteur du couvert, n'est plus entretenue par le fort cisaillement généré au sommet (figure 4.19). Cette variation d'échelle, qui dans le couvert est de l'ordre de $h$, et sur la surface rugueuse de l'ordre de la rugosité, entraîne une décroissance de l'énergie cinétique le long de la clairière à cette hauteur. Nous notons que cette décroissance est correctement simulée, surtout aux dernières positions de la clairière ;

- En bout de clairière, l'allure de $k$ permet d'affirmer que l'écoulement n'est toujours pas en équilibre avec la surface rugueuse de la clairière. Le retour à l'équilibre après une transition rugueux-lisse, requiert que l'énergie cinétique soit dissipée. Ce processus de dissipation est un processus lent et implique de grandes échelles de temps et d'espace pour le retour à l'équilibre. Une conséquence directe de ces transitions est l'apparition en aval d'une zone de sillage définie entre 3 $h$ et 15 $h$, dans laquelle l'intensité turbulente est forte à proximité du sol.

Les résultats numériques du vent et de la turbulence sont dans l'ensemble très satisfaisants pour la transition végétation-clairière.
### Transition clairière-végétation

Les positions qui correspondent à cette transition sont situées en $X = 0 m$, $X = 0.1 m = 2.1 h$, $X = 0.2 m = 4.3 h$, $X = 0.3 m = 6.4 h$ et $X = 0.5 m = 10.6 h$. Les profils de la vitesse et de l'énergie cinétique sont respectivement représentés sur les figures 4.20 et 4.21. Dès que l'écoulement rentre à nouveau dans le couvert :

- Pour $z < h$, la vitesse moyenne diminue fortement en raison de la résistance que le couvert oppose à l'écoulement. En même temps, au-dessus du couvert, l'écoulement est perturbé et accéléré localement en réponse d'un flux turbulent de quantité de mouvement vertical provoqué au contact direct de la canopée. Cette différence de vitesse va être à l'origine d'un point d'inflexion, dont l'évolution permet de retrouver rapidement la variation exponentielle de la vitesse moyenne. L'écoulement simulé est en bon accord avec les résultats expérimentaux à toutes les positions, excepté en $X = 10.6 h$, où la vitesse moyenne n'a pas retrouvé sa forme initiale. L'écoulement dans notre simulation n'est pas totalement développé à cette position et exige une distance supérieure, de l'ordre de $X = 21 h$, pour s'équilibrer à nouveau avec le couvert;

- Dans la partie basse du couvert, les écarts types du vent décroissent brusquement à l'aval de la clairière en $X = 2.1 h$, pour y atteindre de très faibles valeurs en $X = 4.3 h$ et $X = 6.4 h$. Ces valeurs sont remarquablement inférieures aux valeurs d'équilibre de couvert. L'énergie cinétique simulée est quant à elle surestimée, notamment à la position $X = 2.1 h$, puis diminue progressivement. Les profils de $k$ sont toutefois corrects aux dernières positions. Par similitude à la transition rugueux-lisse, un phénomène de « undershoot » brutal des variables turbulentes par rapport à leurs valeurs d'équilibre de couvert, a lieu pour les transitions lisse-rugueux dans la partie basse de la canopée. Cependant, les deux processus sont différents, et le mécanisme physique qui explique cette brusque diminution, implique un flux turbulent de quantité de mouvement important vers la surface la plus rugueuse. Par suite, l'activité turbulente est particulièrement violente au sommet du couvert où l'énergie cinétique se crée localement et rapidement (figure 4.19);

- C'est l'analyse du moment d'ordre 3 normalisé pour $\langle \bar{u} \rangle$ (Skewness<fnref n="3" />) non accessible par notre modèle de turbulence, qui permet d'affirmer que la turbulence à hauteur

<footnote n="3">3. $S_{k_u} = \langle \overline{u'^3} \rangle / \sigma_u^3$ dont la valeur est nulle pour une turbulence isotrope</footnote>
du couvert est très intense et intermittente dans la région comprise entre $X = 4.3 h$ et $X = 6.4 h$. Dès que l'écoulement rencontre le couvert, le fort cisaillement généré en $z = 2 h$ est à l'origine de la création de rafales (avec une probabilité d'occurrence plus grande que dans un couvert homogène) qui se développent le long du couvert et se retrouvent à une distance centrée en $X = 5 h$ à hauteur du couvert. Cette région qui s'étend entre $X = 3 h$ et $X = 7 h$, est le lieu d'une turbulence très intense, accompagnée de rafales particulièrement violentes (plus violentes que pour une turbulence de couvert à l'équilibre) car caractérisées par des vitesses instantanées trois fois plus élevées que la vitesse moyenne du vent [Raupach 87]. Cette région coïncide avec la zone où des dégâts occasionnés par le vent s'observent fréquemment.

- Enfin, en $X = 10.6 h$, les variables turbulentes mesurées expérimentalement augmentent à nouveau pour retrouver leurs valeurs d'équilibre. Cette augmentation s'explique par la trajectoire plongeante des rafales, qui pénètrent dans le couvert à des distances caractéristiques de l'ordre de $X = 5 h$ à $X = 10 h$, pour se dissiper et permettre de maintenir la turbulence profonde de couvert -la turbulence de couvert est principalement transportée et non créée localement. Ce retour à l'équilibre est donc rapide, car il est déterminé par le comportement du mouvement turbulent sur la surface la plus rugueuse, où nous avons vu que l'énergie cinétique moyenne du mouvement turbulent est créée rapidement. Enfin, le profil calculé de $k$ en derrière position $X = 10.6 h$ confirme que l'écoulement simulé n'est pas à l'équilibre et que ce retour à l'équilibre dans notre simulation exige des distances supérieures à $X = 10.6 h$. En fait, il s'agit du phénomène de « overshoot » que nous retrouvons à cette position puisque la valeur de $k$ est bien inférieure à sa valeur d'équilibre et augmente au-delà.

La transition clairière-végétation et le retour à l'équilibre évoque une turbulence violente et très intermittente, à l'aval de la clairière. Dans cette zone, l'apparition d'intenses rafales est accrue, ce qui lui vaut le nom de « enhanced gust zone » [Raupach 87]. Cette turbulence, qui est loin d'être en équilibre, reste très délicate à simuler avec un modèle statistique de type $k - \varepsilon$, dont la faiblesse se traduit ici, par des distances supérieures pour le retour à l'équilibre.
### Profils de l'énergie cinétique pour $L_c = 0.2\,\text{m} = 4.3\,\text{h}$

Nous avons reproduit l'étude numérique, pour une clairière de $L_c = 0.2\,\text{m}$. Sur la figure 4.22, nous avons comparé les résultats expérimentaux et le profil de $k$ aux positions $x = -0.1\,\text{m} = -2.1\,\text{h}$, $x = 0\,\text{m}$ et $X = 0\,\text{m}$ (ou $x = 0.2\,\text{m} = 4.3\,\text{h}$). L'objectif est de voir l'influence de la dimension de la clairière (ou l'influence amont-aval des hétérogénéités de surface) sur l'écoulement. Nous avons donc reporté les valeurs de $k$ obtenues pour

<figure><img src="image_60.png" /><figcaption>Figure 4.22 : Énergie cinétique à différentes positions pour une clairière de dimension $L_c = 4.3\,\text{h}$. Comparaison avec les résultats obtenus pour une clairière de $L_c = 21.3\,\text{h}$ aux positions $x = -2.1\,\text{h}$ et $x = 4.3\,\text{h}$ [Raupach 87]</figcaption></figure>

$L_c = 21.3\,\text{h}$ à ces mêmes positions :

- À l'amont de la transition, les profils de l'énergie cinétique sont identiques (figure 4.22). Ce résultat est confirmé par Raupach, qui déduit que la clairière et donc sa dimension n'a pas d'influence amont et ne perturbe pas l'écoulement dans le couvert. Ceci reste valable pour l'ensemble des variables turbulentes;

- En $x = 0\,\text{h}$ et en $x = 4.3\,\text{h}$, l'énergie cinétique est bien calculée, avec des valeurs légèrement surestimées vers le sommet du couvert. Les profils turbulents pour les deux dimensions de la clairière, montrent que l'écoulement n'est pas vraiment modifié par sa dimension, les profils étant très proches. Donc, le couvert ne propage pas de perturbations significatives sur l'écoulement dans la clairière, si ce n'est en empêchant l'écoulement de s'y développer.
<figure><img src="image_61.png" /><figcaption>Figure 4.23 : *Schématisation des couches internes et des sous-couches internes pour une transition rugueux-lisse et lisse-rugueux* [Raupach 87]</figcaption></figure>

La dimension de la clairière n'apparaît pas comme un paramètre déterminant à petite échelle. Cependant, à plus grande échelle, chaque surface induit des perturbations, qui se propagent bien au-delà de l'hétérogénéité et qui provoquent le développement de couches internes (figure 4.23).

## Conclusion

Nous venons de simuler des transitions végétation-clairière-végétation pour deux dimensions de clairière :

- Pour la transition végétation-clairière, les résultats sont très satisfaisants. Ils permettent de montrer le bon comportement de l'utilisation d'un milieu hétérogène fluide-couvert sur une grille unique, pour représenter l'écoulement turbulent en situation hétérogène, à l'échelle des transitions ;

- Pour la transition clairière-végétation, les résultats obtenus confirment la difficulté des modèles statistiques de type $k - \varepsilon$ à simuler une turbulence très intermittente, pour laquelle les champs turbulents instantanés peuvent s'avérer très différents des champs moyens. Cette caractéristique, inhérente aux modèles statistiques, explique les différences relevées pour cette transition, qui est caractérisée par une zone propice aux phénomènes de rafales, dont l'activité turbulente est particulièrement violente.

Néanmoins, les résultats sont appréciables, et le modèle numérique réagit correctement aux transitions végétation-clairière-végétation et aux variations d'échelle que cela implique.
En particulier, la vitesse moyenne horizontale calculée est satisfaisante pour les deux transitions végétation-clairière-végétation.

## 4.3 Dispersion de spores à l'échelle d'un couvert continu

Nous disposons d'un modèle numérique qui permet de simuler la dispersion atmosphérique des particules biotiques à plusieurs échelles, et pour des configurations hétérogènes. Cependant, il existe peu d'études expérimentales sur la dispersion de particules à l'échelle d'un paysage hétérogène. En revanche, à l'échelle plus locale pour des couverts continus, de nombreuses études sont menées qui s'intéressent à la dispersion à courte distance. Afin de valider le modèle de transport des particules à cette échelle, nous nous sommes basés sur les travaux de [Bainbridge 79]. Cette étude a par ailleurs été reprise par [Legg 79] et par [Aylor 82].

### Description de l'expérience

L'expérience de [Bainbridge 79] porte sur l'étude du transport de spores *Lycopodium Clavatum* dans un champ d'orge continu hétérogène. La surface végétale des différentes parties constituant l'orge a été mesurée à l'aide d'un planimètre optique. Toutes ces surfaces ont été assimilées à des éléments géométriques simples (cônes, cylindre, feuille) afin de calculer les densités de surface végétales projetées suivant leur orientation privilégiée (figure 4.24). Une incertitude de l'ordre de 20% est considérée avec cette approche.

<figure><img src="image_62.png" /><figcaption>Figure 4.24 : Densité de surface projetée des éléments végétaux constituant l'orge [Legg 79]</figcaption></figure>
La source émettrice est une source artificielle linéique, composée de six coupelles contenant les spores (figure 4.25). Les coupelles sont séparées d'une distance de 1 mètre et sont placées à hauteur du sommet du couvert. À la perpendiculaire de la source linéique, et en partant de son centre, des mâts de capture des spores ont été disposés à différentes distances et permettent ainsi de mesurer la concentration à 2, 3, 4, 5 et 7 mètres de la source, à des hauteurs comprises entre 0.3 mètre et 4 mètres du sol. Enfin, le dispositif dans sa totalité est placé de telle sorte qu'il soit perpendiculaire à la direction moyenne du vent. Une déviation de la direction du vent par rapport à la ligne de mesure, ou bien une variation trop importante des conditions météorologiques ont été des critères d'interruption de l'expérience, pour ne pas biaiser les mesures de concentration.

<figure><img src="image_63.png" /><figcaption>Figure 4.25 : *Dispositif expérimental vu de dessus* [Bainbridge 79]</figcaption></figure>

La spore *Lycopodium Clavatum* est une particule assimilable à une sphère de diamètre moyen de 32 µm. Sa vitesse de sédimentation a été mesurée par de nombreux auteurs et sa valeur est généralement comprise entre 1.7 cms⁻¹ et 2.1 cms⁻¹.

## Paramètres de simulation

Les conditions météorologiques mesurées à chaque campagne (4 au total) ont été moyennées pour obtenir un champ moyen représentatif [Legg 79]. Le profil du vent moyen utilisé en condition limite en entrée du domaine suit une loi exponentielle dans le couvert, et logarithmique au-dessus (figure 4.26).

La vitesse de frottement, relevée au sommet du couvert ($u_* = 0.33 \text{ ms}^{-1}$), a été utilisée pour initialiser l'énergie cinétique moyenne et le taux de dissipation. Les conditions aux limites en haut du domaine sont des conditions de couche de surface et d'adhérence au sol. La densité de surface frontale du couvert est définie à partir des éléments tige, tête et barbe, qui sont supposés dirigés uniquement suivant le plan vertical et le coefficient de
traînée du couvert, dont la hauteur est de l'ordre de $h = 1 \text{ m}$, est égale à 0.3, ce qui est
 une valeur courante [Guyot 97].

<figure><img src="image_64.png" /><figcaption>Figure 4.26 : Profil de la vitesse moyenne horizontale dans et au-dessus du couvert utilisé en condition limite à l'entrée du domaine [Legg 79]</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>Figure 4.27 : Profil de la concentration, pris à 2 mètres de la source, utilisé en condition limite à l'entrée du domaine [Legg 79]</figcaption></figure>

Nous avons interpolé le profil de concentration mesuré à 2 mètres de la source (points expérimentaux sur la figure 4.27) pour imposer notre source. La condition à l'entrée du domaine est une condition limite de Dirichlet pour la concentration. Les processus de dépôts ont été pris en compte, en considérant que les spores adhèrent par impact sur les
surfaces verticales, et par sédimentation sur les surfaces horizontales avec une vitesse de sédimentation de $1.8 \text{ cms}^{-1}$, valeur couramment utilisée.

Il est à relever que les angles d'inclinaison des éléments végétaux ne sont pas disponibles dans l'article. Or, les dépôts dans le couvert dépendent de leur orientation. Il est toutefois peu probable que les parties végétales telles que la tête, la barbe, et la tige soient concernées par le dépôt par sédimentation. En revanche, les feuilles peuvent intervenir dans les deux processus de dépôts. Pour étudier l'influence de l'inclinaison des feuilles sur le dépôt dans le couvert, nous avons simulé différents cas : des feuilles complètement perpendiculaires au vent moyen aux feuilles inclinées $-\theta = 0, \theta = \pi/4$ et $\theta = \pi/2$, l'angle $\theta$ étant défini suivant le sens trigonométrique.

## Profil de la vitesse moyenne horizontale

<figure><img src="image_66.png" /><figcaption>Figure 4.28 : Profil de la vitesse moyenne dans et au-dessus du couvert calculé numériquement</figcaption></figure>

Sur la figure 4.28, nous avons reporté la vitesse moyenne horizontale calculée par notre simulation une fois l'écoulement développé, ainsi que les valeurs expérimentales et le profil utilisé par [Legg 79]. Malgré l'incertitude sur les densités de surfaces frontales, la vitesse moyenne horizontale obtenue numériquement est proche des conditions météorologiques relevées. Nous pouvons noter que ce profil de vitesse est caractérisé par un point d'inflexion qui n'est pas situé au sommet du couvert. C'est typiquement la signature d'un couvert de densité modérée, dont l'échelle de cisaillement associée à de forts transferts turbulents est de taille comparable à $0.5 \text{ h}$.
### Profils de concentration

Sur la figure 4.29, nous avons tracé les profils de concentration simulés à 3, 4, 5, et 7 mètres de la source, pour différentes inclinaisons des feuilles. Nous y avons reporté les mesures expérimentales de [Bainbridge 79]. Ces résultats sont obtenus en imposant un nombre de Prandtl-Schmidt turbulent égal à $\sigma_c = 0.7$, constant dans le domaine d'étude. Nous avons donc admis que les particules solides se comportent comme des particules fluides, avec une diffusivité identique à celle d'un gaz (pas d'inertie).

De manière globale, nous constatons sur les figures 4.29a, 4.29b, 4.29c et 4.29d que la concentration en spores décroît rapidement en aval de la source, et cela pour les trois angles d'inclinaison. L'activité des éléments végétaux est donc importante et s'avère bien prise en compte par le modèle numérique.

De plus, nous pouvons clairement distinguer deux zones dans le couvert. La partie supérieure du couvert, qui est définie entre $h/2$ et $h$, où la concentration calculée est pratiquement la même pour les 3 angles d'inclinaison des feuilles. Nous pouvons donc conclure que les processus de dépôt par sédimentation, qui ont lieu sur les surfaces horizontales, ne jouent qu'un rôle mineur dans cette partie du couvert. La décroissance de la concentration est donc principalement induite par les processus de dépôt par impact. Par ailleurs, l'efficacité d'impact calculée à l'aide de l'expression de Aylor donne des résultats satisfaisants. Dans cette région, les profils de concentration calculés aux différentes distances de la source sont en bon accord avec les résultats expérimentaux.

En revanche, la concentration dans la partie inférieure du couvert est très différente suivant les 3 angles d'inclinaison des feuilles :

- Pour un angle d'inclinaison de $\theta = 0$, seul le processus de dépôt par sédimentation intervient. Les profils de concentration calculés dans la partie basse du couvert correspondent de manière satisfaisante aux résultats expérimentaux;

- Lorsque les deux processus de dépôt interviennent, c'est à dire pour $\theta = \pi/4$, la concentration calculée près du sol est surestimée. Cependant, elle tend à diminuer de manière significative, au fur et à mesure que la distance $x$ de la source croît. Le dépôt par impact qui intervient dans la partie basse du couvert joue un rôle non négligeable dans les processus de dépôt;

- Lorsque les feuilles sont inclinées dans le plan vertical ($\theta = \pi/2$) la concentration près du sol y est alors clairement surestimée, lorsque le dépôt par sédimentation est
<figure><img src="image_67.png" /><figcaption>Figure 4.29 : Profils de concentration dans et au-dessus du couvert végétal à différentes positions de la source (x exprimé en mètres) [Bainbridge 79, Legg 79]</figcaption></figure>
ignoré. Dans la partie basse du couvert, les processus par impact ne suffisent plus pour décrire la décroissance de la concentration où le dépôt se fait donc principalement par sédimentation.

Les profils de concentration calculés sont très proches des résultats expérimentaux dans les différentes parties du couvert, lorsque la sédimentation est prise en compte sur les feuilles. Le modèle numérique rend compte de la décroissance rapide de la concentration des particules dans le couvert. Notre simulation est en bon accord avec les conclusions de Legg et d'Aylor, qui ont constaté le rôle prépondérant du dépôt par impact dans la partie supérieure du couvert, alors que la sédimentation joue un rôle important dans la partie basse du couvert, où le vent est suffisamment faible pour ne pas affecter ce processus de dépôt.

### Influence de nombre de Prandtl-Schmidt turbulent

<figure><img src="image_68.png" /><figcaption>Figure 4.30 : Profils de concentration dans et au-dessus du couvert végétal à 5 mètres de la source pour différents nombres de Prandtl-Schmidt turbulents. À gauche $\theta = 0$ à droite $\theta = \pi/4$</figcaption></figure>

A l'issue de cette étude, nous avons testé l'influence du nombre de Prandtl-Schmidt sur la diffusion turbulente des spores. Nous avons évoqué que la diffusivité turbulente des scalaires, à hauteur du couvert végétal, est considérée de l'ordre de 2 à 3 fois supérieure à celle de la quantité de mouvement. Nous avons donc répété ces mêmes simulations, en
fixant le nombre de Prandtl-Schmidt à $\sigma_c = 0.5$ ainsi qu’à $\sigma_c = 0.33$ dans la totalité du domaine de calcul.

Nous avons tracé le profil de concentration à 5 mètres de la source sur la figure 4.30. En comparant avec les simulations effectuées précédemment, nous observons une légère amélioration sur l'estimation du profil de la concentration. Cette amélioration est plus marquée pour un nombre de Prandtl-Schmidt $\sigma_c = 0.33$, en particulier entre $z = h$ et $z = 2h$. Nous restons toutefois sceptiques vis-à-vis de l'analyse de ces résultats, compte tenu des incertitudes expérimentales sur les densités de surface végétale et des relations empiriques que nous utilisons pour modéliser les mécanismes de dépôt dans le couvert. Nous conserverons donc la valeur du nombre de Prandtl-Schmidt turbulent à $\sigma_c = 0.7$.

## Conclusion

Nous avons simulé le transport de particules biotiques dans et au-dessus du couvert, à l'aide de l'équation turbulente de convection-diffusion. L'activité des éléments végétaux est décrite en terme de dépôt par impact et par sédimentation. Nous notons que la prise en compte des mécanismes de dépôt, par des expressions simples, permet d'appréhender le champ de concentration dans le couvert. Les résultats obtenus sont très cohérents et satisfaisants. Par ailleurs, l'utilisation d'une modélisation phénoménologique n'est pas en contradiction avec les résultats expérimentaux.

Nous avons également mis en évidence que les dépôts dans le couvert se font suivant des zones privilégiées. Alors que dans la partie inférieure du couvert, le dépôt des particules a lieu principalement par sédimentation, le dépôt par impact peut intervenir dans les différentes parties du couvert, avec un rôle majeur dans la partie supérieure du couvert, où l'écoulement est le plus turbulent.
Conclusion

Dans ce chapitre, nous venons de présenter les validations effectuées sur notre code de calcul. Le code de calcul du laboratoire a été validé par ailleurs, et nous avons exposé ici les parties que nous avons développées.

Les résultats numériques relatifs à l'étude du vent sur des surfaces hétérogènes simples, ainsi que la modélisation de la couche à flux constant sont très satisfaisants, en 2D comme en 3D. Cette partie constitue les premiers éléments d'un modèle à l'échelle sub-méso.

À l'échelle du couvert, les simulations ont été plus délicates à mener. Dans un premier temps, nous nous sommes attachés à vérifier que le formalisme mathématique défini pour le modèle de canopée permet de mettre en valeur la spécificité des écoulements à cette échelle. Cette étape indispensable a été réalisée en s'appuyant sur les mesures de [Brunet 94], disponibles loin de la discontinuité. Nos résultats numériques ont révélé que le modèle de turbulence est particulièrement sensible au choix de la constante empirique de la dissipation de sillage du taux de dissipation. Cette constante tend à inhiber ou à accentuer les nombreux termes qui gouvernent le mouvement turbulent dans la sous-couche rugueuse. Une étude de sensibilité a été menée pour identifier les sources d'erreurs, et notre modèle ainsi défini est pour l'ensemble des variables satisfaisant.

Puis, nous avons simulé un écoulement pour un cas complexe, faisant intervenir des transitions végétation-clairière-végétation. Les résultats relatifs à la transition végétation-clairière permettent d'affirmer que l'utilisation d'un système d'équations modifiées selon l'échelle d'observation est efficace. La transition clairière-végétation est plus délicate à simuler, car elle évoque une turbulence très intermittente. Les résultats sont toutefois intéressants, et restent conformes aux choix du modèle de turbulence, connu pour ces difficultés à reproduire cette turbulence. Il est à relever que ce modèle de turbulence, de par ces défauts, met en avant les phénomènes physiques qui interviennent dans ces configurations et souligne ainsi les efforts qui doivent être effectués pour améliorer la qualité des résultats attendus.

Enfin, nous avons simulé la dispersion de spores dans et au-dessus d'un couvert hétérogène infini. Les relations empiriques utilisées pour prendre en compte les mécanismes de dépôt dans le couvert se sont révélées efficaces, puisque les champs de concentration calculés sont très proches des résultats expérimentaux. Ce dernier cas de validation confirme le bon comportement général du code de calcul.
<figure><img src="image_69.png" /></figure>
# Chapitre 5

## Études numériques sur la dispersion de particules biotiques

### 5.1 Introduction

La dispersion atmosphérique des particules biotiques dépend de nombreux paramètres fortement corrélés puisqu'elle est la conséquence des mécanismes couplés de la convection et de la diffusion turbulente. Or, le vent, l'intensité turbulente et par conséquent les transferts turbulents dans la couche de surface, sont en grande partie tributaires de la variabilité spatiale de la surface terrestre, à petite comme à grande échelle. Nous pouvons d'emblée concevoir que la structure d'un paysage ait des effets sur la dispersion et sur l'échelle de dispersion.

Ce point que nous abordons là est une des questions essentielles, débattue par un très large public (médias, scientifiques, agriculteurs...), avec l'introduction volontaire dans l'environnement d'organismes génétiquement modifiés. À dire vrai, les préoccupations relatives à la contamination génétique d'un champ à l'autre existent déjà pour les champs de production de semences. Pour garantir la pureté des semences, les principaux moyens de prévention utilisés à l'heure actuelle s'appuient sur deux considérations. Une isolation par la distance des champs de cultures (tableau 5.1) nommée isolation synchronique, stratégie qui considère que le flux de gène est limité par des barrières naturelles, telles que la distance, la direction du vent, le relief... La seconde stratégie utilisée est une isolation diachronique qui se traduit par le choix d'espèces ayant des périodes de floraison différentes dans un même espace. Ces mesures de prévention visant à limiter les flux de gène prennent une dimension accrue, avec la mise en place d'une politique de cultures à grande échelle
d'OGM, et cela, aussi bien pour l'environnement, que pour des filières de production
garanties sans OGM.

<table><caption>Tableau 5.1 : *Distances d'isolement en mètres utilisées en France et en Grande-Bretagne pour garantir la pureté des semences produites. Pour la France, il se distingue les conditions pour la production de semences de base et pour la production de semences certifiées (entre parenthèses). Données issues de GNIS et de Raybould et Cray (1993). Tableau repris de Klein [Klein 00]*</caption>
   <thead>
    <tr>
     <td>
      Espèce
     </td>
     <td>
      France
     </td>
     <td>
      Grande-Bretagne
     </td>
    </tr>
   </thead>
   <tbody>
    <tr>
     <td>
      Betterave (Beta vulgaris)
     </td>
     <td>
      1000 (600)
     </td>
     <td>
      3200
     </td>
    </tr>
    <tr>
     <td>
      Colza (Brassica rapa)
     </td>
     <td>
      1000 (400)
     </td>
     <td>
      900
     </td>
    </tr>
    <tr>
     <td>
      Oignon (Allium cepa)
     </td>
     <td>
      500 (300)
     </td>
     <td>
      900
     </td>
    </tr>
    <tr>
     <td>
      Maïs (Zea mays)
     </td>
     <td>
      400 (200)
     </td>
     <td>
      180
     </td>
    </tr>
    <tr>
     <td>
      Radis (Raphanus sativus)
     </td>
     <td>
      500 (300)
     </td>
     <td>
      270 – 300
     </td>
    </tr>
   </tbody>
  </table>

Ainsi, s'il est évident que la structure d'un paysage influence la dissémination des parti-
cules, ces manifestations ne sont pas faciles à quantifier pour autant. Pour s'en convaincre,
il suffit de comparer les distances d'isolement pratiquées par la Grande-Bretagne et la
France, qui sont pour certaines espèces considérablement différentes (tableau 5.1). La mo-
délisation de la dispersion des particules peut donc constituer une voie d'investigation
pour évaluer et quantifier les potentialités de contamination dans un paysage agricole, et
participer à la mise en place de stratégies d'aménagement durable du territoire.

L'objet de ce chapitre, qui s'articule autour de trois sections, est donc d'apporter des éléments de compréhension supplémentaires sur la dispersion atmosphérique, puis sa relation avec la structure d'un paysage agricole. La première section est consacrée à définir les paramètres qui influencent la dispersion des particules dans la couche de surface. La seconde section porte sur la dispersion de particules de pollen de maïs pour une configuration hétérogène, plus concrète et complexe, comme elle se rencontre dans un paysage agricole. En effet, nous mènerons cette simulation en nous plaçant dans des conditions similaires au dispositif expérimental de l'AGPM<fnref n="1" />, agencé spécifiquement pour étudier

<footnote n="1">1. Association Générale des Producteurs de Maïs</footnote>
l'effet des discontinuités dans la végétation. Enfin, en dernière section, nous extrapolerons
ce dispositif de référence à différentes configurations, faisant intervenir ou non des discon-
tinuités dans la végétation. Ces simulations visent à analyser l'influence de la nature de
la surface sur la dispersion.

## 5.2 La dispersion dans la couche de surface

La dispersion atmosphérique des particules biotiques, et de toutes grandeurs transpor-
tables est influencée par de nombreux paramètres. Afin d'identifier ces paramètres, nous
avons mené une étude de sensibilité du code de calcul en 2 dimensions. Cette étude nu-
mérique porte sur la dispersion à courte distance dans la couche de surface. Précisons que
nous ne disposons pas de référence expérimentale, mais que nous nous sommes inspirés
des travaux numériques réalisés par [McCartney 90].

### 5.2.1 Paramètres de simulation

Nous avons simulé un écoulement turbulent dans la couche de surface homogène, en
stratification neutre. Les forces de Coriolis sont donc négligées. Les caractéristiques de
la surface sont intégrées à l'aide des paramètres aérodynamiques $z_0$ et $d$. Les conditions
limites utilisées s'appuient sur ces considérations et sont, par conséquent, des conditions
classiques de couche de surface.

Pour le transport turbulent des particules pesantes, nous avons imposé à l'entrée du
domaine de calcul, une source émettrice placée à la hauteur $H$. La concentration est dé-
crite par une courbe gaussienne. La condition limite au sol pour le champ de concentration
utilise une vitesse de dépôt prise égale à deux fois la vitesse de sédimentation des par-
ticules, puisque le sol est rugueux [McCartney 91]. Nous avons étudié successivement les
effets découplés de chacun des facteurs suivant sur la dispersion des particules :

- La vitesse de sédimentation des particules $v_s$. À cet effet, différentes classes de particules biotiques sont considérées;

- La hauteur de la source $H$;

- La vitesse de frottement $u_*$;

- La rugosité de la surface $z_0$.
Toutes les simulations numériques réalisées utilisent un maillage identique. Le domaine de calcul s'étend sur 50 mètres, suivant la direction horizontale et 100 mètres, suivant la hauteur. Les paramètres physiques spécifiques à chaque simulation sont mentionnés en plus de la présentation des résultats.

### 5.2.2 Les résultats numériques : les flux et les dépôts

Les résultats de nos simulations sont présentés sous forme de courbes de flux de particules convectées et de flux de particules déposées. Ces flux sont définis par :

- Le flux intégré de particules convectées par l'écoulement moyen, passant à travers chaque plan vertical positionné en $x$, est calculé par l' expression suivante :

$$ \Phi_x(x) = \int_{z_p}^{z_H} \bar{u}(x, z) \bar{C}(x, z) dz $$

Il est d'usage de rapporter $\Phi_x(x)$ [$particule m^{-1} s^{-1}$] au débit de la source, car l'évolution du rapport caractérise l'épuisement du panache initial. Le débit (pour une largeur unitaire) de la source est estimé par la valeur du flux intégré pris à l'entrée du domaine de calcul, ce qui correspond par ailleurs à sa valeur maximale :

$$ \Phi_{\text{émis}} = \int_{z_p}^{z_H} \bar{u}(0, z) \bar{C}(0, z) dz = \text{Max} [\Phi_x(x), 0] \quad (5.1) $$

$$ \Phi_{\text{aéro}}(x) = \frac{\Phi_x(x)}{\Phi_{\text{émis}}} \quad (5.2) $$

$\Phi_{\text{aéro}}(x)$ (en % par rapport à l'émission) permet de renseigner sur la quantité de particules qui sont transportées par convection. Cette grandeur nous donne donc une indication sur la proportion de particules présentes au sein de l'écoulement, susceptibles d'être dispersées au-delà du domaine de calcul et par la même, définit un potentiel de dispersion. Bien que $\Phi_{\text{aéro}}(x)$ donne une idée sur le potentiel de dispersion des particules, elle n'apporte pas d'indication sur la distance de dispersion elle-même. Enfin, signalons également que $\Phi_{\text{aéro}}(x)^{-1}$ qui correspond donc au rapport de la quantité de particules débitées par la source à la quantité de particules passant à travers une section droite, définit l'appauvrissement moyen d'un panache ;

- Pour les dépôts, nous rappelons que le flux de particules déposées [$particule m^{-2} s^{-1}$] est donné par :

$$ F_d(x) = v_d \bar{C}(x, z_p) $$
$v_d$ étant la vitesse de dépôt et $\overline{C}(x, z_p)$ la concentration à proximité du sol. Nous avons adimensionné $F_d$ par l'émission de la source :

$$ \Phi_{d\acute{e}p\acute{o}t}(x) = \frac{F_d(x)}{\Phi_{\acute{e}mis}} \qquad (5.3) $$

$\Phi_{d\acute{e}p\acute{o}t}(x)$ (en % $m^{-1}$) est le flux de matière déposée par mètre par rapport à l'émission de la source. Il rend compte des dépôts au sol qui s'effectuent à chaque pas d'espace $\Delta x$.

D'autre part, nous avons vu que l'équation de transport des particules est basée sur le principe de conservation du champ de concentration. Les particules étant pesantes, le bilan de conservation des flux est vérifié si le dépôt cumulé et le flux de particules convectées par le fluide vérifient à chaque position $x$ du domaine de calcul, la relation :

$$ \int_{0}^{x} v_d \overline{C}(x, z_p) dx = \int_{z_p}^{z_H} \bar{u}(0, z) \overline{C}(0, z) dz - \int_{z_p}^{z_H} \bar{u}(x, z) \overline{C}(x, z) dz \quad (5.4) $$

Ce qui se traduit, en reprenant les grandeurs et les notations introduites plus haut, par :

$$ \int_{0}^{x} \Phi_{d\acute{e}p\acute{o}t}(x) dx = 1 - \Phi_{a\acute{e}ro}(x) \qquad (5.5) $$

### Influence de la vitesse de chute

La position de la source émettrice est placée en $H = 3 m$. Les vitesses de chutes choisies pour cette simulation sont respectivement $v_s = 5 cm s^{-1}$, $v_s = 10 cm s^{-1}$, $v_s = 20 cm s^{-1}$ et enfin $v_s = 40 cm s^{-1}$. Ces valeurs correspondent à des spores, des pollens, et des graines. Les conditions des simulations numériques (maillage, écoulement...) sont identiques pour toutes les classes de particules. Le sol est rugueux de rugosité $z_0 = 0.05 m$ et la vitesse de frottement est $u_* = 0.3 m s^{-1}$.

Sur la figure 5.1, nous avons représenté la fraction de particules présentes dans l'air ainsi que le dépôt par mètre en pourcentage par rapport à l'émission. De manière claire, il apparaît que le pourcentage de particules présentes dans l'écoulement décroît plus rapidement lorsque la vitesse de chute des particules est élevée. En considérant les cas extrêmes, 80% des particules dont la vitesse de chute est de $v_s = 5 cm s^{-1}$ sont encore présentes dans l'écoulement, alors qu'il ne reste plus que 20% de particules dont la vitesse de chute est de $v_s = 40 cm s^{-1}$. En même temps, le dépôt des particules lourdes a lieu très rapidement à l'aval de la source alors qu'il se fait de manière plus étalée pour les particules fines.
La vitesse de sédimentation est une caractéristique inhérente à la particule et demeure un paramètre important dans les modèles de transport des particules pesantes. Son utilisation est justifiée à plus d'un titre, car elle permet simplement de tenir compte de la géométrie parfois complexe des particules biotiques, et de leurs capacités à être transportées efficacement par le vent. Sa valeur fixe donc le potentiel de dispersion d'une particule. Les particules fines ont un potentiel de dispersion particulièrement élevé, puisque leur mouvement est fortement dominé par le mouvement turbulent de l'air. Par conséquent, elles sont susceptibles d'être transportées sur de longues distances.

<figure><img src="image_70.png" /><figcaption>Figure 5.1 : Flux convectifs intégrés et dépôts par mètre (en %). Influence de la vitesse de sédimentation sur la distance de dispersion. Les paramètres numériques sont $u_* = 0.3 \text{ m s}^{-1}$, $z_0 = 0.05 \text{ m}$ et la source est positionnée en $H = 3 \text{ m}$</figcaption></figure>

## Influence de la hauteur de la source

La hauteur de la source émettrice joue également un rôle important, car elle détermine la dispersion initiale d'un panache. Il a été constaté que pour des sources positionnées à haute altitude, la distance horizontale à laquelle la concentration près du sol atteint son maximum, s'accroît. Les dépôts ont donc lieu loin à l'aval de la source, puisque la concentration près du sol fixe l'intensité des dépôts.
<figure><img src="image_71.png" /><figcaption>Figure 5.2 : Flux convectifs intégrés et dépôts par mètre (en %). Influence de la hauteur de la source. Les paramètres numériques sont $u_* = 0.3 \text{ m s}^{-1}$ et $z_0 = 0.05 \text{ m}$, $v_s = 20 \text{ cm s}^{-1}$</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>Figure 5.3 : Flux convectifs intégrés et dépôts par mètre (en %). Influence de la vitesse de frottement. Les paramètres numériques sont $v_s = 20 \text{ cm s}^{-1}$, $z_0 = 0.05 \text{ m}$ et la source est positionnée en $H = 3 \text{ m}$</figcaption></figure>