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

HAL
open science

# Développement d'un modèle numérique de dispersion atmosphérique de particules à l'échelle d'un paysage hétérogène.

Hadjira Foudhil

► To cite this version:

Hadjira Foudhil. Développement d'un modèle numérique de dispersion atmosphérique de particules à l'échelle d'un paysage hétérogène. Sciences de l'ingénieur [physics]. Université Bordeaux 1, 2002. Français. NNT: 2002BOR12525. tel-04486375

HAL Id: tel-04486375

https://theses.hal.science/tel-04486375

Submitted on 1 Mar 2024

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

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

Public Domain
THÈSE

PRÉSENTÉE À

L'UNIVERSITÉ BORDEAUX I

ÉCOLE DOCTORALE
DES SCIENCES PHYSIQUES ET DE L'INGÉNIEUR

par Hadjira FOUDHIL

POUR OBTENIR LE GRADE DE

DOCTEUR

SPÉCIALITÉ : MÉCANIQUE

DÉVELOPPEMENT D'UN MODÈLE NUMÉRIQUE DE
DISPERSION ATMOSPHÉRIQUE DE PARTICULES À
L'ÉCHELLE D'UN PAYSAGE HÉTÉROGÈNE

Soutenue le : 27 juin 2002

Après avis de :

MM. H-C. BOISSON, Directeur de Recherche CNRS, IMFT Toulouse Rapporteur

P. MESTAYER, Directeur de Recherche CNRS, ECN Nantes Rapporteur

Devant la Commission d'examen formée de :

MM. V. MILISIC, Professeur, ISTAB Bordeaux Président

Y. BRUNET, Directeur de Recherche INRA, Bordeaux Directeur

J-P. CALTAGIRONE, Professeur, Université Bordeaux I Directeur

A. DRUILHET, Directeur de Recherche CNRS, OMP Toulouse Examinateur

C. LAVIGNE, Maître de Conférence, Université Paris XI Examinatrice

A. TUZET, Chargée de Recherche INRA, Grignon Examinatrice
<figure><img src="image_2.png" /></figure>
*Pour progresser, il ne suffit pas de vouloir agir, il faut d'abord savoir dans quel sens agir.*

Gustave Le Bon
<figure><img src="image_3.png" /></figure>
Remerciements

Le travail présenté dans ce mémoire est l’aboutissement d’une étroite collaboration entre l’Unité de Bioclimatologie de l’INRA Bordeaux et le laboratoire de Modélisation Avancée des Systèmes Thermiques et des Écoulements Réels (MASTER) de l’École Nationale Supérieure de Chimie Physique de Bordeaux (ENSCPB). Il a été effectué dans le cadre d’une bourse INRA-Région Aquitaine.

En premier lieu, je tiens à manifester ma sincère reconnaissance envers mes deux directeurs de thèse, sans qui ce travail n’aurait jamais vu le jour. Je remercie Yves Brunet, directeur de recherche et directeur de l’Unité de Bioclimatologie de l’INRA Bordeaux, pour l’aide précieuse et la confiance qu’il m’a accordée durant ces années de thèse, et qu’il m’accorde encore aujourd’hui. Je remercie également Jean-Paul Caltagirone, professeur à l’Université Bordeaux I et directeur alors du laboratoire MASTER, pour nous avoir suivis dans cette aventure en m’accueillant et en mettant à ma disposition tous les moyens et les compétences du laboratoire (sans parler d’AQUILON et de ses fameuses sources...).

Je tiens à exprimer toute ma gratitude envers Henri-Claude Boisson, directeur de recherche à l’IMFT, pour avoir accepté d’être rapporteur de cette thèse et pour l’intérêt qu’il a porté à ce travail.

Je remercie vivement Patrice Mestayer, directeur de recherche à l’École Centrale Nantes, pour son analyse fructueuse de ce travail. Ses nombreuses remarques m’auront permis d’améliorer la qualité de ce mémoire. Je lui sais gré d’avoir accepté d’être rapporteur de cette thèse.

J’adresse mes remerciements au professeur Vladan Milisic pour m’avoir fait l’honneur et le plaisir de présider mon jury ; à Aimé Druilhet, directeur de recherche à l’OMP Toulouse, pour la pertinence de ses propos, à Claire Lavigne, maître de conférence à l’Université Paris XI et à Andrée Tuzet, chargée de recherche à l’INRA Grignon, pour avoir bien voulu juger ce travail et participer à cette commission.

Durant ces trois années (et demie), j’ai eu la chance d’être rattachée à deux laboratoires
où j'ai pu côtoyer de nombreuses personnes qui ont toutes participé à ce travail, que ce soit directement par leurs compétences scientifiques ou indirectement par leurs qualités humaines.

Je remercie chaleureusement tous les membres de l'Unité de Bioclimatologie pour l'accueil et pour l'intérêt qu'ils m'ont porté. En particulier, je remercie Éric Lamaud pour avoir participé à mon comité de thèse ; Charles Valancogne et Sylvia Dayau avec qui j'ai eu le plaisir de travailler lors de la campagne « Arbanat» ; Jean-Pierre Lagarde, Dominique Guyon, Paul Berbigier, Dewar Roderick, Jean-Pierre Wigneron et Mark Irvine pour avoir prêté une oreille attentive à mes répétitions et contribué ainsi à parfaire ma présentation de thèse (tant redoutée) ; Gaston Courrier et Jean-Marc Bonnefond pour leur gentillesse, Anne-Marie Bouchon, Josyane Bonnet et Françoise Prud'Homme pour leur serviabilité ; Didier Garrigou et Alain Kruszewski pour les moments de détente au volley.

Je remercie tout aussi chaleureusement les permanents du laboratoire MASTER pour leur aide, Éric Arquis notre nouveau directeur, Stéphane Vincent notre tout récent maître de conférence à l'humour décapant et Jean Toutain pour ses multiples instigations à mes préoccupations d'un moment. Pour l'informatique, je remercie Sandrine et Bernard qui m'a fait découvrir un monde merveilleux, sans bugs ni virus, encore nommé UNIX.

Je n'oublierai pas de remercier également Agnès Calonnec de l'Unité de Santé Végétale à l'INRA Bordeaux, l'instigatrice de la campagne « Arbanat» ; Brigitte Durand et Nathalie Jarosz, dite la végétarienne (on y retourne quand tu veux), toutes deux à l'Unité Environnement et Grandes Cultures à l'INRA Grignon, à qui je me suis jointe l'espace d'un mois pour œuvrer sur une campagne de mesure très... maïs. En plus de l'expérience fort plaisante, j'y ai appris certains choses inattendues, telles que déceler l'origine d'une défaillance d'un digitaliseur magnétique. Surprenant.

Place maintenant à la communauté des thésards du laboratoire MASTER que j'ai fréquentée au quotidien et que j'affectionne sincèrement.

Mes pensées vont tout d'abord vers les anciens qui sont, pour la plupart, partis vers d'autres rivages : l'éternel amoureux Manu, Jeff qui m'a fait découvrir l'Algarve, Rémy à qui j'adresse des salutations très sportives et la douce Naïma qui n'a tout de même pas hésité à m'abandonner parmi tous ces mâles pour s'installer au pays du soleil. Il y a également les expatriés, Hugues en qui j'ai trouvé un véritable ami et à qui je dis hippopotame, hippopotame, et mon Gégé avec qui j'ai partagé bien plus que le bureau S-18. Il reste heureusement le généreux Stéphane, qui m'a couvé et m'a permis rapidement de prendre en main AQUILON. Ses compétences dans ce domaine et son rôle de médiateur
au laboratoire font qu'il mérite amplement son futur poste d'ingénieur de recherche. Je le souhaite pour lui et pour son adorable petite famille.

Pour ce qui est des nouveaux thésards, j'adresse un petit clin d'œil à Guiguit, le multitâche, à défaut de comprendre tous ses jeux de mots, j'ai compris qu'il est loin du canard boiteux qu'il prétend. Il est également important de souligner ses performances oratoires journalières, car je crois bien qu'il parle plus que moi (si, si, ça existe). Cédric, l'hyper sensible, d'une gentillesse à toute épreuve, sache que l'essentiel est d'avoir les pieds qui touchent le sol (et NA pour Guiguit). Taz, un mélange subtile « du bon, de la brute et du truand », qui s'adonne à toutes sortes d'expériences bizarres pourvu qu'elles soient un tantinet douloureuses (maintenant en thèse, tu seras servi). Frédéric, encore un Breton à l'esprit torturé et tortueux mais qui finit toujours par trouver son chemin (et quel chemin !). Boris, mon nouveau co-locataire S-18, surfeur, jogger, planchiste et j'en passe, à qui j'ai dû en faire voir de toutes les couleurs pendant ma rédaction (pardon). Damien, passionné et engagé pour ses principes, à qui il arrive toutes sortes de mésaventures lors de ses périoles au bout du monde. Tout aussi fougueux, Claude l'amateur de bonne cuisine exotique et de jeux vidéo. Grégoire, sympathique mais réservé, que j'espère bien découvrir cette année. Coté miss, l'attachante Cendrille et l'adorable Delphine, les seules thésardes, une espèce en survie au laboratoire, à qui je souhaite beaucoup de réussite.

Bien sûr, j'ai une pensée pour ma famille...

Je remercie toute ma famille pour son soutien moral, ma tante Zohra, mes cousines que j'adore et que je ne citerai pas (Férouze, Assia, Kenza la future maman à qui je souhaite bonne chance, Hassiba et LES El-Gholam Nadia et Rafika... Ouff).

Je remercie ma seconde famille ici, Patricia, Sarah et Jihane qui ont bien voulu m'adop-
ter et qui ont fait de Bordeaux, une ville accueillante et chaleureuse.

Je remercie mes parents et ma sœur, que j'aime profondément. Papa, merci pour ce que tu es, pour ce que tu as fait pour moi, et pour les milliers de fois où tu as su trouver les mots justes face à mes coups de blues (et Dieu sait s'il y en a eu). Maman, je ne peux trouver de mots plus doux. Cette réussite, aujourd'hui, est la tienne, merci. Thania, ma petite sœur chérie, je te remercie énormément pour ton soutien. Je mesure et te félicite pour ton courage, car je sais combien notre séparation a été douloureuse à vivre au quotidien, pour toi comme pour moi.

Frédéric, mon acolyte, je te dis simplement merci d'être là...

Ma dernière pensée ira à mon cher oncle Merzak et à mon ami Krime, qui nous ont quittés trop tôt et à qui je dédie ce mémoire.
<figure><img src="image_4.png" /></figure>
Résumé :

Ce travail de thèse a pour objet le développement d'un modèle numérique de dispersion de particules biotiques, basé sur la résolution des équations de la couche limite atmosphérique. L'accent est porté sur l'élaboration d'un code capable de simuler des écoulements pour des configurations de paysages hétérogènes, car la dispersion des particules est fortement influencée par la structure même du paysage. Le modèle dynamique utilise les équations de Reynolds simplifiées dans le cadre des hypothèses de Boussinesq et la turbulence est traitée par une approche statistique (modèles $k-l$ et $k-\varepsilon$). Le système d'équations mis en place permet de modéliser les principaux processus physiques à deux échelles différentes : la modélisation à grande échelle s'appuie sur la théorie de similitude de Monin et Obukhov, qui permet d'élaborer des lois de paroi pour simuler les transferts turbulents entre la surface apparente et la basse atmosphère ; la modélisation à l'échelle du couvert amène à résoudre explicitement l'écoulement dans la sous-couche rugueuse, définie par une canopée végétale, jusqu'à la surface du sol. La dispersion des particules est traitée par une approche eulérienne, qui considère l'ensemble des particules dans l'écoulement comme un champ de concentration continu, dont l'évolution est décrite par une équation de convection-diffusion. Le code de calcul est d'abord validé à partir de données expérimentales acquises dans diverses configurations, en soufflerie comme in situ : transition de rugosité de surface, turbulence en couvert homogène, alternances couvert-clairière-couvert, diffusion de spores en couvert homogène, diffusion de pollen en paysage discontinu. Une série de simulations est ensuite menée pour analyser l'influence de la structure de divers paysages complexes sur la dispersion de particules.

Mots-clés : dispersion atmosphérique, turbulence, écoulement atmosphérique, simulation numérique, particule biotique, pollen, couvert végétal, paysage, discontinuité

Summary :

The aim of this thesis is to develop a numerical model for the dispersal of biological particles, based on the resolution of the atmospheric boundary layer equations. As the dispersal of biological particles is strongly influenced by the landscape structure, one important point is to elaborate a tool able to simulate flows over heterogeneous landscape configurations. The aerodynamic model uses simplified Reynolds equations with the Boussinesq assumption and turbulence is treated by a statistical approach ($k-l$ and $k-\varepsilon$ models). The set of equations used here allows the main physical processes to be represented at two different scales. Large-scale modelling is based on the Monin-Obukhov similarity theory, which provides wall boundary conditions to model turbulent transfer between the apparent surface and the lower atmosphere. Canopy-scale modelling amounts to solve the flow equations in the entire roughness sublayer down to the soil surface. The model for particle dispersal is based on an Eulerian approach, which considers the set of particles in the flow as a continuous concentration field whose evolution is described by a convection-diffusion equation. In a first step the numerical model is validated against data sets from several in situ and wind-tunnel experiments: smooth-to-rough surface change, turbulence in a homogeneous canopy, forest-clearing sequences, spore diffusion in a homogeneous canopy, pollen diffusion across a discontinuous landscape. A series of simulations are then carried out to analyse the influence of various landscape structures on particle dispersal.

Keywords : atmospheric dispersal, turbulence, atmospheric flow, numerical simulation, biological particles, pollen, plant canopy, landscape, discontinuity
<figure><img src="image_5.png" /></figure>
# Sommaire

Nomenclature 1

Introduction Générale 9

1 Modélisation physique des phénomènes 13

1.1 La couche limite atmosphérique 15

1.1.1 Structure de la couche limite atmosphérique 15

1.1.2 Problème d'échelles. Définition des échelles de l'étude 16

1.1.3 Les équations instantanées simplifiées 19

1.2 Prise en compte du mouvement turbulent 21

1.2.1 Quelques approches pour la modélisation de la turbulence 22

1.2.2 Description statistique de la turbulence 23

1.2.3 Traitement statistique de Reynolds 25

1.3 Les modèles de fermeture 26

1.3.1 Modélisation phénoménologique des corrélations 27

1.3.2 Système d'équations modélisées 29

1.3.3 Modèle de fermeture à une équation de transport 30

1.3.4 Modèle de fermeture à deux équations de transport 32

1.4 Modélisation des échanges turbulents près du sol 35

1.4.1 La couche de surface sur un sol homogène: hypothèses et équations 36

1.4.2 Interaction dynamique de l'écoulement avec une paroi 37

1.4.3 Interaction thermique de l'écoulement avec une paroi 41

1.4.4 Paramètres de rugosités des surfaces naturelles 44

1.5 Les effets de la stratification thermique de l'atmosphère 47

1.5.1 Mesure de la stabilité et nature des transferts convectifs 47

1.5.2 La théorie de similitude de Monin et Obukhov 50
<table>
  <tbody>
    <tr>
      <td>1.5.3</td>
      <td>Relations de similitude dans la couche de surface</td>
      <td>51</td>
    </tr>
    <tr>
      <td>1.5.4</td>
      <td>Les fonctions universelles de Dyer</td>
      <td>52</td>
    </tr>
    <tr>
      <td>1.6</td>
      <td>Modélisation de l'écoulement à l'échelle d'un couvert végétal</td>
      <td>54</td>
    </tr>
    <tr>
      <td>1.6.1</td>
      <td>La turbulence et le vent</td>
      <td>55</td>
    </tr>
    <tr>
      <td>1.6.2</td>
      <td>Moyenne spatio-temporelle</td>
      <td>59</td>
    </tr>
    <tr>
      <td>1.6.3</td>
      <td>Équations du mouvement dans le couvert végétal</td>
      <td>61</td>
    </tr>
    <tr>
      <td>1.6.4</td>
      <td>Modélisation du mouvement turbulent dans le couvert</td>
      <td>64</td>
    </tr>
    <tr>
      <td><b>2</b></td>
      <td><b>Dispersion atmosphérique des particules biotiques</b></td>
      <td><b>71</b></td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Processus de dispersion dans l'atmosphère</td>
      <td>73</td>
    </tr>
    <tr>
      <td>2.1.1</td>
      <td>La libération des particules biotiques</td>
      <td>74</td>
    </tr>
    <tr>
      <td>2.1.2</td>
      <td>Le transport turbulent</td>
      <td>77</td>
    </tr>
    <tr>
      <td>2.1.3</td>
      <td>Sédimentation et impact, les mécanismes de dépôt</td>
      <td>79</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>L'aérodynamique des particules biotiques</td>
      <td>80</td>
    </tr>
    <tr>
      <td>2.2.1</td>
      <td>Vitesse de sédimentation des particules biotiques</td>
      <td>80</td>
    </tr>
    <tr>
      <td>2.2.2</td>
      <td>Nombre de Stokes</td>
      <td>82</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Revue sur les principaux modèles de dispersion</td>
      <td>83</td>
    </tr>
    <tr>
      <td>2.3.1</td>
      <td>Modèles Gaussiens</td>
      <td>84</td>
    </tr>
    <tr>
      <td>2.3.2</td>
      <td>Modèles Lagrangiens</td>
      <td>86</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Approche Eulérienne</td>
      <td>91</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Équation instantanée de convection-diffusion</td>
      <td>92</td>
    </tr>
    <tr>
      <td>2.5.1</td>
      <td>Dispersion dans la couche limite atmosphérique</td>
      <td>93</td>
    </tr>
    <tr>
      <td>2.5.2</td>
      <td>Dispersion dans le couvert végétal</td>
      <td>95</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Les processus de dépôt</td>
      <td>101</td>
    </tr>
    <tr>
      <td>2.6.1</td>
      <td>Les dépôts dans la couche limite atmosphérique</td>
      <td>101</td>
    </tr>
    <tr>
      <td>2.6.2</td>
      <td>Les dépôts dans le couvert végétal</td>
      <td>104</td>
    </tr>
    <tr>
      <td><b>3</b></td>
      <td><b>Méthodologie numérique du code de calcul</b></td>
      <td><b>111</b></td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>D'un problème continu à un problème discret</td>
      <td>113</td>
    </tr>
    <tr>
      <td>3.1.1</td>
      <td>Une grille unique</td>
      <td>113</td>
    </tr>
    <tr>
      <td>3.1.2</td>
      <td>Choix des modèles de turbulence et récapitulatifs du système d'équations final</td>
      <td>114</td>
    </tr>
    <tr>
      <td>3.1.3</td>
      <td>Traitement des conditions aux limites</td>
      <td>117</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Procédure numérique des équations</td>
      <td>118</td>
    </tr>
    <tr>
      <td>3.2.1</td>
      <td>Discrétisation temporelle</td>
      <td>118</td>
    </tr>
  </tbody>
</table>
<table><tr><td>3.2.2</td><td>Découplage et linéarisation des équations</td><td>119</td></tr><tr><td>3.2.3</td><td>Couplage vitesse-pression : méthode du Lagrangien Augmenté</td><td>121</td></tr><tr><td>3.2.4</td><td>Discrétisation spatiale : méthode des volumes finis</td><td>122</td></tr><tr><td>3.2.5</td><td>Résolution numérique des systèmes linéaires</td><td>125</td></tr><tr><td>3.3</td><td>Conditions limites au sol</td><td>126</td></tr><tr><td>3.3.1</td><td>Loi de paroi dynamique</td><td>126</td></tr><tr><td>3.3.2</td><td>Loi de paroi thermique</td><td>129</td></tr><tr><td>3.3.3</td><td>Traitement des variables turbulentes</td><td>130</td></tr><tr><td>3.4</td><td>La couche de surface et les conditions limites qui en découlent</td><td>131</td></tr><tr><td>3.5</td><td>Généralisations des conditions limites en stratification thermique</td><td>132</td></tr><tr><td>3.6</td><td>Les conditions limites pour l'équation de transport des particules</td><td>134</td></tr></table>

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

<table><tr><td>4.1</td><td>Les transitions dans la couche de surface</td><td>143</td></tr><tr><td>4.1.1</td><td>Transition lisse-rugueux</td><td>143</td></tr><tr><td>4.1.2</td><td>Saut de température pariétale</td><td>152</td></tr><tr><td>4.2</td><td>Simulations à l'échelle du couvert végétal</td><td>157</td></tr><tr><td>4.2.1</td><td>Écoulement dans un couvert végétal homogène</td><td>157</td></tr><tr><td>4.2.2</td><td>Écoulement dans des transitions végétation-clairière-végétation</td><td>166</td></tr><tr><td>4.3</td><td>Dispersion de spores à l'échelle d'un couvert continu</td><td>177</td></tr></table>

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

<table><tr><td>5.1</td><td>Introduction</td><td>187</td></tr><tr><td>5.2</td><td>La dispersion dans la couche de surface</td><td>189</td></tr><tr><td>5.2.1</td><td>Paramètres de simulation</td><td>189</td></tr><tr><td>5.2.2</td><td>Les résultats numériques : les flux et les dépôts</td><td>190</td></tr><tr><td>5.2.3</td><td>Discussion</td><td>198</td></tr><tr><td>5.3</td><td>Dispersion de pollen de maïs pour un dispositif discontinu</td><td>200</td></tr><tr><td>5.3.1</td><td>Présentation du dispositif expérimental discontinu</td><td>200</td></tr><tr><td>5.3.2</td><td>Paramètres de simulation</td><td>201</td></tr><tr><td>5.3.3</td><td>Résultats numériques : l'écoulement turbulent</td><td>204</td></tr><tr><td>5.3.4</td><td>Résultats numériques : la dispersion des particules de pollen</td><td>212</td></tr><tr><td>5.3.5</td><td>Comparaison qualitative avec l'expérience</td><td>218</td></tr><tr><td>5.3.6</td><td>Discussion</td><td>220</td></tr><tr><td>5.4</td><td>Dispersion pour des dispositifs discontinus différents</td><td>221</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>5.4.1</td>
      <td>Présentation des dispositifs simulés</td>
      <td>221</td>
    </tr>
    <tr>
      <td>5.4.2</td>
      <td>Résultats numériques: l'écoulement turbulent</td>
      <td>222</td>
    </tr>
    <tr>
      <td>5.4.3</td>
      <td>Résultats numériques: la dispersion des particules</td>
      <td>226</td>
    </tr>
    <tr>
      <td>5.4.4</td>
      <td>Résultats numériques: les vitesses de dépôt.</td>
      <td>232</td>
    </tr>
    <tr>
      <td>5.4.5</td>
      <td>Discussion</td>
      <td>236</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Conclusion</td>
      <td>237</td>
    </tr>
    <tr>
      <td colspan="2">Conclusion Générale et Perspectives</td>
      <td>239</td>
    </tr>
    <tr>
      <td colspan="2">Annexes</td>
      <td>247</td>
    </tr>
    <tr>
      <td colspan="2">A Les équations dans la couche limite atmosphérique</td>
      <td>247</td>
    </tr>
    <tr>
      <td>A.1</td>
      <td>Équations fondamentales en mécanique des fluides</td>
      <td>247</td>
    </tr>
    <tr>
      <td>A.1.1</td>
      <td>L'équation de conservation de masse</td>
      <td>247</td>
    </tr>
    <tr>
      <td>A.1.2</td>
      <td>L'équation d'état d'un gaz parfait</td>
      <td>247</td>
    </tr>
    <tr>
      <td>A.1.3</td>
      <td>L'équation de conservation de quantité de mouvement.</td>
      <td>247</td>
    </tr>
    <tr>
      <td>A.1.4</td>
      <td>L'équation de l'énergie</td>
      <td>247</td>
    </tr>
    <tr>
      <td>A.1.5</td>
      <td>La force de Coriolis</td>
      <td>248</td>
    </tr>
    <tr>
      <td>A.2</td>
      <td>Hypothèses de Boussinesq</td>
      <td>248</td>
    </tr>
    <tr>
      <td>A.2.1</td>
      <td>État de référence de l'atmosphère</td>
      <td>248</td>
    </tr>
    <tr>
      <td>A.2.2</td>
      <td>Énoncé des hypothèses de Boussinesq</td>
      <td>249</td>
    </tr>
    <tr>
      <td>A.2.3</td>
      <td>Simplification de Boussinesq</td>
      <td>250</td>
    </tr>
    <tr>
      <td>A.3</td>
      <td>Mise en place des équations du mouvement turbulent</td>
      <td>251</td>
    </tr>
    <tr>
      <td>A.4</td>
      <td>Équations de l'énergie cinétique moyenne du mouvement turbulent et de sa dissipation</td>
      <td>253</td>
    </tr>
    <tr>
      <td>A.5</td>
      <td>Écoulement dans un couvert végétal.</td>
      <td>256</td>
    </tr>
    <tr>
      <td>A.5.1</td>
      <td>Équation de conservation de quantité de mouvement</td>
      <td>256</td>
    </tr>
    <tr>
      <td>A.5.2</td>
      <td>Équations de l'énergie cinétique moyenne du mouvement turbulent.</td>
      <td>259</td>
    </tr>
    <tr>
      <td colspan="2">B Compléments sur les particules biotiques</td>
      <td>261</td>
    </tr>
    <tr>
      <td>B.1</td>
      <td>Estimation de la vitesse de sédimentation en régime de Stokes</td>
      <td>261</td>
    </tr>
    <tr>
      <td>B.2</td>
      <td>Notion de temps de relaxation de la particule</td>
      <td>263</td>
    </tr>
    <tr>
      <td>B.3</td>
      <td>Glossaire</td>
      <td>265</td>
    </tr>
    <tr>
      <td colspan="2">Bibliographie</td>
      <td>271</td>
    </tr>
  </tbody>
</table>
# Nomenclature

## Lettres Latines

A : densité de surface végétale .................................... m².m⁻³

a : densité de surface végétale perpendiculaire à l'écoulement .... m².m⁻³

C : concentration instantanée, $\overline{C}$ moyenne temporelle, $\langle \overline{C} \rangle$ moyenne spatio-temporelle .................................... particule.s⁻³

C. : constantes du modèle de turbulence $k-l$, $C_{l1}$ $C_{l2}$ ............ adim.

C. : constantes du modèle de turbulence $k-\epsilon$ standard, $C_\mu$ $C_{\epsilon1}$ $C_{\epsilon2}$ adim.

C. : constantes supplémentaires du modèle de turbulence $k-\epsilon$ de canopée, $C_{Pw}$ $C_{\epsilon w}$ $C_{P\epsilon w}$ $C_{D\epsilon w}$ .................................... adim.

C_D : coefficient de traînée de la particule solide .................... adim.

C_d : coefficient de traînée du couvert végétal .................... adim.

C_P : capacité calorifique massique à pression constante ....... m².s⁻².K⁻¹

C_V : capacité calorifique massique à volume constant ......... m².s⁻².K⁻¹

D_\epsilon w : dissipation de sillage de $\epsilon$ ........................ m².s⁻⁴

d : hauteur de déplacement .................................... m

d_p : diamètre moyen de la particule .................................... m

E, $E_\infty$: paramètre pour la vitesse sur une plaque rugueuse ............ adim.

E_i : efficacité d'impact .................................... adim.
$E_s$ : efficacité de sédimentation .................................... adim.

$E_0$ : paramètre pour la vitesse sur une plaque lisse .................... adim.

$F$ : force volumique de traînée .................................... kg.m$^{-2}$.s$^{-2}$

$F_D$ : force de traînée de Newton et de Stokes ........................ kg.m.s$^{-2}$

$F_d$ : flux de particules déposées ................................ particule.m$^{-2}$.s$^{-1}$

$F_{di}$ : composante suivant i de la force volumique totale exercée sur le couvert végétal .................................... kg.m$^{-2}$.s$^{-2}$

$f$ : grandeur aléatoire associée au mouvement turbulent du fluide ........

$f_{Fi}$ : composante suivant i de la force volumique de la traînée de forme .... .................................... kg.m$^{-2}$.s$^{-2}$

$f_M$ : fonction universelle de $\xi$ .................................... adim.

$f_{Vi}$ : composante suivant i de la force volumique de la traînée de frottement. .................................... kg.m$^{-2}$.s$^{-2}$

$f_x, f_z$ : fraction de la densité de surface végétale suivant $z$ et $x$ ........ adim.

$g$ : accélération de la pesanteur .................................... m.s$^{-2}$

$H$ : hauteur de la source ou de la couche limite interne .................... m

$H^+$ : nombre de Reynolds de rugosité .................................... adim.

$h$ : hauteur moyenne du couvert végétal ou de la sous-couche interne .. m

$h_s$ : hauteur moyenne des aspérités de surface ........................ m

$k$ : énergie cinétique moyenne du mouvement turbulent ............ m$^2$.s$^{-2}$

$k_t$ : diffusivité thermique turbulente .................................... m$^2$.s$^{-2}$

$k_\phi$ : diffusivité moléculaire associée à la grandeur $\phi$ .................... m$^2$.s$^{-1}$

$k_{\phi t}$ : diffusivité turbulente associée à la grandeur $\phi$ .................... m$^2$.s$^{-1}$

$L$ : longueur de Monin et Obukhov .................................... m

$L_c$ : longueur de la clairière .................................... m
$L_s$ : échelle de cisaillement pour un couvert végétal .................... m

$L_E$ : macro-échelle spatiale de Taylor .................................... m

$l$ : longueur caractéristique de l'obstacle .................................... m

$l_m$ : longueur de mélange .................................................... m

$l'$ : échelle de longueur représentative du mouvement turbulent ............ m

$m$ : masse ................................................................. kg

$m$ : paramètre de rugosité .................................................... adim.

$P$ : pression instantanée, $\overline{P}$ moyenne temporelle, $\langle \overline{P} \rangle$ moyenne spatio-temporelle .................................................... kg.m$^{-1}$.s$^{-2}$

$P$ : paramètre pour la température sur une plaque rugueuse ............ adim.

$P_e$ : nombre de Peclet .................................................... adim.

$P_h$ : production gravitationnelle de $k$ .................................... m$^2$.s$^{-3}$

$P_k$ : production cinétique de $k$ .................................................... m$^2$.s$^{-3}$

$P_r$ : nombre de Prandtl .................................................... adim.

$P_{rt}$ : nombre de Prandtl turbulent $P_{rt} = \nu_t/k_t$ .................... adim.

$P_w$ : production de sillage de $k$ .................................................... m$^2$.s$^{-3}$

$P_{ew}$ : production de sillage de $\varepsilon$ .................................................... m$^2$.s$^{-4}$

$P_0$ : paramètre pour la température sur une plaque lisse ............ adim.

$Q$ : débit de la source .................................................... particule.s$^{-1}$

$Q_0$ : densité de flux de chaleur sensible .................................... kg.s$^{-3}$

$R$ : constante des gaz parfaits .................................................... adim.

$R_L$ : fonction d'autocorrélation lagrangienne .................................... adim.

$r$ : cadence de dépôt dans le couvert végétal .................................... s$^{-1}$

$R_a$ : nombre de Rayleigh .................................................... adim.
$R_e$ : nombre de Reynolds .................................................... adim.

$R_{ep}$ : nombre de Reynolds particulaire .................................... adim.

$R_f$ : nombre de Richardson .................................................... adim.

$R_i$ : nombre de Richardson de gradient .................................... adim.

$R_{ib}$ : nombre de Richardson volumique .................................... adim.

$S_p$ : terme source ou puits de particules .................... particule.m$^{-3}$.s$^{-1}$

$\langle \overline{S_p} \rangle$ : dépôt de particules par impact et par sédimentation dans le couvert végétal $\langle \overline{S_p} \rangle = \langle \overline{S_i} \rangle + \langle \overline{S_s} \rangle$ .................... particule.m$^{-3}$.s$^{-1}$

$S_t$ : nombre de Stokes .................................................... adim.

$T$ : température instantanée, $\overline{T}$ moyenne temporelle .................... K

$T_L$ : échelle de temps intégrale lagrangienne .................... s

$T_*$ : température de frottement .................................................... K

$t$ : temps .................................................... s

$u'$ : échelle de vitesse représentative du mouvement turbulent ...... m.s$^{-1}$

$u_i$ : composante de la vitesse instantanée suivant i .................... m.s$^{-1}$

$\overline{u}_i$ : composante de la vitesse moyenne temporelle suivant i ...... m.s$^{-1}$

$u'_i$ : fluctuation de la composante de la vitesse suivant i, écart à la vitesse moyenne temporelle .................................................... m.s$^{-1}$

$\langle \overline{u}_i \rangle$ : composante de la vitesse moyenne spatio-temporelle suivant i .. m.s$^{-1}$

$\overline{u}_i''$ : fluctuation de la composante de la vitesse suivant i, écart à la vitesse moyenne spatio-temporelle .................................................... m.s$^{-1}$

$u_*$ : vitesse de frottement .................................................... m.s$^{-1}$

$u_{*t}$ : vitesse de frottement critique .................................................... m.s$^{-1}$

$V$ : volume .................................................... m$^3$

$v_d$ : vitesse de dépôt .................................................... m.s$^{-1}$
<table><tr><td>v<sub>f</sub></td><td>vitesse de fuite</td><td>m.s<sup>-1</sup></td></tr><tr><td>v<sub>s</sub></td><td>vitesse de chute ou de sédimentation</td><td>m.s<sup>-1</sup></td></tr><tr><td>x</td><td>distance horizontale parcourue par une particule</td><td>m</td></tr><tr><td>x<sub>i</sub></td><td>coordonnée cartésienne suivant la direction i x,y,z</td><td>m</td></tr><tr><td>z</td><td>altitude</td><td>m</td></tr><tr><td>z<sub>*</sub></td><td>longueur de frottement</td><td>m</td></tr><tr><td>z<sub>0</sub></td><td>longueur de rugosité dynamique</td><td>m</td></tr><tr><td>z<sub>0h</sub></td><td>longueur de rugosité thermique</td><td>m</td></tr></table>

## Lettres grecques

α : coefficient d'atténuation du vent dans le couvert végétal ou facteur de forme dynamique de la particule adim.

β : coefficient de flottabilité m.s<sup>-2</sup>.K<sup>-1</sup>

ε : taux de dissipation de k m<sup>2</sup>.s<sup>-3</sup>

ε<sub>w</sub> : dissipation de sillage de k m<sup>2</sup>.s<sup>-3</sup>

η : micro-échelle de longueur de Kolmogorov m

θ : température potentielle K

κ : constante de Von Karman adim.

λ : conductivité thermique kg.m.s<sup>-3</sup>.K<sup>-1</sup>

μ : viscosité dynamique kg.m<sup>-1</sup>.s<sup>-1</sup>

μ<sub>t</sub> : viscosité turbulente kg.m<sup>-1</sup>.s<sup>-1</sup>

ν : viscosité cinématique ν = μ/ρ m<sup>2</sup>.s<sup>-1</sup>

ν<sub>t</sub> : viscosité turbulente ν<sub>t</sub> = μ<sub>t</sub>/ρ m<sup>2</sup>.s<sup>-1</sup>
$\xi$ : indice de Monin-Obukhov ou nombre aléatoire .................... adim.

$\rho$ : masse volumique ................................................ kg.m$^{-3}$

$\sigma_\phi$ : écart type dans les modèles gaussiens suivant $x, y, z$ .................... m

$\sigma_\phi$ : nombre de Prandtl-Schmidt turbulent associé à la variable $\phi$ \\
$\sigma_\phi = \nu_t/k_{\phi t}$ .................................................... adim.

$\tau$ : échelle de temps caractéristique .................................... s

$\tau_p$ : temps de relaxation des particules .................................... s

$\tau_0$ : contrainte de cisaillement à la surface .................... kg.m$^{-1}$.s$^{-2}$

$\Phi_x$ : flux de particules convectées par l'écoulement moyen, intégré suivant le plan horizontal .................................... particule.s$^{-1}$

$\Phi_{aéro}$ : flux intégré de particules convectées par l'écoulement moyen, rapporté au débit de la source .................................... %

$\Phi_{cum}$ : flux cumulé de particules déposées .................... particule.s$^{-1}$

$\Phi_{dépôt}$ : flux de particules déposées rapporté au débit de la source .... % m$^{-2}$

$\Phi_{émis}$ : débit estimé de la source .................................... particule.s$^{-1}$

$\phi$ : grandeur aléatoire associée au mouvement turbulent du fluide ........

$\phi_h$ : fonction de similitude flux-gradient pour la température ........ adim.

$\phi_m$ : fonction de similitude flux-gradient pour la vitesse ........ adim.

$\phi_\varepsilon$ : fonction de similitude associée à $\varepsilon$ .................................... adim.

$\psi_h$ : fonction de similitude pour le profil de la température ........ adim.

$\psi_m$ : fonction de similitude pour le profil de la vitesse ........ adim.

$\Omega_i$ : composante du vecteur rotation terrestre suivant $i$ ........ rad. s$^{-1}$

$\omega$ : vitesse angulaire de rotation de la terre .................... rad. s$^{-1}$

## Vecteurs
F : force volumique ................................................ kg.m⁻².s⁻²

g : vecteur gravité ................................................ m.s⁻²

L : terme linéaire de l'équation de Navier-Stokes ...................... s⁻¹

n : normale sortante ................................................ adim.

r : vecteur position ................................................ m

u : vecteur vitesse ................................................ m.s⁻¹

v : vecteur vitesse instantané du fluide, vue par la particule solide . m.s⁻¹

v_p : vecteur vitesse de la particule solide ........................ m.s⁻¹

x : vecteur position ................................................ m

## Opérateurs mathématiques

$f = \bar{f} + f'$ : décomposition d'une variable instantanée suivant l'opérateur moyenne de Reynolds

$f = \langle \bar{f} \rangle + \bar{f}'' + f'$ : décomposition d'une variable instantanée suivant l'opérateur moyenne spatio-temporelle

$\langle \phi \rangle$ : moyenne spatiale $\langle \phi \rangle (\mathbf{x}, t) = \frac{1}{V} \iiint_V \phi (\mathbf{x} + \mathbf{r}, t) \, d\mathbf{r}$

$\frac{\partial}{\partial t}$ : dérivée partielle Eulérienne par rapport au temps

$\frac{\partial}{\partial x_j}$ : dérivée partielle Eulérienne dans la direction $j$

$\frac{d}{dt}$ : dérivée particulière $\frac{d}{dt} = \frac{\partial}{\partial t} + u_j \frac{\partial}{\partial x_j}$

$\delta_{ij}$ : symbole de Kronecker

$\epsilon_{ijk}$ : symbole d'Einstein

$\nabla$ : opérateur gradient

$\nabla$. : opérateur divergence
Indices et exposants

$\square_f$ : fluide

$\square_H$ : altitude du domaine d'étude

$\square_h$ : relatif à la thermique ou au couvert végétal

$\square_{i,j,k}$ : directions suivant $x, y, z$ ou 1, 2, 3

$\square_l$ : écart à l'état de référence, variable fonction de $x,y,z$

$\square_m$ : relatif à la dynamique

$\square_p$ : paroi ou particule

$\square_r$ : état de référence, variable fonction de $z$

$\square_s$ : sol

$\square_t$ : turbulent

$\square_w$ : sillage

$\square_0$ : valeur moyenne de l'état de référence

$\square_\infty$ : valeur à l'infinie

$\square_*$ : frottement

$\square'$ : fluctuation, écart à la moyenne temporelle

$\square''$ : fluctuation, écart à la moyenne spatiale de la moyenne temporelle
# Introduction Générale

La modélisation des écoulements diphasiques est de plus en plus nécessaire et trouve ses applications dans de vastes domaines : en industrie comme par exemple pour le traitement des eaux résiduaires, en génie chimique, mais aussi et surtout dans les écoulements dits naturels, comme les écoulements en rivière ou dans notre cas le transport des particules fines dans l'atmosphère. La présence de particules dans l'atmosphère est en effet un problème important car elle est à l'origine de la détérioration de la qualité de l'air.

De nombreux modèles de pollution atmosphérique visant à améliorer la qualité de l'air ont été élaborés ces dernières années. Ces modèles sont principalement orientés vers la pollution générée par l'activité humaine (circulation automobile, usines...). Paradoxalement, peu de modèles traitent de la pollution atmosphérique issue des particules biotiques telles que les spores, les pollens et les graines. Pourtant, il est maintenant établi qu'une part non négligeable de l'altération de la qualité de l'air est imputable à un grand nombre d'organismes biologiques provocant chez certaines personnes, des manifestations allergiques plus ou moins sévères telles que rhinites, conjonctivites, urticaire...

Si cette forme de pollution reste relativement bénigne, l'intérêt d'élaborer des modèles de dispersion atmosphérique appliqués aux particules biotiques n'en demeure pas mineur pour autant. En effet, le vent joue un rôle spécifique dans la dispersion des organismes biologiques, car cette dispersion altère nos pratiques agricoles et notre environnement. Tout d'abord, le vent est le vecteur principal des agents phytopathogènes, dont la dispersion est souvent à l'origine de l'intensification de foyers infectieux et de la propagation des épidé mies. Ces maladies plus ou moins virulentes peuvent occasionner des pertes économiques non négligeables, et les stratégies de prévention s'appuient souvent sur l'utilisation de produits de traitement chimique comme les pesticides. Enfin, le transport par le vent des particules biotiques entraîne une dissémination des gènes dans l'environnement, causant ainsi une évolution naturelle des espèces. Cette dissémination est également nommée flux de gène.
Dans ce contexte, vient s'ajouter depuis peu la problématique liée à l'introduction
volontaire d'organismes génétiquement modifiés dans l'environnement. En 1983, les bio-
technologies ont « conçu » une première génération de plantes tolérantes vis-à-vis d'un
agent spécifique (ravageur, herbicide...) en modifiant durablement son patrimoine géné-
tique par des techniques de transgénèse. Dès 1987, un maïs transgénique est cultivé pour la
première fois aux États-Unis. Depuis, cette innovation récente est en perpétuelle évolution
et sa banalisation par une culture à grande échelle est amplement sujette à controverse
qui s'inscrit dans une dimension sociale, économique, politique et scientifique. Parmi les
réflexions soulevées, la question des risques sur notre écosystème, liés au flux de gène entre
une culture de plantes transgéniques et une population de plantes sauvages apparentées
à la culture considérée, est récurrente car non évaluée et met en exergue la nécessité de
disposer d'outils pour quantifier ces risques.

<figure><img src="image_6.png" /><figcaption>Figure 1 : *Les évènements qui interviennent dans la dissémination des particules bio-tiques: phase de libération passive, phase de transport et phase de dépôt. Cycle « aérien » d'une épidémie causée par un phytopathogène [Aylor 90]*</figcaption></figure>

La dispersion atmosphérique des particules biotiques peut se décrire par trois évè-
nements importants: une phase de libération de la particule, son transport par le vent
dans l'atmosphère et enfin la phase de dépôt au sol ou dans un couvert végétal (figure 1).
Tous ces processus s'avèrent corrélés car ils dépendent fortement des facteurs climatiques
puisqu'ils ont lieu *in situ* et sont de plus, complexes à décrire car ils font intervenir di-
verses disciplines comme la biologie (structure d'ancrage des particules), la mécanique
des fluides et les transferts thermiques. Si de nombreuses expériences sont menées pour
modéliser et appréhender la dispersion atmosphérique des particules biotiques, elles sont
malheureusement peu transposables d'un site à l'autre. Toutes ces considérations rendent l'élaboration des modèles de dispersion difficiles, dont les processus requièrent une modélisation plus physique. Ces modélisations ont connu un essor significatif ces dernières années. Toutefois, elles s'avèrent souvent insuffisantes car les modèles dynamiques sont basés sur une paramétrisation du vent, de la turbulence et de la stabilité thermique. En effet, ces modèles n'intègrent pas des informations essentielles comme la structure du paysage qui est rarement homogène. Or la dispersion atmosphérique des particules biotiques est étroitement liée à l'écoulement de l'air dont les caractéristiques –intensité du vent moyen, turbulence générée par cisaillement ou d'origine thermique– sont tributaires de la variabilité spatiale du paysage, aussi bien à petite qu'à grande échelle. La prise en compte de l'influence du paysage et des conditions météorologiques qui en découlent passe inévitablement par des outils numériques basés sur la résolution des équations du mouvement du vent, de la turbulence, de la température.

C'est précisément l'objet de notre travail de thèse. Cette étude vise à développer un modèle numérique de dispersion atmosphérique de particules biotiques, basé sur la résolution des équations du mouvement qui gouvernent la couche limite atmosphérique. En particulier, le modèle doit contribuer à l'analyse de l'influence de la structure d'un paysage hétérogène présentant diverses transitions comme des transitions sol nu-végétation sur la dispersion des particules car chaque discontinuité a des effets sur le mouvement de l'air et sur sa capacité à diffuser et à transporter les particules. À terme, et pour répondre à une demande de plus en plus présente avec l'adoption éventuelle d'une politique de culture à grande échelle d'OGM, la modélisation de la dispersion atmosphérique, en étroite association avec des expérimentateurs de terrain, peut constituer une voie d'investigation pour évaluer et quantifier les potentialités de contamination dans un paysage agricole.

Il est évident que notre ambition n'est pas de réaliser un modèle complet, qui couvre l'ensemble des processus physiques qui interviennent, mais d'avancer pas à pas sur chacun de ces aspects en développant des méthodes adaptées aux problèmes qui nous préoccupent. Ce mémoire de thèse a donc pour objectif de présenter la démarche et l'approche menée pour l'élaboration de notre outil numérique, en précisant ses limites, ses possibilités et bien entendu les résultats obtenus. Ce mémoire s'articule autour de cinq chapitres :

- Le premier et le deuxième chapitre présentent l'écriture mathématique des principaux phénomènes physiques qui interviennent dans la dispersion atmosphérique des particules biotiques. Deux échelles d'observation des processus physiques ont été choisies pour modéliser et étudier la dispersion : l'échelle du paysage et l'échelle
du couvert végétal. Si l'échelle du paysage s'appuie sur la théorie de similitude de Monin et Obukhov pour reproduire les transferts turbulents qui s'effectuent entre la surface apparente et l'atmosphère, le modèle à l'échelle du couvert résout expliquement l'écoulement dans la sous-couche rugueuse. Cette étape a nécessité l'introduction d'un opérateur moyenne spatio-temporelle pour la description des variables turbulentes, ce qui a contribué à l'élaboration d'un modèle de canopée végétale;

- Le troisième chapitre présente quant à lui, les techniques numériques qui sont déjà disponibles dans l'outil de calcul choisi, ainsi que les modifications que nous avons apportées pour l'adapter à la dispersion atmosphérique des particules biotiques, à grande et à petite échelle;

- Le quatrième chapitre est axé sur les validations nécessaires du modèle numérique de dispersion atmosphérique. Nous avons reproduit numériquement plusieurs expériences issues de la littérature. Ces simulations visent principalement à vérifier que le modèle défini «ressent» les transitions de surfaces –saut de rugosité de surface et saut de température pariétale– à l'échelle du paysage. Pour le modèle de canopée, après une phase de validation très complète, nous nous sommes attachés à vérifier le bon comportement du modèle pour des transitions végétation-clairière-végétation;

- Le dernier chapitre est consacré à l'analyse et à la détermination des nombreux paramètres qui influencent la dispersion atmosphérique des particules biotiques. Pour permettre une meilleure compréhension des interactions entre la structure du paysage et la dispersion des particules biotiques, nous avons mené diverses simulations numériques en 2 et 3 dimensions pour des dispositifs discontinus maïs-sol nu. Ces dernières simulations constituent les premiers résultats pertinents et mettent en valeur l'intérêt de poursuivre ce travail.
# Chapitre 1

## Modélisation physique des phénomènes

### Introduction

L'objet de ce chapitre est de modéliser les phénomènes physiques qui ont lieu dans la couche limite atmosphérique pour définir un modèle dynamique et thermique.

En premier lieu, nous présenterons la structure de la couche limite atmosphérique et nous spécifierons les échelles qui nous intéressent pour cette étude. Nous avons choisi de travailler à l'échelle du paysage et à une échelle plus locale, à savoir l'échelle d'un couvert.

Que ce soit à l'échelle locale ou à plus grande échelle, nous utiliserons les équations fondamentales de la mécanique des fluides pour modéliser le vent. En raison du caractère turbulent de toutes les variables présentes au sein de cette couche, les équations nécessitent un traitement spécifique pour appréhender le mouvement turbulent. Nous choisirons des méthodes statistiques basées sur une description des propriétés moyennes de l'écoulement. Nous verrons que ces traitements sont différents suivant l'échelle d'observation des processus. Ces traitements statistiques entraînant des problèmes de fermeture, nous proposerons des modèles de fermeture à une ou à deux équations d'évolution.

Les transferts turbulents entre les surfaces naturelles et l'atmosphère sont de types convectifs et sont influencés par la stabilité de l'atmosphère. À l'échelle du paysage, la stabilité thermique sera traitée par la théorie de similitude de Monin et Obukhov définie dans la couche de surface.
<figure><img src="image_7.png" /></figure>
## 1.1 La couche limite atmosphérique

Le sol terrestre a une grande influence sur les mouvements de la basse atmosphère et sur la nature même de la turbulence. L'air étant un fluide visqueux, sa vitesse s'annule au contact du sol. Près de cette zone, les gradients (de vitesse, de température ...) sont intenses et l'air est soumis à de fortes agitations. Cette partie de l'atmosphère s'appelle la couche limite atmosphérique : elle est définie, d'un point de vue dynamique, comme étant la zone de l'atmosphère où le comportement de l'air est directement influencé par le sol. D'un point de vue thermique, c'est la zone de l'atmosphère au voisinage de la surface terrestre dans laquelle la variation diurne du rayonnement solaire est directement perceptible.

La couche limite atmosphérique est fortement turbulente. En effet, si nous considérons une vitesse caractéristique du vent de $u = 10 \, m \, s^{-1}$, que la viscosité cinématique de l'air est de l'ordre $\nu = 10^{-5} \, m^2 \, s^{-1}$ et l'épaisseur de la couche atmosphérique $h_{CLA} = 1000 \, m$, le nombre de Reynolds caractéristique de la couche limite atmosphérique est largement supérieur à la valeur critique de transition d'un écoulement laminaire vers un écoulement turbulent :

$$R_e = \frac{uh_{CLA}}{\nu} = 10^9 \gg R_{crit} \approx 3000$$

De même, les effets thermiques, caractérisés par le nombre de Rayleigh, sont tels que le nombre de Rayleigh de la couche limite atmosphérique est bien supérieur à sa valeur critique. Un rapide calcul permet d'évaluer ce nombre sans dimension en considérant que le coefficient de diffusivité thermique turbulente est $k_t = 2.10^{-5} \, m^2 \, s^{-1}$ que la gravité est $g = 10 \, m \, s^{-2}$ pour une température moyenne de $T = 300 \, K$. Une variation de 1 degré conduit à un nombre de Rayleigh de l'ordre de :

$$R_a = \frac{g h_{CLA}^3 \Delta T}{\nu k_t} = 1.510^{17} \gg R_{crit} \approx 50000$$

### 1.1.1 Structure de la couche limite atmosphérique

La couche limite atmosphérique a une épaisseur $h_{CLA}$ pouvant varier de quelques mètres à quelques kilomètres [Pielke 84], suivant que l'atmosphère est en stratification stable, instable ou neutre. Dans cette couche, il se distingue plusieurs zones, possédant chacunes ses propres caractéristiques :

- La couche limite superficielle ou couche de surface est en contact avec le sol. Son
épaisseur $h_s$ varie typiquement entre 10 m à 100 m [Pielke 84]. Au sein de cette
couche, il est possible d'identifier deux zones : la sous-couche rugueuse et la sous-
couche inertielle. Dans la sous-couche rugueuse, l'écoulement du vent dépend di-
rectement des variables locales de la surface (rugosité, petite topographie, couvert
végétal...). Ces obstacles génèrent de nombreux tourbillons caractérisés par une tur-
bulence très intermittente et instationnaire. C'est dans cette couche qu'ont lieu les
échanges (énergie, quantité de mouvement, polluants...) entre le sol et l'atmosphère.
Juste au-dessus suit la sous-couche inertielle caractérisée par un écoulement de vent
quasi-stationnaire où les forces de Coriolis sont négligeables, ce qui a pour consé-
quence une direction du vent constante avec l'altitude. Les flux verticaux de quan-
tité de mouvement, de masse et de chaleur sont constants. Les caractéristiques de
la surface n'influencent l'écoulement qu'au travers de paramètres macroscopiques :
ces paramètres sont appelés longueurs de rugosité dynamique et thermique. Du fait
de son accessibilité aux mesures, la couche de surface reste de loin la plus connue;

- La couche d'Ekman (ou la couche de transition) correspond à la partie supérieure de la couche limite atmosphérique. Dans cette zone, comprise entre $h_s$ et $h_{CLA}$, la structure du vent est influencée à la fois par le frottement de l'air à la surface (couche bien mélangée par la turbulence à grande échelle) et par les forces de Coriolis (apparition de spirales d'Ekman). Avec l'altitude, les effets du frottement au sol deviennent négligeables devant la force de Coriolis et la direction du vent subit des déviations angulaires (pouvant atteindre 30 à 40°) pour atteint à son sommet celle du vent géostrophique.

Au sommet de la couche limite, il apparaît une couche d'entraînement et/ou une inversion thermique caractérisée par une stabilité thermique, qui bloque les ascensions d'air, et par une brusque variation de la température et du vent. Enfin, à la couche limite succède l'atmosphère libre où l'écoulement est régi uniquement par les forces de pression et de Coriolis et dans laquelle les champs météorologiques sont déterminés par la dynamique à grande échelle.

1.1.2 Problème d’échelles. Définition des échelles de l’étude

Les mouvements de l'air au sein de l'atmosphère s'étendent sur une très large gamme
d'échelles spatiales et temporelles (figure 1.1), échelles qui identifient leur extension hori-
zontale ainsi que leur durée de vie.
<table><caption>Figure 1.1 : *Les diverses échelles de mouvement atmosphérique et les phénomènes météo-rologiques associés avec leurs échelles de temps et d'espace* [Thunis 96]</caption><thead><tr><th>L<sub>H</sub></th><th>Lifetime</th><th>Stull (1988)</th><th>Pielke (1984)</th><th>Orlanski (1975)</th><th>Present</th><th>Atmospheric phenomena</th></tr></thead><tbody><tr><td rowspan="2">10 000 km</td><td rowspan="2">1 month</td><td rowspan="2"></td><td>S<br>y<br>n<br>o<br>p<br>t<br>i<br>c</td><td>Macro-α</td><td>Macro-α</td><td>General circulation, long waves</td></tr><tr><td>R<br>e<br>g<br>i<br>o<br>n<br>a<br>l</td><td>Macro-β</td><td>Macro-β</td><td>Synoptic cyclones</td></tr><tr><td>2000 km</td><td>1 week</td><td>M<br>a<br>c<br>r<br>o</td><td></td><td>Meso-α</td><td>Macro-γ</td><td>Fronts, hurricanes</td></tr><tr><td rowspan="2">200 km</td><td rowspan="2">1 day</td><td rowspan="2"></td><td rowspan="2"></td><td>Meso-β</td><td>Meso-β</td><td>Low-level jets, thunderstorm groups, mountain winds and waves, sea breeze, urban circulations</td></tr><tr><td>Meso-γ</td><td>Meso-γ</td><td>Thunderstorm, clear-air turbulence</td></tr><tr><td>20 km</td><td rowspan="2">1 h</td><td rowspan="2"></td><td rowspan="2"></td><td>Micro-α</td><td>Meso-δ</td><td>Cumulus, tornadoes, katabatic jumps</td></tr><tr><td>2 km</td><td>Micro-β</td><td>Micro-β</td><td>Plumes, wakes, waterspouts, dust devils</td></tr><tr><td rowspan="2">200 m</td><td rowspan="2">30 min</td><td rowspan="2">M<br>i<br>c<br>r<br>o</td><td rowspan="2"></td><td rowspan="2">Micro-γ</td><td>Micro-γ</td><td rowspan="2">Turbulence, sound waves</td></tr><tr><td>Micro-δ</td><td></td></tr><tr><td>20 m</td><td>1 min</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2 m</td><td>1 s</td><td>M<br>i<br>c<br>r<br>o<br>δ</td><td></td><td></td><td></td><td></td></tr></tbody></table>
Les plus grands mouvements atmosphériques définissent l'échelle climatologique, planétaire, synoptique, supérieure à la centaine de kilomètres et d'une durée de vie de plus de 1 jour mais qui peut atteindre un mois. La turbulence associée à ces mouvements est, compte tenu de l'étendue verticale limitée de l'atmosphère, quasi bidimensionnelle : elle diffère donc de la turbulence tridimensionnelle « usuelle » et constitue un domaine privilégié pour de nombreuses études.

Les plus fines échelles des mouvements atmosphériques, imposées par les structures dissipatives, ont une taille inférieure au kilomètre avec une durée de vie de quelques minutes et forment la micro-échelle. La turbulence est tridimensionnelle et l'énergie cinétique est transférée depuis les échelles les plus grandes vers les plus petites échelles de la dissipation. Au cours de ce processus de transfert énergétique, il y a une disparition progressive de la direction privilégiée des gros tourbillons, et à un certain stade, la turbulence devient localement isotrope. C'est le processus de cascade décrit par Kolmogorov (1941).

Enfin, entre ces deux échelles extrêmes, les mouvements de tailles intermédiaires définissent la méso-échelle et assurent la transition. Notons qu'il se distingue également, les échelles sub-meso définies entre 10 mètres à 10 kilomètres [Costes 96].

<figure><img src="image_8.png" /><figcaption>Figure 1.2 : Échelle du paysage (domaine méso-δ)</figcaption></figure>

Les processus atmosphériques couvrent donc une large gamme d'échelles spatiales et temporelles, il en ressort qu'il n'est pas possible de décrire l'ensemble de ces processus par une seule modélisation. L'objectif de cette étude étant de simuler le transport de particules biotiques, la distance de dispersion caractéristique impose directement l'échelle d'observation des processus. Nous avons choisi de travailler sur deux échelles :

- La dispersion à l'échelle sub-méso qui correspond aux écoulements à l'échelle d'un paysage (figure 1.2). À cette échelle, il est important de tenir compte de la structure
du paysage qui est rarement homogène. En effet, les discontinuités comme des brisements ou des couverts végétaux modifient les échanges convectifs dans la couche de surface. Chaque surface provoque un forçage local sur l’écoulement, dont les caractéristiques sont conservées au-delà de la surface. C’est là l’advection locale ou advection de petite échelle [Brunet 97]. Il en résulte que les flux turbulents sont très différents dans la zone de transition, et peuvent considérablement s’écarter de leurs valeurs amont et aval. Par suite, ces effets peuvent modifier la dispersion de polluants, de particules. Un modèle à cette échelle intègre ses informations en décrivant le paysage comme une mosaïque d’unités élémentaires, chaque unité étant considérée comme homogène caractérisée par plusieurs paramètres moyens tels que la rugosité, la hauteur de déplacement;

- La dispersion à courte distance. Cette modélisation s’appuie sur une description du vent et de la turbulence à une échelle plus locale, telle que l’échelle d’un couvert végétal ou une parcelle homogène. À cette échelle, il est indispensable de prendre en compte de manière explicite la sous-couche rugueuse. Cette étape est d’une part nécessaire car la dispersion de particules biotiques est en fait un processus d’échange entre l’atmosphère et la végétation, d’autre part, la modélisation à l’échelle du paysage ne peut être acquise qu’en approfondissant ces mécanismes élémentaires. Compte tenu de la spécificité de la turbulence dans le couvert, cette partie (pour des raisons de simplicité) sera traitée en dernière section de ce chapitre (section 1.6).

### 1.1.3 Les équations instantanées simplifiées

Afin de modéliser le vent dans la couche limite atmosphérique, en se limitant au cas de l'air sec, nous avons choisi d'utiliser les équations fondamentales de la mécanique des fluides et de la thermique<fnref n="1" />, l'air étant assimilé à un fluide newtonien se comportant comme un gaz parfait. Les variables qui interviennent dans la description de l'écoulement dans la couche limite atmosphérique sont :

- Les composantes de la vitesse $u_i$;

- Les variables thermodynamiques : la pression $P$, la température $T$ et la masse volumique $\rho$.

<footnote n="1">1. Équations en annexe A.1</footnote>
Le système d'équations de la couche limite atmosphérique s'avère complexe sous sa forme primaire et peu exploitable. Cependant, l'approximation de Boussinesq <fnref n="2" /> (1903) permet de définir un système d'équations simplifié. Cette approximation considère que l'état thermodynamique de l'atmosphère s'écarte peu d'un état hydrostatique et adiabatique, considéré comme un état de référence connu : ainsi les variables thermodynamiques sont entièrement déterminées. Toute variable thermodynamique admet une décomposition en une variable indicée *r*, qui est l'état de référence évalué au repos et une variable indicée *l*, qui est l'écart à cet état et sur laquelle nous travaillons directement.

Dans le cadre des hypothèses de Boussinesq <fnref n="3" />, les équations instantanées retenues qui gouvernent la couche limite atmosphérique sont :

- L'équation instantanée de conservation de masse sous forme incompressible :

$$ \frac{\partial u_j}{\partial x_j} = 0 \qquad (1.1) $$

- L'équation d'état :

$$ \rho_l = -\rho_r \frac{T_l}{T_r} \qquad (1.2) $$

- L'équation instantanée de conservation de quantité de mouvement :

$$ \rho_r \frac{\partial u_i}{\partial t} + \rho_r u_j \frac{\partial u_i}{\partial x_j} = -\rho_l g \delta_{i3} - \frac{\partial P_l}{\partial x_i} + \frac{\partial}{\partial x_j} \left[ \mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \right] - 2 \rho_r \epsilon_{ijk} \Omega_j u_k \quad (1.3) $$

- L'équation instantanée de la conservation de l'énergie :

$$ \rho_r C_P \left[ \frac{\partial T_l}{\partial t} + u_j \frac{\partial T_l}{\partial x_j} \right] = \frac{\partial}{\partial x_j} \left( \lambda \frac{\partial T_l}{\partial x_j} \right) \qquad (1.4) $$

Pour les écoulements atmosphériques, il est d'usage d'introduire la température potentielle $\theta$, définie comme la température que prendrait une particule à l'état $(P, T)$ si elle était amenée adiabatiquement à la pression $P_{00} = 10^5$ Pa. Cette température, inviante au cours d'une transformation adiabatique, sera transportée lors des mouvements verticaux de l'atmosphère. Pour un gaz parfait, elle s'exprime :

$$ \theta = T \left( \frac{P}{P_{00}} \right)^{-R/C_P} \qquad (1.5) $$

<footnote n="2">2. Énoncé en annexe A.2.2</footnote>

<footnote n="3">3. Calculs présentés en annexe A.2.3</footnote>
La température potentielle est simplement reliée à la température statique [DeMoor 83]
$T_l$ par :

$$
\theta = T_l + T_r(h_0) \tag{1.6}
$$

Dans cette étude, les échelles mises en jeu nous permettent de considérer en première approximation, qu’il est équivalent de travailler sur la température statique $T_l$.

Ce système d'équations possède une solution régulière (en considérant que les conditions aux limites sont connues) tant que les paramètres caractéristiques de l'écoulement, tel que le nombre de Reynolds et le nombre de Rayleigh restent inférieurs à des valeurs critiques. Or nous avons vu que ces paramètres sont largement supérieurs et que la couche limite atmosphérique est fortement turbulente. Du fait de la nature turbulente de l'écoulement, ces équations présentent des solutions fortement instables et nécessitent des traitements spécifiques à l'approche utilisée pour la prise en compte du mouvement turbulent.

1.2 Prise en compte du mouvement turbulent

Il n'existe pas encore aujourd'hui de théorie générale et complète de la turbulence. À défaut de donner donc une définition rigoureuse, il est possible d'indiquer un ensemble de caractéristiques dont la présence simultanée indique que l'écoulement considéré est turbulent :

- Les variables $u_i, P, T, \dots$ décrivant un écoulement turbulent présentent un caractère aléatoire en fonction du temps et de l'espace. Cette caractéristique laisse entrevoir l'utilisation naturelle des méthodes statistiques;

- Les écoulements turbulents présentent un caractère rotationnel : des tourbillons de toutes tailles -les plus gros sont de dimension comparable à celle de l'écoulement, les plus petits correspondent à des nombres de Reynolds proche de l'unité- peuvent être observés. Cette large gamme d'échelle fait que la turbulence reste un phénomène continu ;

- La turbulence est un phénomène fondamentalement non linéaire, intrinsèquement liée aux termes d'inertie dans les équations de Navier qui sont responsables du mécanisme de cascade de l'énergie cinétique du mouvement turbulent qui est produite aux grandes échelles et dissipée en chaleur par viscosité aux petites échelles ;
- La turbulence est un phénomène dissipatif, c'est-à-dire que les processus irréversibles, tels que la dissipation de l'énergie cinétique en chaleur, sont considérablement accrus;

- Enfin, la propriété considérée comme la plus importante est l'aptitude d'un écoulement turbulent à diffuser toute grandeur associée au fluide, avec une efficacité dépassant de plusieurs ordres la diffusion moléculaire.

1.2.1 Quelques approches pour la modélisation de la turbulence

Les caractéristiques des écoulements turbulents sont fondamentalement différentes des
écoulements laminaires. Le rôle de la turbulence au sein de l'atmosphère est essentiel
puisqu'elle régule le fonctionnement de la machine thermique que représente la terre. Elle
constitue l'un des facteurs prédominants, non seulement des mouvements mais encore des
répartitions des propriétés physiques et chimiques dans tout l'espace aérien. Compte tenu
du rôle et des effets de la turbulence dans l'atmosphère, il est capital de tenir compte
de tous ces phénomènes. Il n'existe pas encore de théorie générale pour appréhender
le mouvement turbulent, mais il existe néanmoins différentes approches qui permettent
de résoudre les équations que nous venons de voir. Parmi ces méthodes, nous pouvons
distinguer les méthodes de simulation instationnaire des échelles de la turbulence :

- La simulation numérique directe (DNS) utilise le système d'équations présenté sans aucun traitement spécifique, mais doit « capter » l'écoulement dans sa totalité : cette simulation déterministe doit résoudre toutes les échelles spatiales et temporelles du mouvement. Il est alors nécessaire de disposer d'un maillage de l'ordre de $R_e^{3/4}$ dans chaque direction (soit un total de $10^{81/4}$ points pour un nombre de $R_e = 10^9$) afin de traiter les plus petites échelles dissipatives de la turbulence. Ces contraintes dépassent largement les moyens informatiques actuels pour la modélisation des écoulements atmosphériques ;

- La simulation directe des grandes échelles (LES) permet grâce à une technique de filtrage spatio-temporelle [Sergent 00] de résoudre les grandes échelles de la turbulence et modélise les structures inférieures en extrayant artificiellement de l'énergie pour les échelles dissipatives. Ces méthodes sont en plein essor et connaissent beaucoup de succès (un module de LES est en cours de développement au laboratoire), mais
imposent une bonne connaissance des processus qui interviennent à petite échelle, dont la modélisation reste semi-empirique.

Une approche très différente consiste à modéliser le mouvement turbulent à l'aide de méthodes statistiques que nous détaillerons à la section suivante. Cette approche, qui sera la nôtre, est largement utilisée [Duynkerke 88, Detering 85] et semble bien adaptée pour la simulation des écoulements atmosphériques à grands nombres de Reynolds, pour une turbulence pleinement développée.

### 1.2.2 Description statistique de la turbulence

Hypothèse d'ergodicité: la moyenne dans le temps d'un paramètre quelconque est égale à la moyenne de ce paramètre prise sur un ensemble de systèmes [Ngô 95]. Ce postulat suppose implicitement l'équiprobabilité de micro-états accessibles à un système isolé en équilibre, et que les variables sont statistiquement stationnaires.

Soit un ensemble de N systèmes au temps t. La moyenne de la variable f au temps t, prise sur l'ensemble des N systèmes est définie par :

$$ \langle f \rangle_N = \frac{1}{N} \sum_{i=1}^{N} f(t) \qquad (1.7) $$

La moyenne dans le temps de la même variable f considérée comme statistiquement stationnaire est quant à elle :

$$ \bar{f} = \lim_{t' \to \infty} \frac{1}{t'} \int_0^{t'} f(t') dt' \qquad (1.8) $$

Le postulat d'ergodicité s'écrit :

$$ \bar{f} = \langle f \rangle_N \qquad (1.9) $$

L'intérêt du postulat d'ergodicité apparaît de façon évidente: pour définir une loi de probabilité du phénomène aléatoire étudié, il faudrait considérer N expériences identiques et faire tendre N vers l'infini, ce qui est irréalisable. L'hypothèse énoncée stipule qu'il est équivalent de prendre N échantillons sur une même expérience, et de faire une moyenne sur un temps T assez long.

Devant l'étendue des échelles spatio-temporelles des écoulements atmosphériques, dont nous avons vu que la structure interne de la turbulence peut s'avérer très différente, il se pose une difficulté supplémentaire dans l'approche statistique. En effet, la définition
<figure><img src="image_9.png" /><figcaption>Figure 1.3 : *Spectre d'énergie des mouvements atmosphériques d'après Van der Hoven*
(1957)</figcaption></figure>

d'une moyenne pose des problèmes, et non des moindre, car elle sous-entend la possi-
bilité de dissocier les mouvements atmosphériques « moyens » ou « macroscopiques » du
mouvement qui constitue la turbulence usuelle, dont les caractéristiques nous intéressent
pour cette étude. La solution empirique, à ce problème délicat, repose sur l'analyse du
spectre de la vitesse de vent horizontale donnée par Van der Hoven (1957) (figure 1.3). Ce
spectre présente trois évènements remarquables : plusieurs pics d'énergie pour des basses
fréquences autour de 1 an (variation saisonnière) 4 jours (variation climatique) et 1 jour
(variation journalière) qui sont clairement identifiés aux mouvements atmosphériques de
large amplitude.

Le spectre présente ensuite un minimum accentué pour des périodes comprises entre
quelques minutes et quelques heures : c'est le trou spectral (spectral gap) [Coantic 76,
DeMoor 83]. Enfin, le spectre d'énergie a un maximum aux hautes fréquences, autour de
la minute. Ces périodes sont celles des mouvements à micro-échelle caractérisés par une
turbulence générée par le frottement au sol (rugosité, obstacle) et constituent la turbu-
lence proprement dite, celle qui se rencontre usuellement en laboratoire. En s'appuyant
sur l'analyse de ce spectre, il en ressort qu'il est possible de dissocier les mouvements
atmosphériques à grandes échelles des mouvements à petites échelles puisque le spectre
d'énergie présente un « trou », en opérant des moyennes sur des durées de temps de l'ordre
de quelques dizaines de minutes.

1.2.3 Traitement statistique de Reynolds

Ainsi l'utilisation de l'opérateur moyenne temporelle permet de décrire le phénomène
aléatoire qu'est la turbulence. C'est en s'appuyant sur les hypothèses émises plus haut,
que Reynolds (1878) s'est intéressé non plus aux grandeurs instantanées d'un écoulement
turbulent, mais à ses propriétés moyennes.

Soit $f$ une fonction aléatoire qui représente une grandeur associée au mouvement tur-
bulent du fluide. Le traitement statistique introduit par Reynolds consiste à décomposer
$f$ en une partie moyenne temporelle notée $\bar{f}$ et une partie fluctuante $f'$, qui est l'écart à
cette moyenne. Le traitement statistique de Reynolds s'écrit donc :

$$
f = \bar{f} + f' \tag{1.10}
$$

Reynolds établit les règles de calcul suivantes (g est une variable aléatoire et α une variable non aléatoire) :

$$
\begin{align*}
\overline{f+g} &= \overline{f} + \overline{g} \\
\overline{\bar{f}} &= \overline{f} \\
\overline{\frac{\partial f}{\partial x_i}} &= \overline{\frac{\partial \bar{f}}{\partial x_i}} \\
\overline{\alpha f} &= \alpha \overline{f} \\
\overline{f'} &= 0 \\
\overline{fg} &= \overline{f\bar{g}} + \overline{f'g'} \tag{1.11}
\end{align*}
$$

Une des conséquence inhérente à l'approche statistique est un accroissement du nombre d'inconnues : la turbulence étant un processus non-linéaire, il en résulte que le traitement statistique du produit de deux variables aléatoires (équation (1.11)) conduit à un problème de fermeture. $\overline{f'g'}$ est un terme de corrélation turbulente qui nécessite donc une modélisation physique qui découlera directement de la signification phénoménologique à laquelle $\overline{f'g'}$ se rattache.

En appliquant la décomposition de Reynolds (équation (1.10)) aux variables turbu-
lentes qui interviennent dans le système d'équations défini (équations (1.1) à (1.3)) et
en moyennant<fnref n="4" />, les équations moyennes du mouvement turbulent dans la couche limite

<footnote n="4">4. Procédure et calculs en annexe A.3</footnote>
atmosphérique sont données par :

- L'équation moyenne de conservation de masse :

$$ \frac{\partial \overline{u_j}}{\partial x_j} = 0 \qquad (1.12) $$

- L'équation moyenne de conservation de quantité de mouvement ou équation de Reynolds :

$$ \rho_r \frac{\partial \overline{u_i}}{\partial t} + \rho_r \overline{u}_j \frac{\partial \overline{u_i}}{\partial x_j} = \rho_r \beta \overline{T}_l \delta_{i3} - \frac{\partial \overline{P}}{\partial x_i} + \frac{\partial}{\partial x_j} \left[ \mu \left( \frac{\partial \overline{u_i}}{\partial x_j} + \frac{\partial \overline{u_j}}{\partial x_i} \right) - \rho_r \overline{u'_i u'_j} \right] - 2 \rho_r \epsilon_{ijk} \Omega_j \overline{u_k} \quad (1.13) $$

$\beta$ définit le coefficient de flottabilité $\beta = g/T_r$

- L'équation moyenne de l'énergie :

$$ \rho_r C_P \left[ \frac{\partial \overline{T_l}}{\partial t} + \overline{u}_j \frac{\partial \overline{T_l}}{\partial x_j} \right] = \frac{\partial}{\partial x_j} \left[ \lambda \frac{\partial \overline{T_l}}{\partial x_j} - \rho_r C_P \overline{T' u'_j} \right] \quad (1.14) $$

L'écoulement moyen est aussi à divergence nulle. Le traitement statistique de l'équation instantanée de conservation de quantité de mouvement entraîne l'apparition de $\rho_r \overline{u_i' u_j'}$: ces termes traduisent l'interaction entre le mouvement moyen et le mouvement turbulent. Ils reflètent donc l'effet de la turbulence et restent à priori inconnus bien qu'ils contiennent des termes prépondérants. De la même façon, le traitement statistique de l'équation de l'énergie conduit à un terme inconnu $\rho_r C_P \overline{T' u'_j}$. Ces corrélations doubles rendent donc le système ouvert. Dans la suite de ce travail, les notations $\rho_r, T_l, P_l$ seront volontairement remplacées par $\rho, T, P$ par souci de simplification.

## 1.3 Les modèles de fermeture

Il existe de nombreuses méthodes de fermeture, notamment celles qui consistent à écrire et à résoudre les équations d'évolution des corrélations doubles $\rho \overline{u_i' u_j'}$, $\rho C_P \overline{T' u'_j}$. Cependant le traitement statistique de Reynolds va entraîner l'apparition de corrélations triples, qui requièrent une modélisation à leur tour. En effet, chaque corrélation d'ordre $m$ fait apparaître dans son équation d'évolution de nouvelles corrélations d'ordre $m + 1$. Ces méthodes de fermeture dites à l'ordre $m$ peuvent être retrouvées dans [Schiestel 93, Laurence 97].
L'autre démarche très souvent utilisée, s'appuie sur une modélisation phénoménologique des corrélations doubles en faisant une analogie entre les transferts de types diffusifs par agitation turbulente et par agitation moléculaire. Ces modèles appartiennent à la classe des modèles à concept de viscosité turbulente, que nous avons choisi d'utiliser et que nous allons présenter.

### 1.3.1 Modélisation phénoménologique des corrélations

Soit le flux turbulent moyen d'une variable instantanée $\phi$ transportable, le transfert de la quantité fluctuante $\phi'$ par la vitesse fluctuante $u'_j$ (figure 1.4) est considéré comme proportionnel au gradient de la valeur moyenne $\bar{\phi}$ :

<figure><img src="image_10.png" /><figcaption>Figure 1.4 : Flux turbulent moyen de la quantité $\phi$</figcaption></figure>

$$ 
\begin{aligned}
\phi' &\approx l' \frac{\partial \bar{\phi}}{\partial x_j} \\
u'_j &\approx -u'
\end{aligned}
\qquad (1.15)
$$

où $u'$ et $l'$ sont successivement l'échelle de vitesse et l'échelle de longueur représentatives du mouvement turbulent. Le couple $u'l'$ de dimension $[L^2T^{-1}]$ définit la diffusivité turbulente. En utilisant les relations (1.15), le flux turbulent moyen $\overline{\phi' u'_j}$ s'écrit :

$$ \overline{\phi' u'_j} = -(u'l') \frac{\partial \bar{\phi}}{\partial x_j} \qquad (1.16) $$

Ce transfert s'exprime suivant une relation linéaire entre la diffusivité turbulente $k_{\phi t}$ associée à la variable $\phi$ et le gradient de $\bar{\phi}$. Il s'identifie donc à un transfert par diffusion
turbulente :

$$ \overline{\phi'u'_j} = -k_{\phi t} \frac{\partial \overline{\phi}}{\partial x_j} \quad \text{avec} \quad k_{\phi t} = u'l' \qquad (1.17) $$

Contrairement à la diffusivité moléculaire qui est une grandeur intrinsèque du fluide, $k_{\phi t}$ est une propriété de l'écoulement. Il en résulte que la diffusion turbulente est in-comparablement plus efficace que la diffusion moléculaire (sauf au contact d'une paroi). Considérons pour cela, la quantité de mouvement $\phi = \rho u'_i$ et notons par $\nu_t$ la diffusivité turbulente associée. Le rapport des viscosités cinématiques s'approxime à :

$$ \frac{\nu_t}{\nu} \approx \frac{u'l'}{\nu} \qquad (1.18) $$

$u'$ est d'un ordre de grandeur en-dessous des écarts moyens caractéristiques $u' \approx \Delta \overline{u}/10$ et $l' \approx L/3$ [Coantic 76] $L$ étant une longueur de référence de l'écoulement, l'équation (1.18) est telle que :

$$ \frac{\nu_t}{\nu} \approx \frac{u'l'}{\nu} \approx \frac{\frac{\Delta \overline{u} L}{10}}{\frac{3}{\nu}} \approx \frac{R_e}{30} \qquad (1.19) $$

Le nombre de Reynolds est grand puisque le régime de l'écoulement est turbulent, et en généralisant, la diffusivité moléculaire associée à une variable quelconque est négligeable devant la diffusivité turbulente :

$$ k_{\phi} \ll k_{\phi t} \approx u'l' \qquad (1.20) $$

### Conséquence sur les corrélations $\rho u'_i u'_j$

Les corrélations $\rho u'_i u'_j$ correspondent donc à un flux turbulent moyen de quantité de mouvement dans la direction $j$. Ils sont assimilés à une contrainte turbulente (frottement turbulent) et forment le tenseur de Reynolds dont la première expression a été introduite par Boussinesq (1877) :

$$ \rho \overline{u'_i u'_j} = -\mu_t \left[ \frac{\partial \overline{u_i}}{\partial x_j} + \frac{\partial \overline{u_j}}{\partial x_i} \right] \quad \text{avec} \quad \mu_t = \rho \nu_t \qquad (1.21) $$

$\mu_t$ (ou $\nu_t$ qui est simplement relié à $\mu_t$ par la masse volumique) est appelé viscosité turbulente.

### Conséquence sur les corrélations $\rho C_P T' u'_j$

$\rho C_P T' u'_j$ est le flux de chaleur turbulent dans la direction $j$. Il admet une modélisation phénoménologique et s'écrit :

$$ \rho C_P \overline{T' u'_j} = -\lambda_t \frac{\partial \overline{T}}{\partial x_j} \quad \text{avec} \quad \lambda_t = \rho C_P k_t \qquad (1.22) $$
$k_t$ est la diffusivité thermique turbulente. En s'appuyant sur l'analogie de Reynolds, $\rho C_P \overline{T'u'_j}$ est donné par :

$$ \rho C_P \overline{T'u'_j} = - \frac{C_P \mu_t \partial \overline{T}}{P_{rt} \partial x_j} \quad \text{où} \quad P_{rt} = \frac{\nu_t}{k_t} \quad \text{et} \quad \lambda_t = \frac{C_P \mu_t}{P_{rt}} \quad (1.23) $$

$P_{rt}$ est le nombre de Prandtl turbulent défini comme le rapport du coefficient de transfert turbulent de quantité de mouvement $\nu_t$ à la diffusivité thermique turbulente $k_t$.

Contrairement au nombre de Prandtl $P_r$ qui dépend des propriétés intrinsèques du fluide, dans le cas des écoulements turbulents, c'est le processus dynamique qui fixe l'efficacité de la diffusivité thermique à travers $P_{rt}$. En admettant que les processus turbulents sont similaires pour la dynamique et pour la thermique quel que soit le fluide, la valeur du nombre de Prandtl turbulent est d'un point de vue théorique proche de l'unité. Toutefois, cette hypothèse purement intuitive fait que la valeur de $P_{rt}$ est empreint d'un fort empirisme et explique les différences notoires rencontrées dans la littérature.

### 1.3.2 Système d'équations modélisées

En explicitant les flux turbulents, le système d'équations modélisées de la couche limite atmosphérique est composé de :

- L'équation moyenne de conservation de masse :

$$ \frac{\partial \overline{u_j}}{\partial x_j} = 0 \quad (1.24) $$

- L'équation moyenne de conservation de quantité de mouvement :

$$ \begin{aligned} \rho \frac{\partial \overline{u_i}}{\partial t} + \rho \overline{u_j} \frac{\partial \overline{u_i}}{\partial x_j} &= \rho \beta \overline{T} \delta_{i3} - \frac{\partial \overline{P}}{\partial x_i} \\ &\quad + \frac{\partial}{\partial x_j} \left[ (\mu + \mu_t) \left( \frac{\partial \overline{u_i}}{\partial x_j} + \frac{\partial \overline{u_j}}{\partial x_i} \right) \right] - 2\rho \epsilon_{ijk} \Omega_j \overline{u_k} \end{aligned} \quad (1.25) $$

- L'équation moyenne de l'énergie :

$$ \rho C_P \left[ \frac{\partial \overline{T}}{\partial t} + \overline{u}_j \frac{\partial \overline{T}}{\partial x_j} \right] = \frac{\partial}{\partial x_j} \left[ \left( \lambda + \frac{C_P \mu_t}{P_{rt}} \right) \frac{\partial \overline{T}}{\partial x_j} \right] \quad (1.26) $$

Le problème de fermeture réside maintenant dans la détermination de la viscosité turbulente et de manière générale de la diffusivité turbulente.
La première formulation de $\nu_t$ a été introduite par Prandtl (1925). C'est un schéma de fermeture à zéro équation de transport, très connu sous le nom de longueur de mélange. Ce modèle calcule le tenseur de Reynolds en introduisant une seule échelle de longueur de la turbulence $l_m = l'$ déterminée par des considérations géométriques. Près des parois, Prandtl suggère que la longueur de mélange varie de façon linéaire avec la distance à la paroi, soit $l_m = \kappa z$. La longueur de mélange joue dans le mécanisme de transport un rôle analogue à celui du libre parcours moyen pour la théorie cinétique des gaz. L'échelle de vitesse de la turbulence $u'$ est reliée à $l_m$ par :

$$ u' = l_m \left| \frac{\partial \overline{u_i}}{\partial x_j} + \frac{\partial \overline{u_j}}{\partial x_i} \right| \qquad (1.27) $$

La diffusivité turbulente de la quantité de mouvement s'écrit alors :

$$ \nu_t = u'l' = l_m^2 \left| \frac{\partial \overline{u_i}}{\partial x_j} + \frac{\partial \overline{u_j}}{\partial x_i} \right| \qquad (1.28) $$

Bien que ces modèles à zéro-équation de transport donnent de bons résultats pour des écoulement cisaillés simples, ils souffrent d'un manque d'universalité puisque la longueur $l_m$ est choisie empiriquement et suppose au préalable une connaissance de l'écoulement à calculer. De plus, la relation (1.27) impose à $u'$ une dépendance trop importante vis-à-vis de l'écoulement moyen.

### 1.3.3 Modèle de fermeture à une équation de transport

Afin d'enlever l'empirisme qui réside dans la formulation de la longueur de mélange, l'échelle de vitesse caractéristique de la turbulence, qui est reliée directement à l'intensité de la turbulence, va s'exprimer en fonction de l'énergie cinétique moyenne du mouvement turbulent $k$ :

$$ u' = \sqrt{k} \quad \text{avec} \quad k = \frac{\overline{u_i' u_i'}}{2} \qquad (1.29) $$

Ce modèle est dit à une équation de transport, car il est nécessaire de disposer de l'équation d'évolution de l'énergie cinétique moyenne du mouvement turbulent. Le tenseur de Reynolds défini par Boussinesq (équation (1.21)) est modifié et s'écrit à l'aide du terme isotrope d'énergie cinétique afin d'assurer l'égalité dans le cas de la contraction d'indice –car pour $i = j$, le tenseur de Reynolds est nul pour un écoulement incompressible et impliquerait que $k = 0$, ce qui est une aberration– et attribue à l'énergie cinétique le rôle d'une pression turbulente due aux mouvements d'agitation.

$$ \rho \overline{u_i' u_j'} = -\mu_t \left[ \frac{\partial \overline{u_i}}{\partial x_j} + \frac{\partial \overline{u_j}}{\partial x_i} \right] + \frac{2}{3} \rho \delta_{ij} k \qquad (1.30) $$
L'équation moyenne de conservation de quantité de mouvement s'écrit maintenant :

$$
\rho \frac{\partial \overline{u_i}}{\partial t} + \rho \overline{u}_j \frac{\partial \overline{u_i}}{\partial x_j} = \rho \beta \overline{T} \delta_{i3} - \frac{\partial}{\partial x_i} \left[ \overline{P} + \frac{2}{3} \rho k \right] \\
+ \frac{\partial}{\partial x_j} \left[ (\mu + \mu_t) \left( \frac{\partial \overline{u_i}}{\partial x_j} + \frac{\partial \overline{u_j}}{\partial x_i} \right) \right] - 2\rho \epsilon_{ijk} \Omega_j \overline{u}_k \quad (1.31)
$$

La procédure pour établir l'équation exacte de k se trouve en annexe A.4. Cette équation est constituée par de nombreux termes qui ne sont pas accessibles directement par le calcul. Aussi, modéliser cette équation revient à classer les différents termes dans des « groupes catalogués ». L'équation de l'énergie cinétique moyenne du mouvement turbulent modélisée est donnée par :

$$
\rho \frac{\partial k}{\partial t} + \rho \bar{u}_j \frac{\partial k}{\partial x_j} = \underbrace{\frac{\partial}{\partial x_j} \left[ \left( \mu + \frac{\mu_t}{\sigma_k} \right) \frac{\partial k}{\partial x_j} \right]}_{\text{I}} + \underbrace{\rho P_k + \rho P_h}_{\text{II}} - \underbrace{\rho \varepsilon}_{\text{III}} \quad (1.32)
$$

-I- La variation de l'énergie cinétique due au transport par diffusion sous l'action de la turbulence s'écrit suivant la relation flux-gradient que nous venons de présenter. La prise en compte du transport par pression fluctuante se fait d'une manière artificielle en introduisant un nombre de Prandtl-Schmidt $\sigma_k$;

-II- $P_k$ est le taux de production cinétique (ou par cisaillement) du mouvement turbulent. Il traduit une transformation de l'énergie cinétique du mouvement moyen au profit du mouvement turbulent;

$$
P_k = - \overline{u'_i u'_j} \frac{\partial \overline{u_i}}{\partial x_j} \qquad (1.33)
$$

$P_h$ est la production gravitationnelle qui est directement reliée à la stratification thermique de la couche limite, dont nous parlerons à la section suivante;

$$
P_h = \beta \overline{T' u'_i} \delta_{i3} \quad (1.34)
$$

-III- $\epsilon$ est le taux de dissipation de l'énergie cinétique du mouvement turbulent (processus irréversible) son calcul est donc fondamental. Plusieurs démarches ont été élaborées pour le calcul de $\epsilon$ dont la plus célèbre est sans aucun doute la théorie de Kolmogorov (1941) [Ha-Minh 91]. La turbulence est créée et alimentée en énergie à grande échelle à travers son taux de production, énergie qui va se dissiper au même taux, en chaleur par effet de viscosité. Ce flux d'énergie est égal, à une constante
près, au rapport de l'énergie produite par les gros tourbillons $u'^2$ par le temps ca-
ractéristique des gros tourbillons $\tau = l'/u'$ ($\tau$ étant le temps de recouvrement des
gros tourbillons). Ce flux d'énergie n'est autre que $\varepsilon$, soit :

$$
\varepsilon = C_{l\varepsilon} \frac{u'^2}{\tau} = C_{l\varepsilon} \frac{u'^3}{l'} \quad (1.35)
$$

En considérant la relation (1.29), le taux de dissipation s'écrit :

$$
\epsilon = C_{l\epsilon} \frac{k^{\frac{3}{2}}}{l_m} \qquad (1.36)
$$

Finalement, l'équation d'évolution de l'énergie cinétique moyenne du mouvement turbu-
lent est donnée par :

$$
\rho \frac{\partial k}{\partial t} + \rho \bar{u}_j \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 P_h - \rho C_{l\epsilon} \frac{k^{\frac{3}{2}}}{l_m} \quad (1.37)
$$

<table><caption>Tableau 1.1 : Longueur de mélange et constantes du modèle k − l dans la couche de surface
( z < 200 m), en régime de convection forcée [Pielke 84, Chassaing 00]</caption>
  <tr>
    <td>l<sub>m</sub></td>
    <td>&sigma;<sub>k</sub></td>
    <td>C<sub>l&epsilon;</sub></td>
    <td>C<sub>l&mu;</sub></td>
    <td>P<sub>rt</sub></td>
  </tr>
  <tr>
    <td>&kappa;z</td>
    <td>1.66</td>
    <td>0.16</td>
    <td>0.54</td>
    <td>0.95</td>
  </tr>
</table>

La viscosité turbulente est donc calculée par :

$$
\mu_t = \rho C_{l\mu} l_m \sqrt{k} \quad (1.38)
$$

Le système d'équations (1.24), (1.26) et (1.31) complété des équations (1.37) et (1.38)
est fermé. Les constantes $C_{l\epsilon}$, $C_{l\mu}$, le nombre de Prandtl turbulent $P_{rt}$ ainsi que la lon-
gueur de mélange sont déterminés et évalués empiriquement à travers des expériences.
Les valeurs couramment utilisées sont regroupées dans le tableau 1.1.

1.3.4 Modèle de fermeture à deux équations de transport

Si le modèle *k* – *l* présente l'intérêt d'être simple à mettre en oeuvre, sa principale
difficulté n'en demeure pas moins dans la détermination d'une loi pour l'échelle de lon-
gueur *l*<sub>*m*</sub>: ce modèle souffre du même manque d'universalité que le modèle à zéro-équation
de transport. Pour contourner ce problème, l'idée naturelle qui est apparue est de procéder
de façon identique au modèle $k-l$, à savoir décrire l'évolution de la longueur caractéris-
tique $l_m$ à l'aide d'une équation de transport. Les nombreuses tentatives ont vite avorté car
il semble impossible d'établir de façon exacte, comme cela a été fait pour $k$, une équation
d'évolution de l'échelle de longueur à partir des équations de Navier-Stokes [Ha-Minh 91].

C'est en s'appuyant sur ces constats, qu'une autre approche quasi-identique a émergé,
et qui consiste non pas à rechercher l'évolution de la longueur caractéristique, mais celle
du couple $k^\alpha l_m^\gamma$. En particulier pour $\alpha = 3/2$ et $\gamma = -1$, le couple $k^{3/2}l_m^{-1}$ n'est autre que
le taux de dissipation –en fait, c'est la pseudo-dissipation$^5$– dont l'équation exacte peut
s'établir. Ce modèle de fermeture à deux-équations de transport est le très populaire mo-
dèle de fermeture $k-\epsilon$ [Launder 74]. La pseudo-dissipation qui s'identifie à la dissipation
lorsque le nombre de Reynolds de la turbulence est élevé, s'écrit :

$$
\epsilon = \nu \frac{\partial u_i' \partial u_i'}{\partial x_j \partial x_j} \qquad (1.39)
$$

L'équation de transport de la pseudo-dissipation est explicitée en annexe A.4. Elle se
compose de nombreux termes, qui ne sont malheureusement pas toujours très bien évalués,
aussi bien au niveau de l'interprétation physique qu'au niveau de leur valeur en amplitude
et, n'apporte que des informations globales sur les mécanismes de transfert d'énergie
entre les différentes échelles de la turbulence. Comme il a été fait pour l'équation de k,
cette équation va se diviser en plusieurs termes caractéristiques, avec toutefois des degrés
d'incertitudes amplifiés, ce qui explique les différences notables que l'on rencontre dans la
littérature. De façon générale, l'équation d'évolution modélisée du taux de dissipation a
la forme suivante :

$$
\rho \frac{\partial \varepsilon}{\partial t} + \rho \bar{u}_j \frac{\partial \varepsilon}{\partial x_j} = \underbrace{\frac{\partial}{\partial x_j} \left[ \left( \mu + \frac{\mu_t}{\sigma_\varepsilon} \right) \frac{\partial \varepsilon}{\partial x_j} \right]}_{\text{I}} + \underbrace{C_{\varepsilon 1} \rho \frac{\varepsilon}{k} (P_k + P_h)}_{\text{II}} - \underbrace{C_{\varepsilon 2} \rho \frac{\varepsilon^2}{k}}_{\text{III}} \quad (1.40)
$$

-I- Le transport par diffusion turbulente et visqueuse, ainsi que le transport par pres-
sion fluctuante est englobé dans ce terme. Le nombre de Prandtl-Schmidt turbulent
associé à ε est noté σε ;

-II- Dans des situations de turbulence homogène, Lumley (1972) a montré que les inter-
actions avec le mouvement moyen sont négligeables vis-à-vis des autres contributions
[Ha-Minh 91]. En effet, les termes d'interaction entre les composantes fluctuantes

<footnote n="5">5. Le taux de dissipation est défini par l'équation (A.41) en annexe A.4</footnote>
$P_k$ et $P_h$ représentent la production de la dissipation par étirement des vortex et sont prépondérants pour des nombres de Reynolds élevés;

-III- Terme de dissipation de l'équation de $\epsilon$ par action de la viscosité.

La modélisation des termes de production (II) et de dissipation (III) s'obtient à l'aide du concept d'équilibre spectral, qui stipule que les temps caractéristiques de la production et de la dissipation de $\epsilon$ sont proportionnels à ceux de $k$. En faisant une analyse dimensionnelle :

$$[T] \equiv \frac{[k]}{[P_k + P_h]} \quad \text{temps nécessaire pour produire de l'énergie}$$

$$[T] \equiv \frac{[k]}{[\epsilon]} \quad \text{temps nécessaire pour détruire l'énergie}$$

Soit, si nous notons par $P_\epsilon$ et $D_\epsilon$ les termes de production et de dissipation de $\epsilon$, il en résulte :

$$\begin{align*} \frac{k}{P_k + P_h} &= C_{\epsilon 1} \frac{\epsilon}{P_\epsilon} \Rightarrow P_\epsilon = C_{\epsilon 1} \frac{\epsilon}{k} (P_k + P_h) \\ \frac{k}{\epsilon} &= C_{\epsilon 2} \frac{\epsilon}{D_\epsilon} \Rightarrow D_\epsilon = C_{\epsilon 2} \frac{\epsilon^2}{k} \end{align*}$$

Finalement, la longueur de mélange disparaît de la formulation de la viscosité turbulente, qui devient :

$$\mu_t = \rho C_\mu \frac{k^2}{\epsilon} \qquad (1.41)$$

Ce système composé des équations (1.32), (1.40) et (1.41) forme le modèle $k - \epsilon$ dit standard avec les valeurs des constantes (tableau 1.2) évaluées par Launder & Spalding (1974) et Hanjalic & Launder (1976) [Launder 74]. La constante $C_\mu$ a été déterminée expérimentalement, en considérant une couche de cisaillement où l'équilibre local (production=dissipation) est réalisé. Les autres constantes sont évaluées à partir de résultats d'expériences de décroissance d'une turbulence homogène à l'aval d'une grille [Laurence 97].

Le modèle $k - \epsilon$ est actuellement des plus utilisés, notamment dans les codes de calcul industriel, car il présente l'intérêt d'être simple à mettre en oeuvre et il est connu pour sa robustesse. Cependant, il présente en même temps de nombreuses faiblesses. Les constantes sont ajustées de façon empirique afin de caler au mieux à la physique des écoulements étudiés et le modèle est trop diffusif.

C'est pour remédier à ces principaux défauts que d'autres modèles du types $k-\epsilon$ ont été développés. Ces modèles sont passés en revue dans [Glockner 00] et [Delboulbé 97], et nous
<table><caption>Tableau 1.2 : *Constantes du modèle k − ε standard* [Launder 74]</caption><thead><tr><th>σ<sub>k</sub></th><th>σ<sub>ε</sub></th><th>P<sub>rt</sub></th><th>C<sub>μ</sub><sup>a</sup></th><th>C<sub>ε1</sub></th><th>C<sub>ε2</sub></th></tr></thead><tbody><tr><td>1</td><td>1.3</td><td>0.95</td><td>0.09</td><td>1.44</td><td>1.92</td></tr></tbody></table>

<table><caption>Tableau 1.3 : *Constantes du modèle k − ε étendu* [Chen 87]</caption><thead><tr><th>σ<sub>k</sub></th><th>σ<sub>ε</sub></th><th>P<sub>rt</sub></th><th>C<sub>μ</sub></th><th>C<sub>ε1</sub></th><th>C<sub>ε2</sub></th><th>C<sub>ε3</sub></th></tr></thead><tbody><tr><td>0.75</td><td>1.15</td><td>0.95</td><td>0.09</td><td>1.15</td><td>1.9</td><td>0.25</td></tr></tbody></table>

ne retiendrons que le modèle $k-ε$ étendu, développé par Chen & Kim (1987) qui améliore le caractère trop diffusif du modèle standard par le biais d'un terme supplémentaire (1.42-IV) dans l'équation du taux de dissipation (1.40):

$$ (1.40) + \underbrace{\rho \frac{C_{\varepsilon 3}}{k} (P_k + P_h)^2}_{\text{IV}} \qquad (1.42) $$

Ce terme est calculé en introduisant une autre échelle de temps de production de l'énergie cinétique moyenne du mouvement turbulent et permet au transfert d'énergie de répondre plus efficacement aux distorsions de l'écoulement moyen [Chen 87]. Un nouveau jeu de constantes est utilisé (tableau 1.3).

Qu'il soit sous forme standard [Costes 96, Delboulbé 97], ou encore sous forme adaptée [Detering 85, Duynkerke 88], le modèle de turbulence $k-ε$ est souvent utilisé pour modéliser le vent et la turbulence atmosphérique.

## 1.4 Modélisation des échanges turbulents près du sol

Nous avons défini un modèle dynamique décrit par les équations de Reynolds simplifiées dans le cadre des hypothèses de Boussinesq. Le modèle de turbulence choisi, que ce soit le modèle $k-ε$ ou le modèle $k-l$, se base sur le concept de viscosité turbulente. Ces modèles sont bien adaptés à la modélisation de la couche limite atmosphérique développée dont le nombre de Reynolds est élevé. Près du sol, plus précisément dans la sous-couche

<footnote><sup>a</sup> Les constantes du modèle $k-l$ sont donc définies par $C_{lε} = C_{μ}^{3/4}$ et $C_{lμ} = C_{μ}^{1/4}$</footnote>
rugueuse, ont lieu les échanges par convection de la quantité de mouvement et de la cha
leur entre le sol et l'atmosphère. La résolution de l'écoulement à cette échelle est donc
nécessaire, mais exige l'utilisation de maillages fins, car les gradients sont intenses. De
plus, le vent présente de fortes irrégularités, la turbulence y est très intermittente et ins-
tationnaire : les modèles statistiques donnent des résultats peu satisfaisants pour ce types
d'écoulements.

Pour contourner ce problème et s'affranchir d'une description fine de l'écoulement à cette échelle, nous avons choisi d'utiliser des lois de paroi empiriques [Launder 74] qui permettent de modéliser les flux turbulents et surtout, de reproduire leurs effets sur la couche de surface. L'établissement des lois de parois s'appuient sur la propriété de couche à flux constant de la couche de surface à laquelle nous allons restreindre la discussion aux paragraphes suivants.

1.4.1 La couche de surface sur un sol homogène: hypothèses et
équations

Dans la couche de surface idéale, pour un écoulement développé suivant la direction x,
le frottement total et le flux de chaleur sensible est quasi-constant (la variation suivant
z reste inférieure à 20%) au voisinage immédiat de la paroi. Les influences du gradient
moyen de pression et la force de Coriolis sont négligeables. Dans ces conditions, l'équation
moyenne de conservation de quantité de mouvement et de conservation de l'énergie se
réduisent à :

$$
\mu \frac{\partial \bar{u}}{\partial z} - \rho \overline{u'w'} = \tau_0 \quad (1.43)
$$

$$
-\lambda \frac{\partial \overline{T}}{\partial z} + \rho C_P \overline{T'w'} = Q_0 \quad (1.44)
$$

$\tau_0$ [$kg m^{-1} s^{-2}$] représente le frottement (ou contrainte de cisaillement) total à la surface et $Q_0$ [$W m^{-2}$] est la densité de flux de chaleur sensible.

Au sein de la couche de surface, les transferts sont déterminés par l’interaction avec la
surface et les flux turbulents sont constants et égaux à leur valeur à la surface. Ces consi-
dérations entraînent la définition de deux échelles caractéristiques dont les valeurs sont
constantes dans la couche de surface : l’échelle de vitesse $u_*$, appelée vitesse de frottement
(ou de cisaillement) et l’échelle de température $T_*$ appelée température de frottement.
Ces échelles sont définies par :

$$u_* = \left( \frac{\tau_0}{\rho} \right)^{\frac{1}{2}} \qquad (1.45)$$

$$u_* T_* = -\frac{Q_0}{\rho C_P} \qquad (1.46)$$

En utilisant ces notations, les équations dans la couche à flux constant obéissent à :

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

$$-\lambda \frac{\partial \bar{T}}{\partial z} + \rho C_P \overline{T'w'} = -\rho C_P u_* T_* \qquad (1.48)$$

Soit, en exprimant les flux turbulents de quantité de mouvement et de chaleur, il vient :

$$(\mu + \mu_t) \frac{\partial \bar{u}}{\partial z} = \rho u_*^2 \qquad (1.49)$$

$$-(\lambda + \lambda_t) \frac{\partial \bar{T}}{\partial z} = -\rho C_P u_* T_* \qquad (1.50)$$

De ces équations, il est possible d'expliciter des lois de variations de la vitesse et de la température au voisinage du sol pour des surfaces lisses et rugueuses.

## 1.4.2 Interaction dynamique de l'écoulement avec une paroi

### Paramètres pertinents pour une paroi dynamiquement lisse $h_s < 5z_*$

L'analyse dimensionnelle permet de suggérer des relations (sans toutefois les prouver) et les seuls paramètres pertinents qui décrivent la variation de $\bar{u}(z)$ sont $u_*$, $z$, et $z_*$ appelé longueur de frottement :

$$z_* = \frac{\nu}{u_*} \qquad (1.51)$$

Dans ce cas précis, il existe une relation telle que :

$$\frac{\bar{u}(z)}{u_*} = f_u \left( \frac{z}{z_*} \right) \qquad (1.52)$$

À la paroi, la condition d'adhérence implique que le flux turbulent est nul, le flux moléculaire est donc dominant dans l'équation (1.47) : cette zone mince près de la paroi est nommée sous-couche visqueuse et définie pour $z < 5z_*$. Le frottement est uniforme et la vitesse varie linéairement suivant $z$ :

$$\frac{\bar{u}(z)}{u_*} = \frac{z}{z_*} \quad \text{pour } z < 5z_* \qquad (1.53)$$
En revanche, lorsque $z$ est suffisamment éloigné de la paroi, le flux turbulent n’est plus négligeable, et devient supérieur au flux moléculaire. Cette zone s’appelle la sous-couche inertielle. Dans ces conditions, la relation (1.47) devient :

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

L'équation (1.54) fait apparaître que la variation de $\overline{u}(z)$ ne dépend plus de la viscosité cinématique $\nu$. En faisant l'hypothèse que les échelles de vitesse et de longueur définissant la viscosité turbulente sont $u_*$ et $z$, $\nu_t$ est déterminée par l'expression :

$$ \nu_t = \kappa u_* z \qquad (1.55) $$

Et le gradient de la vitesse $\overline{u}(z)$ s'écrit :

$$ \frac{\partial \overline{u}}{\partial z} = \frac{u_*}{\kappa z} \qquad (1.56) $$

$\kappa$ désigne la constante numérique universelle de Karman, dont la valeur sera prise égale à $\kappa = 0.41$ [Ha-Minh 91]. L'intégration de (1.56) conduit à une loi logarithmique pour le profil de vitesse moyenne :

$$ \overline{u}(z) = \frac{u_*}{\kappa} \ln(z) + A \qquad (1.57) $$

$A$ est une constante d'intégration. En supposant que les deux zones visqueuse et turbulente se raccordent à l'interface $z = z_\nu$ (continuité de la vitesse), la constante $A$ s'écrit :

$$ A = \overline{u}(z_\nu) - \frac{u_*}{\kappa} \ln(z_\nu) \qquad (1.58) $$

Compte tenu de (1.57) et (1.58), il vient :

$$ \overline{u}(z) = \frac{u_*}{\kappa} \ln\left(\frac{z}{z_\nu}\right) + \overline{u}(z_\nu) \qquad (1.59) $$

Puisque $A$ est pris à la limite de validité des deux lois, $\frac{\overline{u}(z_\nu)}{u_*}$ vérifie (1.53) :

$$ \frac{\overline{u}(z)}{u_*} = \frac{1}{\kappa} \ln \left( \frac{z}{z_\nu} \right) + \frac{z_\nu}{z_*} \qquad (1.60) $$

Qui peut se mettre sous la forme de (1.52) :

$$ \frac{\overline{u}(z)}{u_*} = \frac{1}{\kappa} \ln \left( \frac{z}{z_*} \right) \underbrace{- \frac{1}{\kappa} \ln \left( \frac{z_\nu}{z_*} \right) + \frac{z_\nu}{z_*}}_{B^+ \left( \frac{z_\nu}{z_*} \right)} \qquad (1.61) $$
Il est à remarquer que la constante $B^+\left(\frac{z_\nu}{z_*}\right)$ (sans dimension) est fonction de l'épaisseur $z_\nu$. Expérimentalement, elle est évaluée en considérant que l'épaisseur $z_\nu$, qui correspond à la hauteur où les effets de viscosité moléculaire et de turbulence sont approximativement égaux, prend la valeur de $11.2z_*$ [Ha-Minh 91]. Soit :

$$B^+ (11.2) \approx 5.3 \qquad (1.62)$$

Reformulons (1.61),

$$\frac{\bar{u}(z)}{u_*} = \frac{1}{\kappa} \ln \left( \frac{z}{z_*} \exp \left[ \kappa B^+ \left( \frac{z_\nu}{z_*} \right) \right] \right) \qquad (1.63)$$

En posant :

$$E_0 = \exp \left[ \kappa B^+ \left( \frac{z_\nu}{z_*} \right) \right] \quad \text{Avec } B^+ (11.2) \approx 5.3 \quad E_0 \approx 8.8 \qquad (1.64)$$

$$z_0 = \frac{z_*}{E_0} \qquad (1.65)$$

En utilisant ces grandeurs, la variation de la vitesse moyenne dans la sous couche inertielle est donnée par :

$$\bar{u}(z) = \frac{u_*}{\kappa} \ln \left( \frac{z}{z_0} \right) \quad \text{pour} \quad 30z_* < z < 500z_* \qquad (1.66)$$

$z_0$ porte le nom de longueur de rugosité dynamique (1.65). Elle traduit de façon macroscopique l'effet de la rugosité de la surface sur le profil de vitesse. Dans le cas d'une paroi lisse, elle n'est fonction que du fluide, par l'intermédiaire de $\nu$ et de l'écoulement par l'intermédiaire de $u_*$. Le domaine de validité de la sous couche inertielle est environ $30z_* < z < 500z_*$. La sous couche visqueuse et la sous couche inertielle sont séparées par une zone tampon (buffer layer), pour $z$ variant entre $5z_* \le z \le 30z_*$, dans laquelle il y a un passage progressif des deux lois décrivant la vitesse moyenne.

## Paramètres pertinents pour une paroi dynamiquement rugueuse $h_s \ge 5z_*$

Une paroi est aérodynamiquement lisse tant que les rugosités $h_i$ restent dans la sous couche visqueuse, et donc ne perturbent pas l'écoulement. Si $h_s$ est un paramètre de longueur égale à la valeur moyenne des hauteurs des aspérités, cela s'exprime donc par $h_s < 5z_*$. En revanche si $h_s$ est supérieur à $5z_*$, l'écoulement est fortement perturbé par les irrégularités présentes à la surface. Il est donc légitime de rajouter $h_s$ comme paramètre pertinent pour décrire l'écoulement puisqu'il est influencé par la caractéristique de la surface. La relation (1.52) n'est plus valable et l'analyse dimensionnelle conduit à :

$$\frac{\bar{u}(z)}{u_*} = f_u \left( \frac{z}{z_*}, \frac{h_s}{z_*} \right) \qquad (1.67)$$
Connaissant $h_s$, il est d'usage de définir un nombre sans dimension, homogène à un nombre de Reynolds local caractéristique de la rugosité de surface (nombre de Reynolds de rugosité) :

$$H^+ = \frac{h_s}{z_*} = \frac{h_s u_*}{\nu} \qquad (1.68)$$

À une hauteur $z$ prise à une distance suffisamment éloignée de la paroi, le gradient de $\bar{u}(z)$ s'exprime toujours par (1.56), et en intégrant :

$$\frac{\bar{u}(z)}{u_*} = \frac{1}{\kappa} \ln \left( \frac{z}{z_*} \right) + C^+ (H^+) \qquad (1.69)$$

La constante d'intégration $C^+(H^+)$ dépend évidemment de $h_s$. Considérons deux cas extrêmes afin d'expliciter la loi de vitesse pour une paroi rugueuse :

- Lorsque $h_s$ tend vers 0, soit $H^+ < 5$, la paroi est dynamiquement lisse. La constante $C^+(H^+)$ va tendre vers $B^+=5.3$. Nous retrouvons ainsi la loi de paroi établie précédemment (équation (1.63));

- En revanche lorsque $h_s$ est très grand, c'est à dire $H^+ > 60$, l'écoulement est pleinement turbulent, et la zone tampon disparaît. La viscosité $\nu$ n'est plus un paramètre pertinent. La loi de variation de la vitesse moyenne $\bar{u}(z)$ s'écrit si nous explicitons $C^+(H^+)$ :

$$\frac{\bar{u}(z)}{u_*} = \frac{1}{\kappa} \ln \left( \frac{z}{z_*} \right) \underbrace{-\frac{1}{\kappa} \ln \left( \frac{h_s}{z_*} \right) + C'_+ \left( \frac{h_s}{z_*} \right)}_{C^+(H^+)} \qquad (1.70)$$

Les travaux expérimentaux de Nikuradse (1933) ont montré que la constante d'intégration $C'_+ \left(\frac{h_s}{z_*}\right)$ de la loi de vitesse écrite sous cette forme, n'est plus fonction de $H^+$. Sa valeur a été déterminée à :

$$C'_+ \left( \frac{h_s}{z_*} \right) = C^+ \approx 8.5 \qquad (1.71)$$

L'équation (1.70) peut se mettre sous la forme simplifiée :

$$\frac{\bar{u}(z)}{u_*} = \frac{1}{\kappa} \ln \left( \frac{z}{h_s} \right) + C^+ = \frac{1}{\kappa} \ln \left( \frac{z}{h_s} \exp [\kappa C^+] \right) \qquad (1.72)$$

Nous définissons à nouveau la longueur de rugosité $z_0$ qui dans le cas d'une paroi pleinement rugueuse, ne dépend plus que de la géométrie de la surface, caractérisée par $h_s$ :

$$z_0 = \frac{h_s}{30} \quad \text{avec } C^+ = 8.5 \qquad (1.73)$$
Afin de retrouver des notations similaires au cas lisse, et généraliser la loi de variation de la vitesse, nous écrivons :

$$ \frac{\bar{u}(z)}{u_*} = \frac{1}{\kappa} \ln \left( \frac{z}{z_*} \frac{z_*}{h_s} \exp[\kappa C_\infty^+] \right) \quad (1.74) $$

La constante sans dimension équivalente à (1.64) est évaluée dans le cas pleinement rugueux, par :

$$ E_{\infty} = \frac{z_{*}}{h_{s}} \exp[\kappa C_{\infty}^{+}] \quad E_{\infty} \approx \frac{30}{H^{+}} \quad (1.75) $$

Par conséquent, elle vaut :

$$ E_{\infty} = \frac{z_{*}}{z_{0}} \quad (1.76) $$

Finalement, le profil de vitesse pour une paroi rugueuse est donné par :

$$ \bar{u}(z) = \frac{u_{*}}{\kappa} \ln \left( \frac{z}{z_0} \right) \quad \text{pour} \quad z \gg z_0 \quad (1.77) $$

- Si la paroi est légèrement rugueuse, soit $5 \le H^+ \le 60$ la loi logarithmique (1.77) reste valable, cependant $z_0$ est fonction de $H^+$ et de $z_\nu$ à travers $E_0$. Les valeurs retenues pour $z_0$ sont [Guilloteau 99] :

$$ z_0 = \frac{z_*}{E(H^+)} \quad (1.78) $$

$$ E(H^+) = \frac{E_{\infty}}{\sqrt{\alpha + (1 - \alpha) \left(\frac{E_{\infty}}{E_0}\right)}} \quad (1.79) $$

Avec :

$$ E_0 = 8.8 \text{ et } E_\infty = \frac{30}{H^+} $$

$$ \alpha = 1 + x^2 (2x - 3) $$

$$ x = 0.02248 (100 - H^+) (H^+)^{-0.564} $$

### 1.4.3 Interaction thermique de l'écoulement avec une paroi

#### Paramètres pertinents pour une paroi dynamiquement lisse $h_s < 5z_*$

De façon similaire à la dynamique, les interactions de l'écoulement avec la paroi vont faire ressortir différentes zones, dans lesquelles les flux de chaleur moléculaire ainsi que le flux de chaleur turbulent vont plus ou moins gouverner l'équation (1.44). L'analyse
dimensionnelle introduit, compte tenu de l'équation de la chaleur, l'échelle pertinente $T_*$ pour décrire la variation de la température. Très près de la paroi, le flux de chaleur turbulent est négligeable et l'équation (1.50) se réduit à :

$$ -\lambda \frac{\partial \bar{T}}{\partial z} = -\rho C_P u_* T_* \qquad (1.80) $$

En considérant la température au sol $T_s$ et en intégrant l'expression ci-dessus, la température dans la sous-couche conductrice varie suivant une loi linéaire dont la pente est le nombre de Prandtl $P_r$. Ce nombre sans dimension correspond au rapport du coefficient de diffusion moléculaire visqueuse $\nu$ à la diffusivité thermique $k = \frac{\lambda}{\rho C_P}$ (k caractérise la capacité du fluide à conduire la chaleur) et caractérise la distribution de la température par rapport à la distribution des vitesses. Pour l'air, le nombre de Prandtl vaut 0.71.

$$ \frac{\bar{T}(z) - \bar{T}_s}{T_*} = P_r \frac{z}{z_*} \qquad (1.81) $$

En s'éloignant de la paroi, le flux de chaleur turbulent devient prépondérant dans l'équation (1.50), qui s'écrit :

$$ -\lambda_t \frac{\partial \bar{T}}{\partial z} = -\rho C_P u_* T_* \qquad (1.82) $$

En supposant que la diffusivité thermique turbulente s'exprime en fonction de l'échelle de vitesse $u_*$, de l'échelle de longueur $z$ et en introduisant le nombre de Prandtl turbulent $P_{rt}$, $\lambda_t$ s'écrit :

$$ \lambda_t = \rho C_P k_t \text{ avec } k_t = \frac{\nu_t}{P_{rt}} = \frac{\kappa u_* z}{P_{rt}} \qquad (1.83) $$

Le gradient de température est ainsi donné par la relation :

$$ \frac{\partial \bar{T}}{\partial z} = \frac{P_{rt} T_*}{\kappa z} \qquad (1.84) $$

Le gradient de température et le gradient de vitesse s'expriment par des relations de définitions équivalentes, à une constante de proportionnalité près $P_{rt}$. En intégrant la relation flux-gradient, la température est donnée par une loi logarithmique :

$$ \bar{T}(z) = \frac{P_{rt} T_*}{\kappa} \ln(z) + A \qquad (1.85) $$

Où A est une constante d'intégration. Considérons $z_k$ la hauteur prise à la limite de validité de la loi linéaire et de la loi logarithmique de la température. A peut s'écrire :

$$ A = \bar{T}(z_k) - \frac{P_{rt} T_*}{\kappa} \ln(z_k) \qquad (1.86) $$
En remplaçant $A$ dans l'équation (1.85), et en considérant que $\overline{T}(z_k)$ vérifie (1.84), il vient :

$$ \frac{\overline{T}(z) - \overline{T}_s}{T_*} = \frac{P_{rt}}{\kappa} \ln \left( \frac{z}{z_k} \right) + P_r \frac{z_k}{z_*} \quad (1.87) $$

Cette expression est équivalente à :

$$ \frac{\overline{T}(z) - \overline{T}_s}{T_*} = \frac{P_{rt}}{\kappa} \ln \left( \underbrace{\frac{z}{z_*} \frac{z_*}{z_k} \exp \left[ \frac{\kappa}{P_{rt}} \frac{\overline{T}(z_k) - \overline{T}_s}{T_*} \right]}_{B^+ \left( \frac{z_k}{z_*}, P_r, P_{rt} \right)} \right) \quad (1.88) $$

La constante $B^+ \left( \frac{z_k}{z_*}, P_r, P_{rt} \right)$ dépend du nombre de Prandtl turbulent $P_{rt}$, du nombre de Prandtl $P_r$ ainsi que de la hauteur $z_k$ de raccordement des lois. Reformulons l'équation (1.88) :

$$ \frac{\overline{T}(z) - \overline{T}_s}{T_*} = P_{rt} \left( \frac{1}{\kappa} \ln \left( E_0 \frac{z}{z_*} \right) + P_0 \right) \quad (1.89) $$

La constante $P_0 = \kappa^{-1} \ln \left( \frac{B^+}{E_0} \right)$ est due à la résistance thermique introduite par la sous-couche visqueuse au transfert de chaleur. Jayatilleke (1969) [Schiestel 93] estime cette constante dans le cas d'une paroi dynamiquement lisse à :

$$ P_0 = 9.24 \left[ \left( \frac{P_r}{P_{rt}} \right)^{0.75} - 1 \right] \left[ 1 + 0.28 \exp \left( -0.007 \frac{P_r}{P_{rt}} \right) \right] \quad (1.90) $$

En posant :

$$ z_{0h} = \frac{z_*}{E_0} \exp(-\kappa P_0) \quad (1.91) $$

$z_{0h}$ est la longueur de rugosité thermique. Elle caractérise de façon macroscopique l'effet de la rugosité de la surface sur le transfert de chaleur. Finalement, la variation de la température suit une loi logarithmique :

$$ \overline{T}(z) - \overline{T}_s = \frac{P_{rt} T_*}{\kappa} \ln \left( \frac{z}{z_{0h}} \right) \quad (1.92) $$

**Paramètres pertinents pour une paroi dynamiquement rugueuse $h_s \ge 5z_*$**

De la manière identique à la dynamique, les aspérités à la paroi vont avoir pour effet de modifier les échanges thermiques depuis la surface. En effet, la constante $P_0$ dépend de la caractéristique de la surface à travers $E_0$. Or nous avons vu que dans le cas d'une
paroi rugueuse, $E(H^+)$ est fonction du nombre de Reynolds de rugosité $H^+$. Afin d'inclure la résistance au transfert de chaleur induite par une surface rugueuse, Jayatilleke (1969) s'est inspiré des travaux menés pas Dipprey et Sabersky (1963) et donne l'expression de $P_0$ dans le cas rugueux :

$$P(H^+) = 3.15 P_r^{0.695} \left( \frac{1}{E(H^+)} - \frac{1}{E_0} \right)^{0.359} + 0.274 E(H^+)^{0.6} P_0 \quad (1.93)$$

La distinction entre le régime pleinement turbulent et légèrement turbulent ainsi que les expressions de $E(H^+)$ trouvés expérimentalement par Dipprey et Sabersky diffèrent des travaux de Nikuradse. Cependant, dans le cas des écoulements atmosphériques Guilloteau [Guilloteau 99] a montré que l'écart entre les deux formulations est négligeable lorsque $K^+$ est supérieur à $10^2$. Ainsi les constantes $E(H^+)$ que nous utiliserons sont issues des résultats de Nikuradse. Le profil de température moyenne est donné par :

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

La rugosité thermique $z_{0h}$ est, quant à elle, généralisée :

$$z_{0h} = \frac{z_*}{E(H^+)} \exp(-\kappa P(H^+)) = z_0 \exp(-\kappa P(H^+)) \quad (1.95)$$

### 1.4.4 Paramètres de rugosités des surfaces naturelles

Les mesures expérimentales de Nikuradse et Schiller, réalisées en laboratoire, ont permis de dégager une relation entre la hauteur moyenne des aspérités et $z_0$ dans le cas d'une paroi pleinement rugueuse :

$$z_0 = \frac{h_s}{30}$$

Dans les sites naturels, il a été observé que $z_0$ est plus grand que cette valeur proposée. En particulier, de nombreuses études ont été consacrées à la détermination de ce paramètre pour des couverts végétaux (tableau 1.4) et la relation appropriée pour déterminer $z_0$ se rapproche de [Kaimal 94] :

$$z_0 = \frac{h}{10} \quad (1.96)$$

$h$ étant la hauteur moyenne du couvert végétal. Szeicz et al. ont proposé une expression logarithmique [Guyot 97] :

$$\log_{10}(z_0) = 0.997 \log_{10}(h) - 0.883 \quad (1.97)$$
<table><caption>Tableau 1.4 : *Paramètre de rugosité pour différents types de surfaces naturelles*</caption><thead><tr><th>Type de surface</th><th>z<sub>0</sub> [m]</th></tr></thead><tbody><tr><td>Eau libre</td><td>0.002 ←→ 0.006</td></tr><tr><td>Sol nu</td><td>0.005 ←→ 0.02</td></tr><tr><td>Gazon ras (1 cm de haut)</td><td>0.001</td></tr><tr><td>Gazon dense (10 cm de haut)</td><td>0.023</td></tr><tr><td>Prairie (0.5 m de haut)</td><td>0.05 ←→ 0.07</td></tr><tr><td>Blé (1 m de haut)</td><td>0.1 ←→ 0.16</td></tr><tr><td>Maïs</td><td>0.1 ←→ 0.3</td></tr><tr><td>Forêt</td><td>0.3 ←→ 1.5</td></tr></tbody></table>

La zone où la présence de la canopée végétale influence directement la structure de la turbulence est appelée sous-couche rugueuse. Cette couche s'étend jusqu'à une hauteur $z$ égale à environ $2h$ [Raupach 96]. Au-dessus, le profil du vent est décrit par la loi logarithmique :

$$ 
\begin{aligned}
\bar{u}(z) &= \frac{u_*}{\kappa} \ln \left( \frac{z}{z_0} \right) \quad \text{pour } z \ge 2h \\
z_0 &= \frac{h}{10} \qquad (1.98)
\end{aligned}
$$

Afin de prolonger la zone de validité de la loi logarithmique à la hauteur $h$ du couvert, et donc représenter le vent à hauteur du couvert végétal, il est nécessaire de corriger le gradient de vitesse du vent [Kaimal 94]. Pour cela, l'origine du sol est surélevée d'une hauteur $d$ appelée hauteur de déplacement : $d$ est l'altitude d'une paroi fictive qui représente le « plancher » ressenti par l'écoulement [Costes 96]; elle est assimilée à la hauteur moyenne d'absorption de la quantité de mouvement par le couvert végétal [Kaimal 94]. La hauteur de déplacement est nécessairement comprise entre 0 et $h$. Plus précisément, sa valeur est voisine de :

$$ d = 0.7h \qquad (1.99) $$

Rosenberg (1974) précise l'expression de $d$ suivant la relation (1.100), en faisant des ob-
servations de vent au-dessus de différents types de parcelles agricoles [Pielke 84] :

$$log_{10}(d) = 0.979 log_{10}(h) - 0.154 \qquad (1.100)$$

Il est à rajouter que les paramètres $z_0$ et $d$ ne sont pas uniquement reliés à la hauteur du couvert, mais dépendent également de la structure et de la distribution géométrique des éléments du couvert végétal ainsi que de la densité des surfaces foliaires.

En utilisant la hauteur de déplacement $d$, les relations flux-gradients deviennent :

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

$$\frac{\partial \bar{T}}{\partial z} = \frac{P_{rt} T_*}{\kappa (z - d)} \qquad (1.102)$$

qui s'intègrent entre $z$ et $z_0$ pour la vitesse et $z$ et $z_{0h}$ pour la température. Finalement les profils de vitesse du vent et de température s'écrivent dans le cas d'une surface homogène en stratification thermique neutre :

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

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

Toutes ces expressions sont très utiles car elles donnent moyen de modéliser globalement les échanges turbulents entre l'écoulement et une surface homogène (mais qui peut être localement hétérogène) de paramètres moyens $z_0$ et $d$ sans résolution explicite à cette échelle. De surcroît, ces relations, couplées au modèle numérique, permettent de traiter un paysage hétérogène décrit comme une mosaïque de surfaces homogènes telle un « patchwork », dont les tailles et l'agencement spatial perturbent les basses couches de l'atmosphère en modifiant la structure de la couche de surface. En effet, chaque surface homogène provoque un forçage local sur l'écoulement dont les caractéristiques sont convectées en aval et sont diffusées dans l'atmosphère. Il était donc impératif de décrire ces effets de surface dans le modèle numérique puisque l'intensité turbulente, tributaire de la structure du paysage, est directement reliée à la capacité de l'écoulement du vent à transporter et à disperser des particules dans l'atmosphère.
### 1.5 Les effets de la stratification thermique de l'atmosphère

Lorsque l'atmosphère est en stratification neutre, c'est à dire que la température se comporte comme un scalaire passif, la turbulence se produit suivant un processus ciné-tique qui se rencontre notamment dans les parties de l'atmosphère où le cisaillement est important (à proximité du sol terrestre...). Ainsi, la turbulence est purement d'origine mécanique et dépend uniquement de la distribution verticale de la contrainte de cisaillement du sol. En revanche, lorsque l'atmosphère est en stratification thermique, la force d'Archimède va avoir un rôle non négligeable sur la turbulence.

Considérons une particule fluide, de densité différente de l'air qui l'entoure, que la turbulence a déplacée verticalement. Si la densité de l'air diminue avec l'altitude, la force d'Archimède va ramener cette particule fluide à son niveau initial (étant plus lourde). Ainsi les structures turbulentes vont effectuer en permanence un travail contre la pesanceur : l'atmosphère est dite en stratification stable, et la force d'Archimède a pour effet d'atténuer la turbulence. Dans le cas où la densité de l'air augmente avec l'altitude, la force d'Archimède va éloigner la particule fluide de sa position initiale (étant plus légère) et tend ainsi à accroître l'énergie de la turbulence. L'atmosphère est dite alors en stratification instable et la turbulence est renforcée. Dans les deux cas considérés, les processus turbulents sont d'origine dynamique et thermique [Coantic 76].

Donc la stabilité de l'atmosphère influence fortement les transferts turbulents. Il en ressort que les profils verticaux de température et de vitesse doivent être corrigés pour pouvoir représenter correctement les échanges turbulents dans la couche limite atmosphérique. Afin de prendre en compte la stratification thermique de l'atmosphère, nous avons choisi d'utiliser la théorie de similitude de Monin et Obukhov.

#### 1.5.1 Mesure de la stabilité et nature des transferts convectifs

En turbulence géophysique, l'influence de la stratification thermique se fait à travers le nombre de Richardson de flux $R_f$, nombre sans dimension introduit par Richardson en 1920 [Pielke 84]. $R_f$ est défini comme le rapport du terme de production d'origine thermique ou production gravitationnelle (1.34) sur le terme de production dynamique de l'énergie cinétique moyenne du mouvement turbulent ou production cinétique (1.33). Le signe '–' a été introduit par convention, et la première définition du nombre de Richardson
tenait compte de la vapeur d'eau dans les variations de densité [Dumas 74]. Dans la littérature, le nombre de Richardson de flux s'écrit simplement :

$$R_f = -\frac{P_h}{P_k} = -\beta \frac{\overline{T'w'}}{-\overline{u'w'} \frac{\partial \overline{u}}{\partial z} - \overline{v'w'} \frac{\partial \overline{v}}{\partial z}} \quad (1.105)$$

Le terme de production cinétique étant toujours positif, il en ressort que pour $\overline{T'w'} = -k_t \frac{\partial \overline{T}}{\partial z} > 0$, la stratification est instable puisque qu'elle correspond à un gradient thermique négatif (profil thermique suradiabatique représenté sur la figure 1.5). Le nombre de Richardson est dans ce cas négatif.

<figure><img src="image_11.png" /><figcaption>Figure 1.5 : Représentation schématique des profils thermiques diurne et nocturne [Guyot 97]</figcaption></figure>

Lorsque le gradient de température est positif (profil en condition d'inversion illustré sur la figure 1.5) $\overline{T'w'} < 0$, la stratification est stable et le nombre de Richardson est positif. La production gravitationnelle se comporte alors comme un terme puits vis-à-vis de la turbulence ($P_h < 0$), qui ne pourra se maintenir que si $P_d > -P_h$, ce qui est équivalent à $R_f \le R_{fc} < 1$. La dissipation étant également un terme puits, il apparaît que la valeur critique $R_{fc}$ est plus petite que l'unité; expérimentalement elle est évaluée à 0.25 [Kaimal 94, DeMoor 83]. Le régime neutre est atteint lorsque $\overline{T'w'} = 0$.

La stratification modifie donc l'origine de la turbulence, ainsi que sa structure et par conséquent la nature même des échanges convectifs entre le sol et l'atmosphère (figure 1.6). Pour $|R_f| \ll 1$, le régime de convection est de type convection forcée. La température se comporte comme un scalaire passif, et ces effets ne sont pas appréciables. Lorsque $R_f \to -\infty$, le régime est en convection libre, marqué par l'apparition de cellules convectives de plus en plus complexes, au détriment des mouvements horizontaux d'origine mécanique.
<figure><img src="image_12.png" /><figcaption>Figure 1.6 : *Effets du nombre de Reynolds et du nombre de Richardson sur la turbulence dans l'atmosphère* (Woods 1969)</figcaption></figure>

Entre ces limites, la convection est mixte, ce qui est le cas le plus fréquent dans la couche limite atmosphérique.

Un autre nombre de Richardson $R_i$, appelé nombre de Richardson de gradient (équation (1.106)), est couramment utilisé en l'approximant par le nombre de Richardson volumique (bulk Richardson number) $R_{ib}$ entre les hauteurs $z$ et $z_0$:

$$R_i = \beta \frac{\frac{\partial \bar{T}}{\partial z}}{\left(\frac{\partial \bar{u}}{\partial z}\right)^2 + \left(\frac{\partial \bar{v}}{\partial z}\right)^2} \quad (1.106)$$

$$R_{ib} = \beta \frac{(\bar{T} - \bar{T}_{z0h}) (z - z_0)^2}{(z - z_{0h}) (\bar{u} - \bar{u}_{z0})^2} \quad (1.107)$$

Le nombre de Richardson volumique est souvent utilisé pour la détermination de l'indice de Monin et Obukhov. Il est toutefois à souligner qu'il existe dans la littérature différentes expressions du nombre de Richardson volumique avec ou sans les rugosités dynamiques et
thermiques. Enfin, le critère de passage d'un régime de convection à l'autre est [Guyot 97] :

$R_i > -0.03$ : Convection forcée

$-1 \le R_i \le -0.03$ : Convection mixte

$R_i < -1$ : Convection libre

## 1.5.2 La théorie de similitude de Monin et Obukhov

L'analyse dimensionnelle effectuée par Monin et Obukhov (1954) dans la couche de surface fait ressortir cinq paramètres fondamentaux dont la combinaison à permis de définir des échelles caractéristiques.

- Les paramètres qui apparaissent dans les équations fondamentales, à savoir la masse volumique $\rho$ et le coefficient de flottabilité $\beta = \frac{g}{T_r}$;

- Le flux de chaleur au sol $Q_0$ et le cisaillement $\tau_0$ qui caractérisent l'interaction de l'écoulement avec la surface;

- La hauteur $z - d$, $d$ étant ajoutée afin de traiter les surfaces recouvertes par un couvert végétal. La longueur de rugosité $z_0$ est exclue de cette liste en considérant que $z \gg z_0$, et que dans ces conditions, $z_0$ n'affecte pas la variation des champs moyens ni les flux turbulents [DeMoor 83, Guyot 97].

De ces paramètres, il en ressort :

- Une échelle de vitesse $u_*$ et de température $T_*$;

- Une échelle de longueur $L$, qui définit la longueur de Monin et Obukhov et désigne la stratification de l'atmosphère :

$$L = - \frac{u_*^3}{\kappa \beta (T'w')} = - \frac{u_*^3}{\kappa \beta} \left( \frac{\rho C_P}{Q_0} \right) = \frac{u_*^2}{\kappa \beta T_*} \quad (1.108)$$

Des trois échelles, Monin et Obukhov introduisent une variable indépendante sans dimension appelée indice de Monin-Obukhov et définie par :

$$\xi = \frac{z-d}{L} \quad (1.109)$$
La valeur de ξ va caractériser la stabilité de la couche limite atmosphérique :

$\xi > 0$ (ou $L > 0$) Couche stable

$\xi = 0$ (ou $L \to \infty$) Couche neutre

$\xi < 0$ (ou $L < 0$) Couche instable

L'hypothèse de similitude s'appuie sur le principe que toutes les dérivées locales verticales des variables moyennes turbulentes adimensionnées par les échelles (ou une combinaison) $u_*$, $T_*$, $L$ ne sont fonction que de l'indice $\xi$. Ce dernier point est l'essence même de la théorie de Monin et Obukhov, puisqu'elle considère que les caractéristiques des processus turbulents ne dépendent plus que de $\xi$. Le principe de similitude s'écrit :

$$ \frac{M(z)}{M_*} = f_M(\xi) \qquad (1.110) $$

$f_M(\xi)$ est une fonction universelle de $\xi$.

### 1.5.3 Relations de similitude dans la couche de surface

Les fonctions universelles de similitude tiennent leur nom du fait qu'elles sont les mêmes, quelle que soit la couche limite de surface. Elles sont déterminées expérimentalement en s'appuyant sur leurs comportements asymptotiques qui sont relativement bien connus. La théorie de similitude (1.110) mène alors à des relations flux-gradient entre les flux turbulents et les gradients des grandeurs moyennes dans la couche de surface. Ces relations sont d'un grand intérêt pour la compréhension des phénomènes qui interviennent dans la couche limite.

Les flux turbulents et la longueur de Monin et Obukhov au sein de cette couche sont constants. Les gradients s'expriment alors en fonction de l'indice de Monin-Obukhov à travers les fonctions de similitude :

$$ \frac{\partial \bar{u}}{\partial z} = \frac{u_*}{\kappa (z-d)} \phi_m \left( \frac{z-d}{L} \right) \qquad (1.111) $$

$$ \frac{\partial \bar{T}}{\partial z} = \frac{T_*}{\kappa (z-d)} \phi_h \left( \frac{z-d}{L} \right) \qquad (1.112) $$

En intégrant ces relations flux-gradient entre $z-d$ et $z_0$ ($z_{0h}$ pour la température), les
profils corrigés du vent horizontal moyen et de la température statique sont donnés par :

$$ \bar{u}(z) = \frac{u_*}{\kappa} \left[ \ln\left(\frac{z-d}{z_0}\right) - \psi_m\left(\frac{z-d}{L}\right) + \psi_m\left(\frac{z_0}{L}\right) \right] \quad (1.113) $$

$$ \bar{T}(z) - T_s = \frac{P_{rt}T_*}{\kappa} \left[ \ln\left(\frac{z-d}{z_{0h}}\right) - \psi_h\left(\frac{z-d}{L}\right) + \psi_h\left(\frac{z_{0h}}{L}\right) \right] \quad (1.114) $$

Les fonctions de similitude $\psi_m$ et $\psi_h$ étant définies par :

$$ \psi_m(\xi) = \int_{z_0}^{z-d} \frac{1 - \phi_m(\xi)}{\xi} d\xi $$

$$ \psi_h(\xi) = \int_{z_{0h}}^{z-d} \frac{1 - \phi_h(\xi) / P_{rt}}{\xi} d\xi $$

Une fonction de similitude associée à la dissipation $\epsilon$, $\phi_\epsilon$ est introduite, puisque la stabilité thermique modifie l'équilibre local à proximité des parois, où la dissipation doit contrebalancer la production cinétique et la production gravitationnelle [Kaimal 94, Duynkerke 88] :

$$ \epsilon = \frac{u_*^3}{\kappa(z-d)} \phi_\epsilon \left(\frac{z-d}{L}\right) \quad (1.115) $$

Les fonctions de similitude et les nombres de Richardson $R_f$ et $R_i$ sont reliés par ces expressions :

$$ R_f = \xi \frac{k_t \phi_h(\xi)}{\nu_t \phi_m^2(\xi)} = \frac{\xi}{\phi_m(\xi)} \quad (1.116) $$

$$ R_i = \frac{\nu_t}{k_t} R_f = \xi \frac{\phi_h(\xi)}{\phi_m^2(\xi)} \quad (1.117) $$

La diffusivité thermique turbulente est définie à travers le nombre de Prandtl turbulent dont l'expression est donnée en fonction de $\xi$ par :

$$ P_{rt} = \frac{\nu_t}{k_t} = \frac{\phi_m(\xi)}{\phi_h(\xi)} \quad (1.118) $$

### 1.5.4 Les fonctions universelles de Dyer

La détermination des fonctions de similitude a fait l'objet de nombreuses études expérimentales. En se basant sur un grand nombre de données obtenues sur sol nu ou prairie,
[Dyer 74] propose des expressions dont la précision semble suffisante dans la majorité des cas :

- La fonction $\phi_m$ pour les flux turbulents de quantité de mouvement :

$$ \phi_m \left( \frac{z-d}{L} \right) = \left( 1 + 16 \left| \frac{z-d}{L} \right| \right)^{-\frac{1}{4}} \quad -2 \le \frac{z-d}{L} \le 0 $$

$$ \phi_m \left( \frac{z-d}{L} \right) = 1 + 5 \left( \frac{z-d}{L} \right) \quad 0 \le \frac{z-d}{L} \le 1 $$

- La fonction $\phi_h$ de correction des flux turbulents de chaleur :

$$ \phi_h \left( \frac{z-d}{L} \right) = \left( 1 + 16 \left| \frac{z-d}{L} \right| \right)^{-\frac{1}{2}} \quad -2 \le \frac{z-d}{L} \le 0 $$

$$ \phi_h \left( \frac{z-d}{L} \right) = 1 + 5 \left( \frac{z-d}{L} \right) \quad 0 \le \frac{z-d}{L} \le 1 $$

- La fonction $\phi_\epsilon$ pour le taux de dissipation :

$$ \phi_\epsilon \left( \frac{z-d}{L} \right) = \left( 1 + 0.5 \left| \frac{z-d}{L} \right|^{\frac{2}{3}} \right)^{\frac{3}{2}} \quad -2 \le \frac{z-d}{L} \le 0 $$

$$ \phi_\epsilon \left( \frac{z-d}{L} \right) = 1 + 5 \left( \frac{z-d}{L} \right) \quad 0 \le \frac{z-d}{L} \le 1 $$

Les fonctions de correction des profils turbulents du vent et de la température sont données par l'intégration de $\phi_m$ et $\phi_h$ et s'écrivent :

- La fonction $\psi_m$ pour la vitesse horizontale moyenne :

$$ \psi_m(\xi) = \ln \frac{[1 + x(\xi)^2][1 + x(\xi)]^2}{8} - 2 \arctan[x(\xi)] + \frac{\pi}{2} \quad -2 \le \frac{z-d}{L} \le 0 $$

$$ \psi_m(\xi) = -5\xi \quad 0 \le \frac{z-d}{L} \le 1 $$

- La fonction $\psi_h$ pour la température moyenne :

$$ \psi_h(\xi) = 2 \ln \frac{[1 + y(\xi)]}{2} \quad -2 \le \frac{z-d}{L} \le 0 $$

$$ \psi_h(\xi) = -5\xi \quad 0 \le \frac{z-d}{L} \le 1 $$

où les fonctions $x(\xi)$ et $y(\xi)$ sont :

$$ x(\xi) = (1 + 16 |\xi|)^{\frac{1}{4}} \quad \text{et} \quad y(\xi) = (1 + 16 |\xi|)^{\frac{1}{2}} $$
Lorsque la couche de surface tend vers une stabilité neutre ou quasi-neutre (typiquement en régime de convection forcée), le comportement asymptotique des fonctions de similitude est alors :

$$ \phi_m (0) = 1 $$

$$ \phi_h (0) = P_{rt} $$

$$ \phi_{\epsilon} (0) = 1 $$

Ceci permet de retrouver les relations dans le cas d'un écoulement en régime de convection forcée ou en stratification neutre, établies précédemment.

La stabilité thermique de l'atmosphère est intégrée au modèle dynamique et turbulent à l'aide de la théorie de similitude de Monin et Obukhov. Elle permet également de définir une diffusivité thermique turbulente. Dans la couche de surface en stratification thermique, les relations flux-gradient et les profils sont corrigés par des fonctions de similitude. Nous utiliserons les fonctions de correction proposées par Dyer pour généraliser les conditions aux limites au sol et au sommet de la couche de surface. Nous limiterons toutefois l'étude pour des écoulements en faible stratification thermique, pour des valeurs de ξ comprises entre -2 et 1, en accord avec la plage définie par [Dyer 74].

## 1.6 Modélisation de l'écoulement à l'échelle d'un couvert végétal

La libération, le transport et le dépôt des particules biotiques représentent un processus d'échange entre l'atmosphère et la végétation. Il est donc important de comprendre les transferts turbulents à l'échelle d'un couvert végétal, et cela dans le but d'intégrer des informations à plus grande échelle. Nous avons donc «élargi» notre discussion à une échelle plus locale, celle d'un couvert végétal continu. Modéliser le vent dans le couvert végétal passe par une description fine du mouvement et de la turbulence dans la sous-couche rugueuse, que nous n'avons pas considérée jusqu'à présent.

Dans la sous-couche rugueuse (figure 1.7), qui s'étend du sol à une distance $z_* = 2h$ [Raupach 96], les caractéristiques de l'écoulement s'écartent de la couche inertielle. En effet, l'écoulement à cette échelle est largement influencé par la géométrie des éléments présents dans la sous-couche rugueuse. Les relations de similitude établies précédemment
<figure><img src="image_13.png" /><figcaption>Figure 1.7 : Profil schématique du vent moyen au-dessus et à l'intérieur d'un couvert végétal plat et homogène [Ogée 00]</figcaption></figure>

sont alors mises en défaut, et les flux turbulents de quantité de mouvement, de chaleur et de masse ne sont plus constants.

Nous allons présenter le modèle que nous avons mis en place pour la modélisation de l'écoulement à cette échelle, en nous limitant à un écoulement en stratification neutre.

1.6.1 La turbulence et le vent

Un couvert végétal est composé de nombreux éléments tels que des feuilles, des tiges et des branches. L'ensemble peut alors former dans l'espace une structure complexe. Lorsque ces éléments interagissent avec l'écoulement d'air, de nombreux processus aérodynamiques ont lieu :

- Le flux de quantité de mouvement n'est plus conservatif car il est absorbé, par pression et viscosité, par le couvert végétal. Au-dessus du couvert, le flux adimensionné (figure 1.8) est constant et décroît rapidement au fur et à mesure qu'on pénètre dans le couvert, pour atteindre pratiquement une valeur nulle, près du sol;

- Les éléments de la végétation répondent à la sollicitation mécanique du vent en générant des sillages (wake) turbulents. Ces sillages vont favoriser le transfert de l'énergie cinétique du mouvement moyen au profit de l'énergie cinétique du mouvement turbulent, dont la longueur caractéristique est celle des éléments végétaux [Raupach 81];
<figure><img src="image_14.png" /><figcaption>Figure 1.8 : *Familles de profils de vent et de flux turbulents de quantité de mouvement (adimensionnés) dans différents couverts végétaux* [Raupach 96]</figcaption></figure>

- La dispersion verticale des scalaires est accrue compte tenu du brassage intense qui s'exerce entre les éléments de la végétation et l'air;

- L'écoulement est fortement tridimensionnel dans et au-dessus de la canopée puisqu'il est influencé par la répartition spatiale des éléments du couvert.

Caractéristiques de la turbulence de couvert

Plusieurs études se sont consacrées à analyser la turbulence dans des couverts végétaux, que ce soit en soufflerie [Brunet 94] ou en site naturel [Raynor 71]. Ces études ont permis de dégager un certain nombre de caractéristiques propres à la turbulence de couvert mais qui se retrouvent d'un couvert à un autre, ce qui leur confère un caractère universel. La turbulence de couvert est caractérisée par de grosses structures cohérentes dont la dimension est de l'ordre de la hauteur du couvert *h*. Ces structures s'expliquent par l'existence d'un point d'inflexion dans le profil du vent moyen au voisinage du sommet du couvert (figure 1.8). L'évolution de ce point d'inflexion peut générer des instabilités dynamiques qui sous l'influence du cisaillement tendent vers ces structures caractéristiques [Raupach 96].

Ces processus de génération des larges structures cohérentes dans le couvert végétal (figure 1.9) sont très similaires à ceux de la couche de mélange (interface entre deux couches ayant une vitesse différente) près de l'interface couvert-atmosphère, ce qui lui vaut l'analogie « Canopy-Mixing layer » [Raupach 96]. Cette turbulence est très intermittente
<figure><img src="image_15.png" /><figcaption>Figure 1.9 : Étapes du développement de l'instabilité dans la sous-couche rugueuse
[Finnigan 00]</figcaption></figure>
et énergique et assure une part importante des transferts turbulents.

La compréhension des processus qui apparaissent lors de l'interaction de l'écoulement avec les éléments végétaux passe également par la considération des structures fines de la turbulence de couvert. En effet, comme nous le verrons par la suite, ces interactions vont modifier le bilan de l'énergie cinétique moyenne du mouvement turbulent, la forme de son spectre, la dissipation visqueuse ainsi que la cascade énergétique des tourbillons décrite par Kolmogorov [Finnigan 00].

## Profil du vent moyen dans le couvert

Les premières modélisations du vent dans un couvert végétal s'appuient sur la théorie de longueur de mélange de Prandtl-Von Karman [Raupach 81] empruntée à la théorie des couches limites. Soit $F$ la force aérodynamique qui s'exerce dans une tranche de végétation :

$$F(z) = \frac{1}{2} \rho C_d a(z) \bar{u}^2(z) \qquad (1.119)$$

Cette force volumique qui n'est autre que la force de traînée, dépend du coefficient de traînée effectif $C_d$ des éléments végétaux constituant la tranche de végétation et de la surface foliaire frontale (perpendiculaire à l'écoulement) par unité de volume de couvert $a(z)$ (une seule face est considérée). $F$ est reliée au flux net de quantité de mouvement $\tau$ par :

$$\frac{\partial \tau}{\partial z} = F(z) = \frac{1}{2} \rho C_d a(z) \bar{u}^2(z) \qquad (1.120)$$

Or, le flux de quantité de mouvement peut s'exprimer en fonction du flux turbulent, si l'on néglige le flux moléculaire :

$$-\rho \overline{u'w'} = \tau = \mu_t \frac{\partial \overline{u}}{\partial z} \qquad (1.121)$$

En remplaçant (1.121) dans l'équation (1.120) et en dérivant, nous obtenons une équation différentielle du second ordre dont l'inconnue est la vitesse moyenne $\bar{u}(z)$ :

$$\frac{\partial \mu_t}{\partial z} \frac{\partial \bar{u}}{\partial z} + \mu_t \frac{\partial^2 \bar{u}}{\partial z^2} = \frac{1}{2} \rho C_d a(z) \bar{u}^2(z) \qquad (1.122)$$

La viscosité turbulente est modélisée à l'aide de la longueur de mélange $l_m$ caractéristique du niveau $z$ dans le couvert végétal. Cette équation admet une solution analytique dans le cas d'un couvert uniforme ($a(z) = \text{constante}$ et $l_m$ indépendant de $z$) et dont le coefficient de traînée est constant. Dans ce cas particulier, la vitesse moyenne dans le couvert végétal
suit une loi exponentielle donnée par :

$$
\bar{u}(z) = \bar{u}(h) \exp \left[ \alpha \left( \frac{z}{h} - 1 \right) \right] \qquad (1.123)
$$

$h$ représente la hauteur du couvert et $\alpha$ est le coefficient d'atténuation [Raupach 81] du vent qui dépend du coefficient de traînée $C_d$, de la densité de surface végétale frontale $a$ et de la longueur de mélange dans le couvert.

$$
\alpha = \left( \frac{C_d a h^3}{2 l_m^2} \right)^{1/3} \qquad (1.124)
$$

Cette loi de variation de la vitesse horizontale, bien que basée sur des hypothèses fortement restrictives, permet de reproduire le vent dans les couches supérieures d'un couvert de maïs de façon assez satisfaisante [Loubet 00].

Ces modèles restreints aux couverts végétaux homogènes de type cultures, n’offrent pas de possibilités d’applications dans les cas de couverts hétérogènes et discontinus (forêt, couverts en rangs…). De plus ils n’apportent pas ou peu d’informations sur les caractéristiques de la turbulence de couvert. Tous les efforts vont dans ce sens, et se tournent vers une description du vent et de la turbulence par les équations de conservation de quantités de mouvement et par la modélisation du mouvement turbulent. Nous allons présenter les hypothèses supplémentaires à apporter dans le système de la couche limite atmosphérique pour la modélisation du vent dans les couverts végétaux ainsi que le modèle retenu pour la couche limite en stratification neutre.

1.6.2 Moyenne spatio-temporelle

Les variables turbulentes varient fortement dans l'espace car les éléments végétaux
jouent le rôle de sources/puits de matière et de quantité de mouvement. Ceci conduit à
introduire en plus de la moyenne temporelle classique (décomposition de Reynolds équa-
tion (1.10)) une moyenne spatiale pour décrire le microclimat moyen autour des éléments
végétaux [Raupach 82, Finnigan 85, Li 85, Brunet 94].

Soit $f$ une variable turbulente qui admet une décomposition de Reynolds. La tur-
bulence dans un couvert végétal ne vérifie plus le principe d'ergodicité, et la moyenne
temporelle $\bar{f}$ s'écrit à son tour comme étant la somme d'une moyenne spatiale $\langle \bar{f} \rangle$ qui est
la moyenne spatiale de la moyenne temporelle, et de sa fluctuation (écart à cette moyenne)
notée $\bar{f}''$ :

$$
\bar{f} = \langle \bar{f} \rangle + \bar{f}'' \tag{1.125}
$$
L'utilisation d'une moyenne spatiale en plus de la moyenne temporelle est assez récente. Elle a été introduite par Wilson & Shaw (1977) et Raupach & Shaw (1982) [Raupach 82]. Depuis, la moyenne spatio-temporelle est de plus en plus usitée. D'un point de vue mathématique, la définition formelle de la moyenne spatiale est [Raupach 82]:

$$ \langle \phi \rangle (\mathbf{x}, t) = \frac{1}{V} \iiint_V \phi (\mathbf{x} + \mathbf{r}, t) \, d\mathbf{r} \qquad (1.126) $$

Où $\phi$ est un champ turbulent défini autour des éléments végétaux. V est le volume d'intégration duquel sont exclus les éléments végétaux. Il doit être pris suffisamment large suivant le plan horizontal afin d'éliminer les variations de la structure du couvert, mais suffisamment petit suivant $z$ pour préserver les caractéristiques des variations des grandeurs dans la canopée. V est multiplement connexe du fait de la présence des éléments du couvert (figure 1.10).

<figure><img src="image_16.png" /><figcaption>Figure 1.10 : Vue schématique du volume d'intégration dans une forêt. Les parties solides du tronc sont exclues du volume d'intégration et font que V est multiplement connexe [Kaimal 94]</figcaption></figure>

Contrairement à l'opérateur de Reynolds pour lequel la moyenne temporelle et l'opérateur différentiation commutent –ceci du fait du principe d'ergodicité que vérifient les grandeurs turbulentes– l'opérateur moyenne spatiale ne commute pas, sauf si la variable considérée est constante à l'interface air-organes végétaux [Raupach 82].

$$ \left\langle \frac{\partial \phi}{\partial x_i} \right\rangle = \frac{\partial \langle \phi \rangle}{\partial x_i} - \frac{1}{V} \iint_{S_I} \phi n_i dS \qquad (1.127) $$

$S_I$ est la surface qui coïncide avec la surface des éléments végétaux, et $n_i$ est le vecteur normal des éléments végétaux.
La fluctuation spatiale vérifie la propriété suivante :

$$
\langle \phi'' \rangle = 0 \tag{1.128}
$$

Si nous posons $\phi = \bar{f}$ dans l'équation (1.127), la moyenne spatiale de la dérivée de la
moyenne temporelle $\bar{f}$ s'écrit :

$$
\left\langle \frac{\partial \bar{f}}{\partial x_i} \right\rangle = \frac{\partial \langle \bar{f} \rangle}{\partial x_i} - \frac{1}{V} \iint_{S_I} \bar{f} \, n_i dS \quad (1.129)
$$

qui doit être rigoureusement identique à (en décomposant $\bar{f}$ en une partie moyenne et une
partie fluctuante) :

$$
\left\langle \frac{\partial (\langle \overline{f} \rangle + \overline{f}'' )}{\partial x_i} \right\rangle = \left\langle \frac{\partial \langle \overline{f} \rangle}{\partial x_i} \right\rangle + \left\langle \frac{\partial \overline{f}''}{\partial x_i} \right\rangle \quad (1.130)
$$

Par conséquent, nous pouvons en déduire que :

$$
\left\langle \frac{\partial \langle \overline{f} \rangle}{\partial x_i} \right\rangle = \frac{\partial \langle \overline{f} \rangle}{\partial x_i} \qquad (1.131)
$$

$$
\left\langle \frac{\partial \overline{f}''}{\partial x_i} \right\rangle = -\frac{1}{V} \iint_{S_I} \overline{f} \, n_i dS \quad (1.132)
$$

L'opérateur moyenne spatial et l'opérateur différentiation sont donc non-commutatifs, et
nous notons en particulier que :

$$
\left\langle \frac{\partial \overline{f}''}{\partial x_i} \right\rangle \neq \frac{\partial \left\langle \overline{f}'' \right\rangle}{\partial x_i} \quad \text{puisque} \quad \left\langle \overline{f}'' \right\rangle = 0 \qquad (1.133)
$$

1.6.3 Équations du mouvement dans le couvert végétal

Le vent instantané et la pression dans le couvert admettent une décomposition spatio-
temporelle, et les grandeurs instantanées turbulentes sont données par :

$$
u_i = \langle \overline{u_i} \rangle + \overline{u_i}'' + u_i'
$$

$$
P = \langle \overline{P} \rangle + \overline{P}'' + P'
$$

En moyennant les équations instantanées avec l'opérateur de Reynolds (moyenne tem-
porelle) puis en décomposant les variables turbulentes suivant (1.125), et en moyennant
une seconde fois avec la moyenne spatiale$^6$, les équations de conservation de masse et de quantité de mouvement dans le couvert végétal s'écrivent :

$$ \frac{\partial \langle \overline{u}_j \rangle}{\partial x_j} = 0 \qquad (1.134) $$

$$ \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 \langle \overline{P} \rangle}{\partial x_i} + \underbrace{\frac{\partial}{\partial x_j} \left[ \mu \left( \frac{\partial \langle \overline{u}_i \rangle}{\partial x_j} + \frac{\partial \langle \overline{u}_j \rangle}{\partial x_i} \right) \right] - \frac{\partial}{\partial x_j} \left[ \rho \langle \overline{u}_i' \overline{u}_j' \rangle + \rho \langle \overline{u}_i'' \overline{u}_j'' \rangle \right]}_{\text{I}} - \underbrace{\left\langle \frac{\partial \overline{P}''}{\partial x_i} \right\rangle + \left\langle \frac{\partial}{\partial x_j} \left[ \mu \left( \frac{\partial \overline{u}_i''}{\partial x_j} + \frac{\partial \overline{u}_j''}{\partial x_i} \right) \right] \right\rangle}_{\text{II}} \qquad (1.135) $$

-I- Le tenseur des contraintes se compose du tenseur des contraintes visqueuses, du tenseur de Reynolds et d'un tenseur supplémentaire résultant de la corrélation spatiale des quantités ramenées à la moyenne temporelle qui dépendent de la position. Ces termes de corrélations, ou extra tenseur de Reynolds, forment le tenseur dispersif;

-II- La présence de ces termes sont une conséquence directe de la non-commutativité des opérateurs de moyenne et de différentiation spatiale. Compte tenu de l'équation (1.132), ces termes s'écrivent [Raupach 86, Brunet 94]:

$$ f_{F_i} = - \left\langle \frac{\partial \overline{P}''}{\partial x_i} \right\rangle = \frac{1}{V} \iint_{S_I} \overline{P} n_i dS \qquad (1.136) $$

$$ f_{V_i} = \left\langle \frac{\partial}{\partial x_j} \left[ \mu \left( \frac{\partial \overline{u}_i''}{\partial x_j} + \frac{\partial \overline{u}_j''}{\partial x_i} \right) \right] \right\rangle = - \frac{\mu}{V} \iint_{S_I} \frac{\partial \overline{u}_i}{\partial n} dS \qquad (1.137) $$

$f_{F_i}$ et $f_{V_i}$ sont respectivement les forces de traînée de forme et de traînée de frottement qui s'exercent dans le volume d'intégration [Brunet 94]. Elles sont modélisées par :

$$ F_{d_i} = f_{F_i} + f_{V_i} = -\rho C_{d_i} A_{L_i} \sqrt{\langle \overline{u}_j \rangle \langle \overline{u}_j \rangle} \langle \overline{u}_i \rangle \qquad (1.138) $$

$F_{d_i}$ est la force volumique totale due à l'absorption (par pression et viscosité) de quantité de mouvement par les éléments du couvert de densité de surface végétale

<footnote n="6">6. Les calculs sont présentés en annexe A.5.1</footnote>
$A_{L_i}$ qui représente la surface des végétaux (feuilles, tige ...) dans un volume élémentaire autour du point considéré. Par ailleurs, à l'échelle d'une feuille élémentaire, $f_{F_i}$ s'exerce sur la surface normale à l'écoulement (surface frontale) alors que $f_{V_i}$ va s'exercer sur la surface mouillée par l'écoulement (surface totale de la feuille) [Kaimal 94, Finnigan 00].

Ces considérations rendent cette modélisation difficile et inutilisable car il n'est pas possible d'intégrer ces forces sur chaque élément végétal. Cependant, les nombres de Reynolds typiques des écoulements dans la partie supérieure du couvert où la quantité de mouvement est principalement absorbée, montrent que les forces de traînées dues à la pression sont 4 fois supérieures à la traînée de frottement [Kaimal 94]. Finalement la force volumique $F_{d_i}$ est modélisée par :

$$F_{d_i} = -\rho C_d a \sqrt{\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle \langle \overline{u_i} \rangle} \quad (1.139)$$

$F_{d_x}$ représente alors la force de traînée, $F_{d_y}$ et $F_{d_z}$ les composantes de la portance. $C_d$ est le coefficient de traînée effectif au niveau $z$ des éléments végétaux. Il dépend de l'orientation des feuilles, de leurs mutuelles interférences ainsi que du nombre de Reynolds de l'écoulement [Brutsaert 82]. De façon général, il est plus faible que celui d'un élément isolé dans l'espace avec un écoulement contrôlé [Ogée 00]. Enfin, $a$ [$m^{-1}$] est la densité de surface végétale normale à l'écoulement.

Ici encore, le problème de fermeture se pose. L'utilisation de la moyenne spatiale temporelle fait apparaître des termes de corrélation double $\rho \langle u'_i u'_j \rangle$ et $\rho \langle \overline{u_i''} \overline{u_j''} \rangle$ qui représentent respectivement les flux turbulents et les flux dispersifs dans la végétation. Ces flux dispersifs sont néanmoins négligeables devant les flux turbulents, en particulier vers le haut du couvert végétal [Raupach 82, Liu 96].

Sous ces conditions, les flux dispersifs et les flux turbulents sont modélisés par le tenseur de Reynolds, soit :

$$\rho \langle \overline{u'_i u'_j} \rangle + \rho \langle \overline{u_i'' \overline{u_j''}} \rangle = -\mu_t \left[ \frac{\partial \langle \overline{u_i} \rangle}{\partial x_j} + \frac{\partial \langle \overline{u_j} \rangle}{\partial x_i} \right] + \frac{2}{3} \rho \delta_{ij} k \quad (1.140)$$

La modélisation phénoménologique des flux turbulents par cette relation (1.140) dans le cas du couvert végétal est très critiquable pour des raisons détaillées au chapitre 2 section 2.5.2.

Enfin, la viscosité turbulente est évaluée à l'aide de la connaissance de l'énergie cinétique $k$ et du taux de dissipation $\varepsilon$ (par simplification, nous omettrons les crochets) qui
s'écrivent compte tenu de la moyenne spatio-temporelle :

$$k = \frac{1}{2} \langle \overline{u'_i u'_i} \rangle \qquad (1.141)$$

$$\varepsilon = \nu \left\langle \frac{\partial u'_i}{\partial x_j} \frac{\partial u'_i}{\partial x_j} \right\rangle \qquad (1.142)$$

Il est donc nécessaire d'expliciter l'équation d'évolution de l'énergie cinétique moyenne du mouvement turbulent et de sa dissipation dans le couvert végétal.

### 1.6.4 Modélisation du mouvement turbulent dans le couvert

En plus de la force de traînée, l'interaction aérodynamique avec les éléments végétaux induit une transformation profonde de la structure de la turbulence. Ces effets vont agir sur le mouvement turbulent et cela en modifiant le bilan de l'énergie cinétique moyenne du mouvement turbulent $k$, comme nous pouvons le voir sur la figure 1.11. Les diverses contributions des termes de $k$ se distinguent de façon claire :

- Pour $z \approx 2h$ et au-delà, les termes dominants dans l'équation d'évolution de $k$ sont la production par cisaillement et le taux de dissipation. La turbulence est en équilibre locale ($P_k = \varepsilon$), et l'énergie cinétique moyenne du mouvement turbulent est transférée depuis les grandes échelles vers les plus petites. Ces structures de dissipation sont caractérisées par la micro-échelle de longueur de Kolmogorov :

$$\eta = \left( \frac{\nu^3}{\varepsilon} \right)^{\frac{1}{4}} \qquad (1.143)$$

Cette caractéristique de la turbulence est la signature de la couche inertielle;

- Pour $h < z < 2h$, le transport turbulent $T_t$ de $k$ augmente de manière significative au fur et à mesure que $z$ se rapproche de $h$. Ce terme est négatif. Le transport par pression $T_p$ fluctuante intervient également dans le bilan de $k$, cependant il est difficilement mesurable et sujet à caution. La production n'est plus égale à la dissipation et laisse entendre que la majeure partie de la turbulence n'est pas localement produite mais transportée [Raupach 96];

- Pour $h < z$ le transport turbulent de $k$ reste important, mais correspond à un gain (positif). Enfin, le bilan de $k$ est largement dominé par le terme de production de sillage $P_w$, qui dépasse la production par cisaillement. Cette turbulence est rapidement dissipée causant une extra-dissipation. De nouveau l'équilibre local n'est
pas vérifiée. C'est la principale caractéristique de la turbulence dans la sous-couche rugueuse.

<figure><img src="image_17.png" /><figcaption>Figure 1.11 : Contribution des différents termes de l'équation de l'énergie cinétique moyenne du mouvement turbulent dans et au-dessus de la sous-couche rugueuse pour une maquette de couvert en soufflerie [Brunet 94]. $P_s$ = production par cisaillement, $\varepsilon$ = taux de dissipation, $P_w$ = production par sillage, $T_t$ = transport par le mouvement turbulent et $T_p$ = transport par pression fluctuante</figcaption></figure>

Explicitons les équations de transport de $k$ et de $\varepsilon$. Le traitement statistique par la moyenne spatio-temporelle de l'équation de transport du mouvement fluctuant, multipliée par la vitesse fluctuante, permet d'écrire l'équation de transport de l'énergie cinétique moyenne du mouvement$^7$. Cette équation (1.144) se présente sous sa forme habituelle à l'exception d'un terme supplémentaire qui apparaît après le double traitement statistique. Par conséquent, cette modification se répercute sur l'équation (1.145) du taux de dissipation.
Les équations de $k$ et $\varepsilon$ sont données dans le couvert végétal par :

$$ \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 + \underbrace{\rho S_k}_{I} \quad (1.144) $$

$$ \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} + \underbrace{\rho S_\varepsilon}_{II} \quad (1.145) $$

-I- $S_k$ modifie l'équation d'évolution de $k$ par le biais d'un terme source et d'un terme puits :

$$ S_k = P_w - \varepsilon_w $$

Le terme source $P_w$ est le terme de production de sillage de $k$ (wake production). Il traduit une accentuation de la transformation de l'énergie cinétique du mouvement moyen où $P_w$ apparaît avec le signe "+" au profit de l'énergie cinétique moyenne du mouvement turbulent, où il apparaît avec le signe "-" induite par le couvert [Raupach 82]. La production de sillage est conventionellement modélisée comme le produit du vent moyen et de la force de traînée [Green 92] :

$$ P_w = - \left\langle \overline{u'_i u''_j} \frac{\partial \overline{u''_i}}{\partial x_j} \right\rangle = C_d a \left[ \langle \overline{u_j} \rangle \langle \overline{u_j} \rangle \right]^{3/2} \quad (1.146) $$

En même temps, le processus de la cascade énergétique de Kolmogorov des larges échelles vers les petites échelles est accéléré (figure 1.12) par la présence des éléments végétaux, qui activent la dissipation de l'énergie par frottement visqueux. Ce terme puits $\varepsilon_w$, qui est le taux de dissipation de sillage, est proportionnel au vent moyen et à l'énergie cinétique $k$ [Green 92, Wilson 96] :

$$ \varepsilon_w = 4C_d a \sqrt{\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle} k \quad (1.147) $$

-II- L'équation de bilan de $\varepsilon$ est modifiée par le terme $S_\varepsilon$ :

$$ S_\varepsilon = P_{\varepsilon w} - D_{\varepsilon w} $$

Il est bien évident que l'interprétation physique de tous les termes de l'équation exacte de $\varepsilon$ (notamment dans le couvert végétal) présente de nombreuses difficultés. L'expression de $S_\varepsilon$ est ainsi déduite de la relation de Kolmogorov qui relie $\varepsilon$ à $k$ et en dérivant, il en résulte :

$$ \varepsilon = C_{l\varepsilon} \frac{k^{\frac{3}{2}}}{l_m} \Rightarrow d\varepsilon = \frac{3}{2} \frac{\varepsilon}{k} dk $$
Cette relation, qui considère que la longueur de mélange est localement invariante, permet d'évaluer de façon empirique la production de sillage $P_{\varepsilon_w}$ et la dissipation de sillage $D_{\varepsilon_w}$ du taux de dissipation [Liu 96] :

$$P_{\varepsilon_w} = \frac{3}{2} \frac{\varepsilon}{k} P_w = \frac{3}{2} \frac{\varepsilon}{k} C_d a [\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle]^{3/2} \quad (1.148)$$

$$D_{\varepsilon_w} = \frac{3}{2} \frac{\varepsilon}{k} \varepsilon_w = 6 C_d a \sqrt{\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle} \quad (1.149)$$

En réalité, la longueur de mélange est loin d'être constante dans le couvert, ce qui implique que les relations (1.148) et (1.149) ne sont que des approximations. Il apparaît notamment que la constante "6" de la dissipation de sillage $D_{\varepsilon_w}$ joue un rôle capital dans l'équation du taux de dissipation [Liu 96]. En effet, cette constante qui s'appuie sur la relation de Kolmogorov, ne tient pas compte de l'augmentation du taux de dissipation due aux frottements visqueux sur les éléments végétaux, et sa valeur égalée à "6" décrit simplement une diminution du taux de dissipation qui a lieu parallèlement à celle de l'énergie cinétique $k$ dans le couvert. L'utilisation de la valeur "6" dans le modèle de turbulence entraîne une surestimation de l'énergie cinétique.

<figure><img src="image_18.png" /><figcaption>Figure 1.12 : Coupure de la cascade de Kolmogorov induite par les éléments d'un couvert végétal [Finnigan 00]</figcaption></figure>

L'analyse spectrale de $k$ permet d'étudier la répartition en fréquences de l'énergie cinétique moyenne du mouvement turbulent. La forme du spectre de $k$ (figure 1.12) permet
de visualiser le second processus de destruction de l'énergie cinétique du mouvement tur-
bulent. Cette dissipation est générée par les éléments du couvert qui sont à l'origine d'une
activité visqueuse à des échelles plus grandes que l'échelle de Kolmogorov : ces tourbillons
dissipatifs sont donc de tailles comparables et supérieures aux éléments végétaux. Par
ailleurs, cette extra-dissipation n'est pas isotrope car elle est imposée par la géométrie du
couvert [Finnigan 00]. Toutes ces remarques seront mises en valeur lors de la validation
du modèle où nous présenterons l'influence du terme puits de la dissipation $D_{\epsilon w}$ ainsi que
la valeur retenue de ce coefficient que nous noterons $C_{D_{\epsilon w}}$.

Finalement, le modèle de turbulence de canopée de type $k-\epsilon$ complété du jeu de constantes synthétisé dans le tableau 1.5 s'écrit :

$$
\begin{equation}
\begin{split}
\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 \epsilon \\
& + C_{P_w} \rho C_d a \left( \langle \overline{u_j} \rangle \langle \overline{u_j} \rangle \right)^{3/2} - C_{\epsilon w} \rho C_d a \sqrt{\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle} k
\end{split}
\tag{1.150}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
\rho \frac{\partial \epsilon}{\partial t} + \rho \langle \overline{u_j} \rangle \frac{\partial \epsilon}{\partial x_j} ={}& \frac{\partial}{\partial x_j} \left[ \left( \mu + \frac{\mu_t}{\sigma_\epsilon} \right) \frac{\partial \epsilon}{\partial x_j} \right] + C_{\epsilon 1} \rho \frac{\epsilon}{k} P_k - C_{\epsilon 2} \rho \frac{\epsilon^2}{k} \\
& + C_{P_{\epsilon w}} \rho \frac{\epsilon}{k} C_d a \left( \langle \overline{u_j} \rangle \langle \overline{u_j} \rangle \right)^{3/2} - C_{D_{\epsilon w}} \rho C_d a \sqrt{\langle \overline{u_j} \rangle \langle \overline{u_j} \rangle} \epsilon
\end{split}
\tag{1.151}
\end{equation}
$$

La viscosité est calculée par la relation classique :

$$
\mu_t = \rho C_\mu \frac{k^2}{\epsilon} \tag{1.152}
$$

<table><caption>Tableau 1.5 : *Constantes du modèle de turbulence de canopée de type k − ε* [Liu 96]</caption>
  <thead>
    <tr>
      <th>&sigma;<sub>k</sub></th>
      <th>&sigma;<sub>&epsilon;</sub></th>
      <th>C<sub>&mu;</sub></th>
      <th>C<sub>&epsilon;1</sub></th>
      <th>C<sub>&epsilon;2</sub></th>
      <th>C<sub>P<sub>w</sub></sub></th>
      <th>C<sub>&epsilon;<sub>w</sub></sub></th>
      <th>C<sub>P<sub>&epsilon;<sub>w</sub></sub></sub></th>
      <th>C<sub>D<sub>&epsilon;<sub>w</sub></sub></sub></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>1</td>
      <td>1.3</td>
      <td>0.09</td>
      <td>1.44</td>
      <td>1.92</td>
      <td>1</td>
      <td>4</td>
      <td>1.5</td>
      <td>6</td>
    </tr>
  </tbody>
</table>

Nous disposons maintenant d'un modèle qui permet de résoudre l'écoulement dans le
couvert végétal. La plupart des modèles rencontrés dans la littérature utilisent des lois de
variation du vent. Ces modèles basés sur une paramétrisation de l'écoulement turbulent,
ne peuvent pas étudier des configurations hétérogènes discontinues comme par exemple
des transitions forêt-culture. Notre modèle peut simuler ces situations, et cela en spécifiant
les paramètres des différents couverts étudiés: la hauteur du couvert h, le coefficient de
 traînée C<sub>d</sub> et la densité de surface foliaire a. Une transition couvert-sol nu est traitée en
 imposant une densité de surface foliaire nulle.
Conclusion

Ce chapitre présente l'écriture mathématique des phénomènes physiques qui
entrent en compte dans la couche limite atmosphérique. Pour cela, nous avons dis-
socié deux échelles différentes d'étude : l'échelle du paysage et l'échelle du couvert.

À l'échelle du paysage, le modèle dynamique de la couche limite atmosphérique utilise les équations de Navier-Stokes simplifiées dans le cadre des hypothèses de Boussinesq. Les effets thermiques sont pris en compte par la résolution de l'équation de conservation de l'énergie. L'écoulement étant fortement turbulent, les modèles de turbulence choisis pour fermer le système d'équations de la couche limite atmosphérique sont du type k − l ou k − ε. Afin de tenir compte des échanges entre le sol et la basse atmosphère, nous avons utilisé les propriétés de la couche à flux constant pour définir les conditions aux limites.

Dans les basses couches de l'atmosphère, la température joue un rôle important
dans les transferts turbulents. La stabilité de l'atmosphère a été traitée par la théorie
de similitude de Monin et Obukhov. Elle permet en particulier de généraliser les
conditions aux limites.

À l'échelle du couvert, la turbulence est très différente de la couche de surface où les relations flux-gradients sont mises en défaut, les flux turbulents étant absorbés par les éléments du couvert. Par suite, la turbulence de couvert est caractérisée par de grosses structures cohérentes dont la présence s'explique par l'existence d'un point d'inflexion au sommet du couvert. Cette turbulence est très intermittente et énergique.

Du fait de la présence des éléments végétaux, les variables turbulentes varient fortement dans l'espace. Par conséquent, le traitement statistique de Reynolds ne suffit plus à décrire le microclimat autour des éléments végétaux. Nous avons alors introduit un opérateur spatio-temporel pour modéliser l'écoulement dans le couvert végétal. Nous proposons notamment un modèle de turbulence de type k − ε adapté aux écoulements dans le couvert végétal. Ce modèle nécessitera une attention particulière lors de la validation.

Enfin, nous notons que le modèle de canopée ainsi défini, permet de traiter des
configurations hétérogènes telles que des transitions sol nu-couvert végétal.
<figure><img src="image_19.png" /></figure>
# Chapitre 2

## Dispersion atmosphérique des particules biotiques

### Introduction

Dans le chapitre précédent, nous avons mis en place un modèle permettant de décrire la dynamique et la thermique des écoulements atmosphériques dans la couche limite. Nous avons étendu notre modèle à l'échelle d'un couvert végétal. Cette étape nous est apparue indispensable pour l'étude de la dispersion des particules biotiques dans l'atmosphère, puisque les processus de transferts turbulents passent par des échanges végétation-atmosphère.

L'objet de ce chapitre est de mettre en évidence les phénomènes physiques qui interviennent dans la dispersion atmosphérique des particules biotiques. Dans un premier temps, nous introduirons des notions nécessaires à la compréhension de la dispersion atmosphérique. Puis, nous passerons en revue les différents modèles rencontrés dans la littérature. Enfin, nous présenterons notre modèle de dispersion atmosphérique.

Notre modèle repose sur une approche eulérienne qui s'appuie sur l'équation de conservation de la concentration des particules considérées comme un champ scalaire passif. Cette approche est satisfaisante dans la couche limite [McCartney 85]. En revanche, le modèle de dispersion doit être complété pour traiter les échanges turbulents au sein du couvert. Nous présenterons le modèle retenu pour la dispersion des particules à l'échelle locale.
<figure><img src="image_20.png" /></figure>
### 2.1 Processus de dispersion dans l'atmosphère

Le vent joue un rôle important dans la dispersion atmosphérique de nombreux organismes biologiques tels que le pollen, les spores ou les graines<fnref n="1" /> (tableau 2.1). Cette dispersion implique plusieurs évènements interdépendants :

- La libération de la particule;

- Le transport par le vent dans l'atmosphère;

- Le dépôt sur le sol ou dans le couvert végétal, avec ou sans remise en suspension.

<table><caption>Tableau 2.1 : Plage des dimensions des particules biotiques dans l'atmosphère</caption><thead><tr><td>Particules</td><td>Taille [µm]</td></tr></thead><tbody><tr><td>Pollen Myosotis Arvensis</td><td>2.5</td></tr><tr><td>Pollen de Cucurbita</td><td>200</td></tr><tr><td>Spore Fomes Annasus</td><td>5</td></tr><tr><td>Spore Uncinula Necator</td><td>32 x 15</td></tr><tr><td>Graine « ailée » Alsomitra Macrocarpa</td><td>120000</td></tr></tbody></table>

Comme le suggère [Aylor 78], la dispersion des particules dans l'atmosphère peut se résumer par ces équations :

$$ \bar{C}(x_i, t) = Q(t) T(x_i, t) \qquad (2.1) $$

$$ F_d(x_i, t) = v_d(x_i, t) \bar{C}(x_i, t) \qquad (2.2) $$

La première équation exprime que la concentration locale moyenne $\bar{C}(x_i, t)$ [particule $m^{-3}$] dans l'atmosphère est directement reliée au nombre de particules émises à chaque instant par une source $Q(t)$ [particule $s^{-1}$] et aux caractéristiques du transport dans l'atmosphère représentées ici par $T(x_i, t)$ [$s m^{-3}$] ($T$ exprime la répartition de la concentration selon une courbe gaussienne par effet de dilution et de diffusion turbulente (équation (2.12))). En même temps, la seconde équation émet que le nombre de particules déposées par unité de surface et de temps à une distance $x$ de la source et à un instant $t$, $F_d(x_i, t)$

<footnote n="1">1. Complément sur les particules citées en annexe B.3</footnote>
[particule $m^{-2} s^{-1}$] est connue si la concentration ainsi que la vitesse de dépôt $v_d(x_i, t)$ [$m s^{-1}$] le sont. Cette mathématisation simple permet de synthétiser le problème de la dispersion atmosphérique et laisse en même temps entrevoir la difficulté d'appréhender tous ces évènements qui sont très corrélés. L'objectif principal des modèles de dispersion est d'apporter des éléments de compréhension sur le transport proprement dit, en analysant et en intégrant les paramètres fondamentaux qui tendent à influencer la diffusion atmosphérique.

### 2.1.1 La libération des particules biotiques

Bagnold (1941) a été un des précurseurs dans l'étude des phénomènes de mise en suspension dans l'air pour le sable. Il a mis en évidence que pour une contrainte de cisaillement supérieure à une valeur critique, les particules sont entraînées par le fluide [Bagnold 41]. Il a donc introduit la notion de seuil de transport ou de mise en mouvement dont le critère s'appuie sur une vitesse de frottement critique qui s'exprime en fonction des caractéristiques de la particule et de la vitesse de frottement du vent :

$$u_{*t} = A \sqrt{\frac{\rho_p - \rho}{\rho}} d_p g \quad (2.3)$$

$d_p$ est le diamètre moyen de la particule, $\rho_p$ sa masse volumique et $A$ une constante définie comme étant la racine carrée de la fonction de Shields [Nordstrom 90], égale à 0.1 pour un nombre de Reynolds particulière, basé sur la vitesse de frottement, supérieur à $R_{ep} > 3.5$.

<figure><img src="image_21.png" /><figcaption>Figure 2.1 : Deux mécanismes différents de libération sèche des particules dans l'atmosphère [Aylor 90]</figcaption></figure>

Pour les particules biotiques, le processus qui permet à la particule d'être arrachée de sa surface émettrice et d'être mise en suspension dans l'atmosphère est appelé mécanisme
de libération. Cette libération peut être active (figure 2.1a), indépendante ou faiblement dépendante des facteurs climatiques, ou, à l'opposé passive (figure 2.1b) dans le cas où les facteurs climatiques sont à l'origine de la libération –on parle alors d'induction de la libération. C'est l'observation de pics de concentration régulièrement atteints en milieu de journée, lorsque l'écoulement de l'air est très turbulent, qui a mis en évidence l'influence du vent dans la libération de certaines particules. Cette observation permet donc de confirmer que le vent est l'agent principal dans le mécanisme de libération. Cependant, les tentatives pour évaluer le débit Q de la source émettrice par une mesure directe de la concentration dans l'air se sont avérées peu satisfaisantes du fait de la mauvaise connaissance du transport turbulent.

<figure><img src="image_22.png" /><figcaption>Figure 2.2 : *Pourcentage de conidies H. maydis libérées depuis le conidiophore en fonction de la force centrifuge appliquée, normalisée par le poids de la spore. En pointillés, approximation idéale de la force normalisée par le poids et la taille de la spore* [Aylor 78]</figcaption></figure>

Pour comprendre le mécanisme de libération passive, les premiers travaux ont été consacrés à l'évaluation des forces nécessaires pour arracher les particules biotiques de leur structure d'ancrage. Il a été ainsi mis en évidence que certaines particules (*Erysiphe Graminis*) nécessitent des forces d'arrachage très faibles, alors que d'autres exigent des forces équivalentes à 2000 fois leur propre poids (*Helminthosporium maydis* figure 2.2) [Aylor 78]. Ceci révèle que les forces nécessaires pour la libération des particules biotiques couvrent un large éventail. En revanche, pour une même espèce, la force requise semble uniforme, bien qu'elle dépende des caractéristiques propres à chaque particule. En normalisant cette force par le poids et la taille des spores arrachées, [Aylor 75] a confirmé l'existence d'une force seuil qui induit la libération (courbe en pointillés sur la figure 2.2)
dépendante seulement de la force de cohésion des particules à leur support. Il conclut que l'analyse de la structure de connexion des conidies, et en général de la structure d'ancrage des particules biotiques, permettrait certainement une meilleure compréhension des mécanismes de libération dans l'atmosphère.

<figure><img src="image_23.png" /><figcaption>Figure 2.3 : Analyses en soufflerie de certains paramètres pertinents dans le mécanisme de libération des urédospores de Puccinia recondita et de Puccinia striiformis [Geagea 97]</figcaption></figure>

De nombreuses études expérimentales ont été menées, en soufflerie et centrifugeuse [Geagea 97], pour déterminer cette fois-ci les vitesses requises pour la libération de particules, associées par [Aylor 75] à une force seuil. Il en ressort que les particules ne sont pas libérées de façon continue avec l'accroissement de la vitesse de l'écoulement, mais plutôt pour une vitesse supérieure à une valeur seuil, caractéristique à chaque espèce (figure 2.3a). Pour les conidies *H. maydis*, la vitesse seuil de libération, mesurée au niveau des conidies en conditions contrôlées, est égale à $5 \text{ ms}^{-1}$ [Aylor 78]. Or, en milieu naturel, pour que la vitesse de l'écoulement atteigne cette valeur à hauteur des conidies, un vent moyen de $25 \text{ ms}^{-1}$ est nécessaire<fnref n="2" /> [Aylor 78]. Cependant, l'auteur a observé que des vents modérés, inférieurs à $25 \text{ ms}^{-1}$, suffisent à libérer ces particules, conclusion large-ment confirmée par d'autres auteurs et pour d'autres particules. Ce résultat surprenant a

<footnote n="2">2. Ces calculs issus de Schlichting (1960) sont effectués en considérant une couche limite établie</footnote>
été dénommé «paradoxe» de la libération des particules biotiques et a suscité un grand intérêt auprès de la communauté [Grace 76, Aylor 78, Jackson 99].

La figure 2.3b illustre la complexité du mécanisme de libération passive. Alors que les urédospores *Puccinia striiformis* exposées pendant 10 secondes ne sont pas libérées pour un vent inférieur à $0.7 \, ms^{-1}$, plus de 80% sont mises en suspension par une vitesse de vent identique, mais pour un temps d'exposition de 60 secondes. En revanche, pour des vitesses supérieures ($u = 2.3 \, ms^{-1}$, $u = 2.8 \, ms^{-1}$ et $u = 3.3 \, ms^{-1}$) une exposition de 10 secondes suffit à libérer le même pourcentage de particules. Ces résultats ont été interprétés en terme de stress mécanique : les particules sont libérées, suite à certain stress mécanique, qu'il soit bref et énergique ou long avec des vitesses de vent plus faibles, pourvu qu'une rupture de la croissance de la couche limite sur la feuille se produise [Geagea 97].

En conclusion, il apparaît que les bouffées turbulentes énergiques qui ont lieu dans les couverts végétaux ont un rôle majeur dans l'arrachement des particules de leur support. Lors de l'apparition de ces rafales, la couche limite au niveau des feuilles met un temps approximatif de l'ordre de $10^{-3}$ à $10^{-4}$ secondes pour s'établir<fnref n="3" /> [Aylor 78]. Ce bref délai suffit à mettre en suspension les particules. Les mêmes phénomènes interviennent dans la remise en suspension des particules. De ce fait, une approche de libération et de remise en suspension, telle que la suggère Bagnold, en fonction d'une contrainte seuil appliquée sur le vent moyen ou sur la vitesse de frottement (équation (2.3)) n'est pas appropriée pour les particules biotiques. Dans le cadre de ce travail, principalement orienté vers le transport et le dépôt, cet aspect ne sera pas pris en compte et nous nous contenterons d'imposer une densité de source.

Enfin, il est à noter aussi que la libération des particules peut également avoir lieu, suite à l'agitation des feuilles par le vent [Aylor 78] et par l'impact des gouttes de pluie qui n'induiront toutefois que de faibles distances de dispersion [Geagea 00].

### 2.1.2 Le transport turbulent

La dispersion dans l'atmosphère correspond aux mécanismes couplés du transport par l'écoulement moyen et de la diffusion turbulente. Les transferts verticaux sont principalement des processus turbulents et sont générés par les mouvements à petite échelle qui brassent l'air et permettent la dilution et le mélange des particules. Les mouvements à plus grande échelle vont en assurer le transport.

<footnote n="3">3. Ces calculs sont issus des formules de Schlichting (1960)</footnote>
Les distances de dispersion dans l'atmosphère sont très étendues, de quelques mètres à la taille d'un continent (pour s'en convaincre, il suffit de considérer la catastrophe de Chernobyl le 26 avril 1986). Ces distance sont en partie imposées par les considérations aérodynamiques de la particule, qui dépendent de sa forme, de sa taille et de sa masse volumique. Le transport des particules dans l'atmosphère est communément décomposé selon 3 modes, dont les critères sont basés sur la taille des particules [Bagnold 41, Greeley 87]:

- Le transport par suspension où les particules ont un comportement qui s'apparente à celui des gaz et ne sont pratiquement jamais en contact avec le sol. Ce mode de transport concerne les particules de petites tailles dont le diamètre $d_s$ est inférieur à $80\mu m$;

- La saltation est le mode de transport qui concerne les particules de tailles comprises entre $80\mu m$ et $2000\mu m$. Les particules mises en mouvement sont transportées par l'écoulement mais restent fréquemment en contact avec la surface au sol, où elles peuvent être soit emprisonnées soit remises en suspension;

- Enfin, la reptation est le mode de transport des grosses particules, dont la taille est supérieure à $2000\mu m$ et qui roulent sur le sol. Ces cas ne présentent aucun intérêt pour cette étude.

Le mode de transport est donc étroitement lié aux caractéristiques de la particule et dépend bien évidemment du vent. La première expression utilisée pour évaluer la distance $x$ parcourue par une particule, de vitesse de sédimentation $v_s$, libérée à la hauteur $H$ et transportée par un vent, supposé stationnaire, est [McCartney 90]:

$$x = \frac{\overline{uH}}{v_s} \qquad (2.4)$$

Bien que cette équation ne soit qu'une estimation grossière de la distance de dispersion, elle permet toutefois de mettre en évidence les paramètres importants. La vitesse de sédimentation va caractériser le dépôt des particules par effet de gravité et va définir le potentiel de dispersion de la particule. Plus la vitesse de sédimentation est petite plus le potentiel de dispersion est grand, lorsque le vent est faible. En revanche, pour des écoulements turbulents, c'est bien la turbulence qui fixe la distance de dispersion, et reste l'agent dominant du transport éolien. Enfin, la hauteur de la source $H$ influence fortement la dispersion et reflète l'importance de la position de la surface émettrice [Belot 76].
La dispersion des particules biotiques s'effectue essentiellement dans les basses couches de l'atmosphère. Elle dépend principalement de la turbulence, de la stabilité thermique de l'atmosphère et des caractéristiques de la surface terrestre. Toutefois, des échantillonnages (par avion) à haute altitude ont révélé la présence de particules de pollen dans le haut de la couche limite atmosphérique (à 2.4 kilomètres [Raynor 74]). Dans la sous-couche rugueuse, en l'occurrence dans les couverts végétaux, le transport des particules va dépendre en plus de la turbulence dans le couvert végétal où les processus de transfert sont différents de la couche de surface.

### 2.1.3 Sédimentation et impact, les mécanismes de dépôt

[Chamberlain 67] s'est intéressé aux dépôts des particules biotiques sur des surfaces naturelles rugueuses. Dans ce contexte, il a montré qu'une particule se dépose à une vitesse égale à sa vitesse de sédimentation, lorsque le vent est faible ou pour de faibles cisaillements. En revanche, au fur et à mesure que la vitesse de l'écoulement augmente, cette vitesse devient progressivement plus importante que la vitesse de sédimentation : c'est le mécanisme d'impact inertiel qui va accentuer le dépôt des particules. Les particules chutent par sédimentation gravitationnelle et par impact sur les éléments de rugosité de la surface. On définit alors la vitesse de dépôt comme étant la vitesse intégrale à laquelle elles se déposent par sédimentation et par impact.

La compréhension de ce phénomène est largement facilitée dans le cas des couverts végétaux, où le mécanisme de dépôt a été bien compris, contrairement à la libération des particules. À l'exclusion des particules submicroniques, pour lesquelles le mouvement brownien est dominant, le dépôt des particules plus lourdes est gouverné par deux processus. Les particules peuvent se déposer par sédimentation, processus dû à la gravité et chuter par impact sur les éléments du couvert végétal, tels que le tronc, la surface des feuilles. Le dépôt par impact se conçoit aisément dans la sous-couche rugueuse, où nous avons vu que les éléments végétaux jouent un rôle de sources/puits de matière : l'impact est donc un terme puits induit par la présence du couvert. Cependant toutes les particules déposées, en particulier à la surface des feuilles, ne restent pas nécessairement jusqu'à germer ou sporifier. Si la surface de la feuille est sèche et non collante, la particule peut être entraînée à nouveau par l'écoulement. [Chamberlain 67] considère ces bonds successifs de la particule comme le transport par saltation.

Ces observations confirment que l'analyse des processus de dépôt constitue un élément
supplémentaire à la compréhension de la dispersion dans l'atmosphère. Il semble que les caractéristiques de la surface jouent un rôle essentiel dans le dépôt des particules. Par conséquent, des discontinuités de surface peuvent modifier la dispersion des particules biotiques. L'objectif final de notre travail est de permettre une meilleure compréhension de l'influence de la structure d'un paysage sur le transport turbulent et le dépôt de ces particules.

## 2.2 L'aérodynamique des particules biotiques

Le potentiel de dispersion d'une particule est déterminé par sa vitesse de sédimentation dans l'air. Bien que les particules biotiques aient souvent des formes complexes, la vitesse de sédimentation de ces particules s'appuie sur l'estimation de la vitesse limite de chute, accessible en régime de Stokes [Chamberlain 67, Crane 86].

### 2.2.1 Vitesse de sédimentation des particules biotiques

La vitesse de sédimentation<fnref n="4" /> d'une particule sphérique de masse volumique $\rho_p$, supérieure à celle du fluide et de diamètre $d_p > 1\mu m$, est donnée en régime de Stokes, par :

$$v_s = \frac{\rho_p d_p^2 g}{18 \mu} \quad \text{pour } R_{ep} < 1 \qquad (2.5)$$

Cette expression de $v_s$ est largement utilisée pour les particules biotiques de diamètre caractéristique $d_p$ compris entre 1 $\mu m$ et 70 $\mu m$ et de masse volumique proche de celle de l'eau [Di-Giovanni 95]. Cependant le critère de passage du régime de Stokes au régime transitoire est plutôt fixé à $R_{ep} < 0.5$ pour des particules sphériques et $R_{ep} < 0.4$ pour des particules cylindriques.

Pour les particules plus fines, $v_s$ est sous estimée, et une correction doit être effectuée (correction de Cunningham). Pour des particules plus grosses, l'équation (2.5) tend à surestimer $v_s$, car les forces d'inertie du fluide deviennent dominantes et le coefficient de traînée donné par Stokes n'est plus valable. Son expression doit alors être corrigée. Cependant, ces corrections ne sont pas nécessaires pour le plus grand nombre de particules biotiques [McCartney 90].

Les particules présentes dans l'atmosphère ont des formes plus ou moins complexes, qui s'éloignent d'une sphère. Les particules compactes en forme de cylindrique ou d'ellipsoïde

<footnote n="4">4. Les grandeurs et les procédures de calcul sont présentées en annexe B.1</footnote>
<table><caption>Tableau 2.2 : *Facteur de forme dynamique pour des particules compactes d'après [Chamberlain 75] et [McCartney 90]*</caption><thead><tr><td>Forme des particules</td><td>Facteur de forme dynamique α [-]</td></tr></thead><tbody><tr><td>Sable</td><td>1.57</td></tr><tr><td>Ellipsoïde rapport 1 : 4</td><td>1.28</td></tr><tr><td>Cylindre rapport 1 : 1</td><td>1.06</td></tr><tr><td>Cylindre rapport 1 : 4</td><td>1.32</td></tr><tr><td>Quatre sphères en contact</td><td>1.57</td></tr></tbody></table>

ont une vitesse de sédimentation inférieure à celle d'une particule sphérique, car elles ont tendance à sédimenter de telle sorte que l'axe le plus grand (longueur du cylindre) soit à l'horizontale, ce qui augmente la force de traînée. Cette prise en compte de l'augmentation de la traînée peut se faire à l'aide du facteur de forme dynamique, défini comme le rapport de la force de traînée réelle, qui agit sur la particule non sphérique, et de la force de traînée agissant sur une sphère de même volume et à la même vitesse. Pour un écoulement en régime de Stokes, la vitesse de sédimentation d'une particule compacte non sphérique s'écrit :

$$v_s = \frac{\rho_p d_p^2 g}{\alpha 18 \mu} \qquad (2.6)$$

α est le facteur de forme dynamique qui ne dépend pas de la dimension de la particule [Chamberlain 75]. α est évalué théoriquement pour des formes relativement simples (tableau 2.2).

Certaines particules, notamment les spores fongiques, sont dispersées dans l'atmosphère de façon isolée ou en clusters. Ferrandino et Aylor [McCartney 90] ont déterminé la vitesse de sédimentation d'un cluster, composé de n spores de *Lycopodium* et de pollen d'*Ambrosia elatior* L. et précisent que la vitesse de sédimentation semble non corrélée de la forme du cluster :

$$v_{sn} = 0.98 v_s n^{0.53} \qquad (2.7)$$

Les particules dont la structure est à « plume » ou à « ailes », telles que la graine *Also-mitra macrocarpa*, posent des problèmes très complexes pour l'estimation de la vitesse de chute, car les forces de portance ainsi que les forces de Magnus ne peuvent plus être néglig-
gées. Il n'existe pas de formule pour prédire $v_s$ dont le seul recours reste la détermination expérimentale.

De façon générale, l'expression de $v_s$ est utilisée pour les particules simples, mais de nombreuses expérimentations sont menées car elles donnent des résultats plus satisfaisants. Enfin, l'utilisation de la vitesse de chute permet de prendre en compte toutes les caractéristiques de la particule au travers d'un seul paramètre mesurable.

### 2.2.2 Nombre de Stokes

Pour les écoulements chargés en particules solides, il est nécessaire d'introduire le nombre de Stokes. Le nombre de Stokes est une mesure sans dimension du temps de réponse des particules aux mouvements du fluide porteur. $S_t$ est défini comme le rapport :

$$S_t = \frac{\tau_p}{\tau_f} \qquad (2.8)$$

- $\tau_p$ est le temps de réponse ou de relaxation des particules qui sont soumises aux variations de l'écoulement. Ce temps caractéristique dépend seulement de la masse et de la mobilité de la particule. Il traduit donc les effets d'inertie de la particule <fnref n="5" />;

- $\tau_f$ est une échelle de temps caractéristique de l'écoulement du fluide porteur.

$S_t$ est donc un indicateur du comportement des particules solides vis-à-vis de celui du fluide. D'un point de vue pratique, un nombre de Stokes très inférieur à l'unité, signifie que le mouvement de la particule est très sensible aux fluctuations du fluide porteur, et s'adapte de manière quasi instantanée aux modifications de la vitesse. En revanche, lorsque le nombre de Stokes est grand devant l'unité, la particule ne suit pas toutes les fluctuations du fluide. C'est essentiellement son inertie qui gouverne son mouvement. Enfin, entre ces deux cas extrêmes, la particule ne répond qu'à certaines fluctuations du fluide porteur. Il a été démontré que pour un nombre de Stokes proche de l'unité, le mouvement de la particule est couplé aux mouvements des gros tourbillons, présents dans l'écoulement. Dans ce cas, les particules sont très efficacement dispersées et peuvent même être dispersées plus vite que les éléments fluides eux-mêmes (figure 2.4) [Crowe 93].

Différentes expressions du nombre de Stokes sont utilisées dans la littérature car le temps caractéristique de l'écoulement du fluide porteur $\tau_f$ dépend du problème physique étudié. Par exemple, $\tau_f$ est donné par l'échelle de temps des gros tourbillons $\tau_f = l'/u'$

<footnote n="5">5. Ces notions sont présentées en annexe B.2</footnote>
<figure><img src="image_24.png" /><figcaption>Figure 2.4 : Influence de nombre de Stokes sur la dispersion des particules dans des structures cohérentes [Crowe 93]</figcaption></figure>

($l'$ et $u'$ sont respectivement une échelle de longueur et une vitesse caractéristiques de ces gros tourbillons) [Crowe 93]. Une autre échelle de temps appropriée est l'échelle de temps de la turbulence ($C_L \frac{k}{\epsilon}$, $C_L$ est une constante à déterminer expérimentalement).

Le nombre de Stokes est également un indicateur de l'interaction particules-turbulence. Des études expérimentales et numériques s'intéressent à et s'interrogent sur la modification de la turbulence de l'écoulement due à la présence des particules solides et aux interactions particules-fluide [Gore 89, Elghobashi 94]. Il semble que pour des nombres de Stokes inférieurs à l'unité, la phase particulaire, dont le diamètre est inférieur à la micro-échelle de Kolmogorov, n'affecte que très localement le fluide et ces effets sont négligeables sur la turbulence. En revanche, pour des nombres de Stokes plus grands, la présence des particules peut amplifier la turbulence du fluide porteur. Dans le cadre de ce travail, l'interaction des particules sur la turbulence de l'écoulement n'est pas traitée.

## 2.3 Revue sur les principaux modèles de dispersion

Il existe de nombreux modèles pour la description mathématique de la dispersion appliquée aux particules biotiques. Il se distingue les modèles empiriques, qui se basent essentiellement sur des observations expérimentales, des modèles physiques, qui intègrent de manière plus ou moins complexes des informations telles que le mouvement turbulent du vent, la diffusion turbulente, la structure du paysage, le dépôt des particules. Dans cette partie, nous présenterons de manière succincte, les modèles physiques couramment utilisés pour la dispersion atmosphérique, en excluant l'approche eulérienne, car nous la
détaillerons dans la section suivante.

### 2.3.1 Modèles Gaussiens

Les modèles gaussiens (modèles de panache) sont des modèles eulériens statistiques, couramment utilisés pour la prédiction de la dispersion des polluants atmosphériques. De plus en plus, ils sont adaptés et appliqués à la dispersion des particules biotiques, en particulier les spores [McCartney 85] et les graines [Klein 00]. Ces modèles supposent que la distribution moyenne dans l'espace de la concentration émise par une source ponctuelle est décrite par une courbe gaussienne.

La formulation gaussienne découle par analogie directe de la loi de Fick, qui donne une relation entre la répartition d'une grandeur extensive dans l'espace et le flux (modélisation phénoménologique). Soit $C$, une concentration locale quelconque définie par la masse du soluté par unité de volume du fluide [kg m⁻³]. En établissant le bilan de masse de $C$ sous forme locale, l'équation instantanée de conservation de masse n'est autre que son équation de transport, appelée équation de convection-diffusion<fnref n="6" />. Dans un milieu et un écoulement homogène, dans le cas où seul le transport par diffusion est considéré suivant une seule direction, par exemple $x$, l'équation de diffusion admet une solution [Graf 96]:

$$C(x,t) = \frac{M_1}{\sqrt{4\pi k_c t}} \exp \left( -\frac{x^2}{4k_c t} \right) \quad (2.9)$$

Cette équation exprime que l'évolution de la concentration dans le temps et dans l'espace est représentée par une courbe gaussienne, qui décrit l'étalement de la matière $M_0$ [kg] initialement introduite de façon ponctuelle. $M_1 = M_0/S$ est la source surfacique [kg m⁻²]. En posant :

$$\sigma_x^2(t) = 2k_c t \quad (2.10)$$

$\sigma_x^2$ est la variance de la distribution de la répartition de la matière selon $x$ et $\sigma_x$ correspond à une échelle d'étalement de la matière. La variance et l'écart type de cette distribution sont donc des indicateurs de dispersion. En introduisant ces grandeurs, l'équation (2.9) devient :

$$C(x,t) = \frac{M_1}{\sqrt{2\pi}\sigma_x} \exp \left( -\frac{x^2}{2\sigma_x^2} \right) \quad (2.11)$$

De façon similaire, la concentration moyenne locale $\bar{C}(x, y, z)$ [particule m⁻³], à l'aval d'une source $Q$ continue placée à la hauteur $H$, est représentée par une courbe gaussienne

<footnote n="6">6. Équation (2.24) à la section 2.5</footnote>
de paramètres $\sigma_y(x)$ et $\sigma_z(x)$:

$$ \bar{C}(x, y, z) = QT = \frac{Q}{\bar{u}} \frac{1}{\sqrt{2\pi}\sigma_y(x)} \exp\left(-\frac{y^2}{2\sigma_y^2(x)}\right) \\ \frac{1}{\sqrt{2\pi}\sigma_z(x)} \left[ \exp\left(-\frac{(H-z)^2}{2\sigma_z^2(x)}\right) + \exp\left(-\frac{(H+z)^2}{2\sigma_z^2(x)}\right) \right] \quad (2.12) $$

$\bar{u}$ est la vitesse moyenne à la hauteur $H$ dans la direction de l'écoulement, considérée comme constante. $\bar{u}^{-1}$ peut être interprété comme un facteur de dilution puisque le nuage se déplace à la vitesse $\bar{u}$. Une source imaginaire est placée en $-H$ pour une paroi imperméable et réfléchissante. La différence fondamentale entre cette description du champ moyen de concentration et l'équation (2.11) est dans la signification physique de $\sigma_y(x)$ et $\sigma_z(x)$, qui expriment ici la diffusion turbulente du nuage suivant les directions $y$ et $z$. Ils dépendent donc des caractéristiques de l'écoulement, de la turbulence, de la stabilité de l'atmosphère et de la rugosité de la surface. Aussi, la difficulté majeure de ces modèles réside dans la formulation des expressions des écarts types du panache en fonction des conditions locales de l'écoulement, notamment dans les couverts végétaux.

Le formalisme de l'équation (2.12) suppose une turbulence homogène et stationnaire, puisque la vitesse moyenne et les écarts types du panache sont pris constants pour une configuration donnée. De surcroît, elle tend à surestimer le nuage particulaire car le dépôt des particules est négligé. Cette prise en compte du dépôt peut être partiellement envisagée en réduisant le terme source $Q$. Ce terme d'épuisement ne modifie pas la forme du profil de concentration, qui demeure gaussien, et considère que la force de la source $Q$ tend à diminuer à mesure que le panache s'éloigne de la source.

### Prise en compte du dépôt dans la couche de surface

Le dépôt des particules est influencé par les caractéristiques du sol et principalement par la rugosité. [Chamberlain 67] a montré que le flux de particules déposées $F_d$ est proportionnel à la concentration $\bar{C}$ prise à une hauteur $z_p$ près du sol et à une vitesse dite vitesse de dépôt $v_d$:

$$ F_d(x, y) = v_d \bar{C}(x, y, z_p) \quad (2.13) $$

Horst (1977) a utilisé la notion de vitesse de dépôt pour exprimer l'épuisement de la force $Q$ [McCartney 85].
## Prise en compte du dépôt dans le couvert végétal

Dans les couverts végétaux, le dépôt des particules est influencé par les éléments végétaux. Un terme puits $S_p$ [particule $m^{-3}s^{-1}$] est introduit et correspond au dépôt des particules par unité de volume et de temps (sa valeur est donc négative) :

$$S_p = r\bar{C}(x, y, z) \qquad (2.14)$$

où $r$ [$s^{-1}$] est la cadence de dépôt dans le couvert, qui dépend de la structure du couvert, de la vitesse dans le couvert ainsi que de la taille des particules. $r$ traduit les mécanismes de dépôt par sédimentation et par impact. En introduisant la cadence de dépôt dans $Q$, les modèles gaussiens peuvent modéliser la dispersion dans le cas de couverts homogènes [McCartney 85].

Les modèles gaussiens utilisent des hypothèses restrictives, telles que des vents uniformes et des cadences de dépôts constantes. Dans des situations complexes, ces modèles surestiment les champs de concentration et les hypothèses d'homogénéités pèsent lourde-ment sur l'exactitude des résultats. En revanche, ils sont efficaces pour la prédiction de la dispersion dans une couche de surface homogène, en utilisant une vitesse de dépôt ap-propriée. Ils présentent alors l'avantage d'être simples d'utilisation et peuvent s'appuyer sur des informations météorologiques, ce qui explique leur utilisation courante pour des études de contrôle de la qualité de l'air.

### 2.3.2 Modèles Lagrangiens

Il ne s'agit plus ici de considérer un nuage de particules dans son ensemble mais chaque particule est prise individuellement. Cette approche consiste donc à suivre, au sens lagran-gien, une particule individuelle dans son mouvement et à reconstituer sa trajectoire en tenant compte de l'interaction avec l'écoulement du fluide. Enfin, le traitement statistique d'un grand nombre de trajectoires permet d'obtenir des informations sur la diffusion, la concentration et la distance de dispersion en valeur moyenne.

Le point de départ des modèles lagrangiens est l'équation du mouvement de la particule solide. Bien que cette équation contienne de nombreux termes [Laslandes 98], elle peut néanmoins se simplifier dans le cas de particules lourdes $\rho_p \gg \rho$, pour un écoulement diphasique dilué. L'équation du mouvement se réduit alors à [Pozorski 98]:

$$\frac{d\mathbf{v}_p}{dt} = \frac{\mathbf{v}[\mathbf{x}_p(t), t] - \mathbf{v}_p}{\tau_p} + \mathbf{g} \qquad (2.15)$$