# L'analyse dans le sol à Météo-France. Partie 1 :
## Évaluation et perspectives à l'échelle locale

Jean-François Mahfouf

► To cite this version:

Jean-Francois Mahfouf. L'analyse dans le sol à Météo-France. Partie 1 : Évaluation et perspectives à l'échelle locale. Note de travail du GMME N°84, Centre National de Recherches Météorologiques; Groupe de Météorologie à Moyenne Echelle. 2007. hal-04390802

HAL Id: hal-04390802

https://hal.science/hal-04390802v1

Submitted on 16 Jan 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
NOTE
DE CENTRE

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

N° 84
GROUPE DE METEOROLOGIE A MOYENNE ECHELLE
L'ANALYSE DANS LE SOL A METEO-FRANCE
PARTIE 1 : EVALUATION ET PERSPECTIVES
A L'ECHELLE LOCALE

JUIN 2007

Jean-François MAHFOUF

CENTRE NATIONAL
DE RECHERCHES MÉTÉOROLOGIQUES

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

METEO FRANCE
<figure><img src="image_4.png" /></figure>

METEO-FRANCE

METEO-FRANCE

CENTRE NATIONAL
DE RECHERCHES METEOROLOGIQUES

NOTE DE TRAVAIL
DU GROUPE DE METEOROLOGIE A MOYENNE ECHELLE

N°84

L'ANALYSE DANS LE SOL A METEO-FRANCE
PARTIE 1 : EVALUATION ET PERSPECTIVES
A L'ECHELLE LOCALE

par

Jean-François MAHFOUF

Cette note n'engage pas la responsabilité de Météo-France

JUIN 2007
## RESUME

Je présente dans cette note une évaluation à l'échelle locale et en mode forcé de l'analyse dans le sol utilisée opérationnellement à Météo-France pour le modèle ARPEGE. Cette analyse, basée sur une approche de type "interpolation optimale", utilise les observations de température et humidité à 2 mètres. Elle est comparée à une méthode très similaire utilisée opérationnellement au Centre Européen (CEPMMT), ainsi qu'à une nouvelle approche de type "2D-Var simplifié" qui pourrait la remplacer dans un futur proche. J'examine également une analyse de type "filtre de Kalman d'ensemble" (EnKF) actuellement très en vogue dans la communauté hydrologique américaine. Dans une deuxième partie, j'assimile des observations d'eau du sol superficiel avec une répétitivité de trois jours par 3 méthodes (2D-Var, Filtre de Kalman étendu, EnKF). Je montre qu'il n'est pas nécessaire d'utiliser une fenêtre d'assimilation de 10 jours dans un système d'analyse cyclée. L'impact de ces observations n'est significatif que si leur erreur est comparable à celle du champ d'essai. Finalement, je regarde la possibilité de combiner ces observations avec les données à 2 mètres disponibles plus fréquemment. L'ensemble de ces tests utilise des observations simulées (expériences de type "jumeaux identiques"). Les résultats obtenues sont extrêmement encourageants et plaident en faveur du développement d'une assimilation de surface "multi-observations" dans le système externalisé SURFEX au moyen d'une méthode de type "2D-Var simplifié".
# Sommaire

1 Introduction 2

2 Définition des expériences 5

3 Les méthodes d'analyse 7

3.1 Interpolation optimale 7

a Formulation Météo-France 7

b Formulation CEPMMT 8

3.2 Le 2D-Var simplifié 8

3.3 Le filtre de Kalman d'ensemble 11

4 Présentation des résultats 14

4.1 Assimilation de $T_{2m}$ et $HU_{2m}$ 14

a Les méthodes de type OI 14

b Le filtre de Kalman d'ensemble 21

4.2 Assimilation de $w_g$ 29

4.3 Assimilation combinée de $T_{2m}$, $HU_{2m}$ et $w_g$ 36

5 Discussion et conclusions 43
# Chapitre 1

## Introduction

L'analyse des variables dans le sol (contenu en eau, température) pour la prévision numérique (PN) du temps a fait l'objet de nombreuses études depuis une quinzaine d'anneées. La connaissance de l'état hydrique des sols est aussi utile pour d'autres applications comme l'hydrologie ou la prévision saisonnière. Les problèmes que posent ces variables sont leur échelle de temps longue par rapport aux phénomènes atmosphériques et l'absence de mesures directes. D'un point de vue opérationnel, l'analyse se base essentiellement sur l'utilisation des mesures du réseau d'observation de surface pour la température et l'humidité relative afin de corriger les variables dans le sol (Météo-France, CEPMMT, DWD, SMC, HIRLAM). La technique employée est une approche de type "interpolation optimale" avec des coefficients analytiques calibrés (Bouttier et al., 1993) ou estimés dynamiquement (Hess 2001). Dans un avenir proche, des observations d'origine satellitaire informatives sur l'eau du sol superficiel vont pouvoir être disponibles en temps quasi-réel pour des études de faisabilité et éventuellement des applications opérationnelles (ASCAT sur MetOP, MIRAS sur SMOS).

L'objectif de ce travail est de commencer à définir les bases du prochain système d'analyse des variables dans le sol pour la PN à Météo-France. Deux limitations existent actuellement dans la pratique opérationnelle à Météo-France: le système d'analyse dans le sol n'a été développé que pour le modèle global ARPEGE. Les autres modèles (ALADIN, ALARO, AROME) ne font qu'interpoler spatialement cette analyse à leur résolution. L'interpolation optimale (OI) utilise des coefficients empiriques calibrés qui ne permettent pas d'assimiler d'autres observations que les variables à 2 mètres. C'est ce second problème qui est le plus gênant pour l'évolution de l'analyse actuelle, mais nous verrons que la solution proposée conduit aussi à résoudre le premier. Il a été résolu pour l'atmosphère en développant un formalisme variationnel du système d'analyse (3D-Var et 4D-Var). Divers auteurs ont montré que le formalisme variationnel peut aussi s'appliquer à l'analyse des variables du sol (Mahfouf, 1991; Bouyssel et al., 1999; Rhodin et al., 1999), mais cette approche a été considérée trop coûteuse pour une utilisation opérationnelle tout en posant le problème du couplage avec l'analyse atmosphérique. L'approche de type "2D-Var simplifié"(que l'on peut aussi nommer "OI dynamique") initialement proposée par Hess (2001) et testée dans le modèle ALADIN par Balsamo et al. (2004) est séduisante dans la mesure où elle ne requiert qu'une intégration supplémentaire du modèle de PN pour une estimation du Jacobien de l'opérateur d'observation (en traitant les colonnes dans le sol de manière indépendante). Comme le formalisme variationnel, elle permet d'assimiler des observations diverses de manière simultanée sur une fenêtre temporelle. L'étude de Seuffert et al. (2004) a en effet montré l'intérêt de cette approche pour l'assimilation combinée de températures de brillance microondes et d'observations conventionnelles au niveau anémométrique. Toutefois, à Météo-France le
gain n'a pas été jugé suffisant par rapport au coût de la méthode pour qu'elle devienne opérationnelle.

Parallèllement, toujours à Météo-France, des études de faisabilité sur l'utilité d'observations du contenu en eau superficiel pour retrouver une information sur le contenu en eau de la zone racinaire (variable influençant le plus fortement les flux d'eau à l'interface) étaient menées par Calvet et al. (1998) et Calvet et Noilhan (2000). Contrairement aux études sur l'assimilation des variables à deux mètres qui utilisent un champ d'ébauche issu d'une prévision atmosphérique, ces travaux ont été menés en mode "forcé" (forçage atmosphérique prescrit à partir d'observations ou d'analyses). Ce cadre de travail est similaire à celui de la chaîne d'analyse hydrologique SIM (Habets et al. 2007) ou du "North American Land Data Assimilation System" (Mitchell et al., 2004). Les conclusions de ces travaux montrent qu'il est possible, grâce à une série temporelle d'eau du sol superficiel tous les 3 jours sur une période d'environ 10 jours, de restituer une valeur réaliste de l'eau du sol profond. La méthode variationnelle proposée ne contenait pas de terme d'ébauche. Muñoz-Sabater et al. (2007) semblent confirmer la viabilité de cet échantillonnage avec des techniques d'assimilation moins rudimentaires (EKF, 2DVAR, EnKF). Même si les forçages ont été considérés comme parfaitement connus, ils pourraient être entachés d'erreurs sans empêcher la technique de fonctionner (voir les résultats présentés dans cette note). Cette approche en mode "forcé" a aussi été utilisée aux Etats-Unis permettant à plusieurs auteurs de montrer l'intérêt d'un filtre de Kalman d'ensemble pour assimiler des températures de brillance en bande L (Crow and Wood, 2003; Reichle et al., 2002).

Comment concilier ces deux approches ? D'une part l'approche qui consiste à forcer les variables atmosphériques près de la surface ne permet pas de les assimiler (puisqu'imposées). D'autre part la nécessité d'avoir un cycle d'assimilation long pour restituer l'eau du sol profond semble incompatible avec l'analyse atmosphérique et les besoins de la PN qui requièrent une estimation de l'eau du sol plus fréquemment qu'une fois tous les 10 jours.

Balsamo et al. (2007) ont proposé une approche en mode "forcé" mais où les paramètres de la couche limite de surface peuvent être prévus en plaçant le niveau de forçage de la valeur classique de 2 m vers un niveau plus élevé (environ 50 m). Leur étude (qui va être ici confirmée) montre que pour les variables à 2 mètres, les Jacobiens estimés en mode forcé se comparent très favorablement à ceux obtenus en perturbant les conditions initiales du modèle complet atmosphère/surface. De plus, les coefficients optimaux dynamiques sont proches des formulations analytiques utilisées opérationnellement à Météo-France et au Service Méteorologique du Canada (SMC) pour ISBA. On peut donc envisager une approche de type "2D-Var simplifié" qui ne soit pas trop coûteuse informatiquement, puisque seule la surface forcée va être utilisée pour calculer des Jacobiens par une méthode de perturbations en différences finies. Elle va permettre d'assimiler à la fois des observations de surface (comme une température de brillance en bande L ou un contenu en eau superficiel) mais aussi les observations conventionnelles à 2 mètres ($T_{2m}$, $HU_{2m}$). Le CNRM ayant développé un module de surface externalisé SURFEX utilisable à terme par l'ensemble des modèles atmosphériques, en y intégrant le module d'assimilation il sera possible d'avoir un système d'analyse spécifique à chaque modèle.

L'objectif de cette note est double. Dans un premier temps, je souhaite confirmer les résultats de Balsamo et al. (2007) à l'échelle locale en comparant les méthodes d'analyse dans le sol opérationnelles à Météo-France (Giard et Bazile, 2000) et au CEPMMT (Douville et al., 2000) avec les futures approches qu'on voudrait leur substituer (EKF, 2D-Var). Je ferai une petite incursion du côté du filtre de Kalman d'ensemble à cause de sa popularité aux Etats-Unis, mais sans en approfondir les résultats. Dans un second temps, je montrerai la possibilité assimiler des observations du contenu en eau superficiel. Cette étude sera comparée à celle déjà publiée avec des données réelles par Muñoz-
Sabater et al. (2007). Finalement je présenterai une première méthodologie permettant de combiner les observations atmosphériques disponibles fréquemment (6 heures ou moins) avec des mesures de l'eau du sol superficiel ayant une répétitivité plus faible (3 jours). Deux questions seront abordées : est-il nécessaire d'avoir une fenêtre temporelle longue (10 jours) pour assimiler l'eau du sol superficiel ? quelle précision sur la mesure de l'eau du sol superficiel permet d'extraire de l'information sur le sol profond ? Nous verrons qu'elles n'ont pas le même sens avec les méthodes présentées ici qu'avec celle proposée par Calvet et al. (1998).

Je décris tout d'abord la réalisation des expériences dans le contexte d'observations simulées (Chapitre 2). Sont ensuite présentées les méthodes d'analyse examinées dans cette note et leur mise en œuvre pratique (Chapitre 3). Après avoir présenté les résultats (Chapitre 4), je dresse un ensemble de conclusions et de perspectives sur ce travail (Chapitre 5). En annexes, je décris la constitution du forçage atmosphérique, une comparaison du schéma ISBA propre à cette étude avec la version SURFEX, ainsi qu'une étude de sensibilité à la spécification de l'erreur d'observation pour l'eau du sol superficiel.
# Chapitre 2

## Définition des expériences

J'ai tout d'abord développé une version simplifiée du schéma ISBA en me basant sur les équations décrites dans Noilhan et Mahfouf (1996) avec comme seules différences la non-prise en compte du réservoir d'interception et la description de l'évaporation du sol nu où j'ai retenu la formulation du CEPMMT qui utilise le concept de résistance de sol. Comme je ne considère ici qu'une période estivale aux moyennes latitudes (situation la plus favorable pour un contrôle de l'évaporation par l'eau du sol), le gel dans le sol et la neige ne sont pas modélisés. Je présente dans l'Annexe 2 une comparaison de ce schéma épuré avec la version ISBA-2L disponible dans SURFEX. J'ai considéré une période de 2 mois (Juillet-Août 2002) sur un site de prairie aux Etats-Unis (correspondant à la localisation géographique de l'expérience FIFE). Les données de forçage sont issues de la réanalyse de 40 ans (ERA40) du CEPMMT (disponibles en ligne sur leur site web sous forme de fichiers au format GRIB). Ces données ne sont disponibles que toutes les 6 heures et certains paramètres du forçage ont été reconstruits. J'explique en Annexe 1 comment l'échantillonnage temporel a été augmenté à une heure et comment les paramètres non disponibles ont été recalculés. Le Tableau 1 rassemble les paramètres d'entrée du schéma ISBA caractérisant la surface (sol, végétation). Ils sont associés à une végétation de type prairie dont les valeurs tabulées par Mahfouf et al. (1995) et Giard et Bazile (2000) proviennent des validations locales sur le site de FIFE réalisées par Germain (1990), notamment l'inclusion d'une dépendance de la résistance de surface vis à vis du déficit de saturation.
<table><caption>Tableau 2.1: Paramètres d'entrée du schéma ISBA pour l'ensemble des expériences présentées en utilisant les notations de Noilhan et Mahfouf (1996)</caption><thead><tr><td>Paramètre (unité)</td><td>Valeur</td></tr></thead><tbody><tr><td>Z<sub>ref</sub> (m)</td><td>50</td></tr><tr><td>Z<sub>0</sub> (cm)</td><td>10</td></tr><tr><td>Z<sub>0H</sub> (cm)</td><td>1</td></tr><tr><td>ε</td><td>0.97</td></tr><tr><td>α</td><td>0.20</td></tr><tr><td>R<sub>smin</sub>(s/m)</td><td>40</td></tr><tr><td>LAI</td><td>1</td></tr><tr><td>veg</td><td>0.85</td></tr><tr><td>d<sub>2</sub> (m)</td><td>1</td></tr><tr><td>d<sub>1</sub> (cm)</td><td>1</td></tr><tr><td>γ (kg<sup>-1</sup>kg)</td><td>20</td></tr><tr><td>R<sub>GL</sub> (W/m<sup>2</sup>)</td><td>100</td></tr><tr><td>C<sub>v</sub> (J<sup>-1</sup>m<sup>2</sup>K)</td><td>2E-5</td></tr><tr><td>% sable</td><td>50</td></tr><tr><td>% argile</td><td>33</td></tr><tr><td>Latitude (degrés)</td><td>40 N</td></tr><tr><td>Longitude (degrés)</td><td>275 E</td></tr></tbody></table>
# Chapitre 3

## Les méthodes d'analyse

Je présente ici succinctement les quatre méthodes qui seront comparées par la suite. On cherchera à analyser les quatre variables pronostiques d'ISBA ($T_s$, $T_2$, $w_g$, $w_2$) même si par la suite c'est surtout le comportement de la variable $w_2$ (eau du sol dans la zone racinaire) qui sera détaillé.

### 3.1. Interpolation optimale

Cette méthode décrite dans Mahfouf (1991) permet de relier les corrections dans le sol aux innovations de température et humidité à 2 mètres. Les coefficients optimaux ont été calculés à partir d'intégrations de type "Monte-Carlo" avec un modèle colonne. Sa mise en œuvre opérationnelle à Météo-France est décrite dans Giard et Bazile (2000) et au CEPMMT dans Mahfouf et al. (2000). Les variables dans le sol sont corrigées toutes les 6 heures en utilisant les incréments $\Delta T_{2m}$ et $\Delta HU_{2m}$ d'une analyse des variables à 2 mètres (dans la pratique tant à Météo-France qu'au CEPMMT il s'agit interpolation optimale univariée utilisant les observations SYNOP). Les équations d'analyse s'écrivent :

$$ \Delta w_g = \alpha_1 \Delta T_{2m} + \alpha_2 \Delta HU_{2m} \qquad (3.1) $$

$$ \Delta w_2 = \beta_1 \Delta T_{2m} + \beta_2 \Delta HU_{2m} \qquad (3.2) $$

$$ \Delta T_g = \mu_1 \Delta T_{2m} \qquad (3.3) $$

$$ \Delta T_2 = \mu_2 \Delta T_{2m} \qquad (3.4) $$

#### a. Formulation Météo-France

Les coefficients analytiques $\alpha_1$, $\alpha_2$, $\beta_1$ et $\beta_2$ ont été obtenus par Giard et Bazile (2000) à partir d'une formulation initiale de Bouttier et al. (2003). Il dépendent principalement de l'heure locale et de la fraction de végétation veg. D'autres dépendances (plus faibles) existent avec le rapport $R_{smin}/LAI$ et la texture du sol. Ils ont été calibrés avec une spécification des erreurs d'observations de 10 % pour l'humidité relative et de 1 K pour la température. Lorsque ces erreurs ont des valeurs différentes un facteur correctif est appliqué (ce n'est pas le cas pour la pratique opérationnelle). Plus récemment, une dépendance vis à vis de l'angle zénithal a été introduite (Bouyssel, 2007, communication personnelle) suivant la proposition de Douville et al. (2000). Les coefficients $\beta_i$ ont été réduits d'un facteur 6 par rapport à ceux décrits dans Giard et Bazile (2000). L'origine de ce facteur vient de ce que, dans la méthode de Monte-Carlo, les perturbations initiales d'eau du sol imposées avaient une dispersion
trop importante (de l'ordre de 100 mm pour un mètre de sol) par rapport aux erreurs possibles du réservoir profond sur une période de 6 heures (cet aspect est discuté plus avant dans Douville et al. (2000)). Finalement, les coefficients $\alpha_i$ et $\beta_i$ sont mis à zéro dans des conditions où les erreurs sur les paramètres à 2 mètres ne viennent pas de l'eau du sol (accumulation de précipitation sur 6-h supérieure à 0.6 mm, température de l'air négative, module du vent supérieur à 10 m/s, neige au sol). Les coefficients pour la température sont fixés de manière arbitraire à $\mu_1 = 1$ et $\mu_2 = 1/2\pi$ en accord avec l'analyse dans le sol développée avant l'implémentation opérationnelle du schéma ISBA.

## b. Formulation CEPMMT

Le schéma de surface du CEPMMT (TESSEL) étant différent de ISBA, les coefficients optimaux ont été redérivés pour ce modèle toujours par une technique de type Monte-Carlo avec une version unidimensionnelle. Toutefois, afin d'éviter la surestimation des variances d'erreurs de l'ébauche pour l'eau du sol, on a gardé la formulation explicite des coefficients optimaux en terme d'écarts-types d'erreurs et de correlations. Ainsi les écarts-types d'erreurs d'observations apparaîssent explicitement et les écarts-types d'erreurs de l'ébauche (prévision 6-h) ont été spécifiés empiriquement à 0.01 m³/m³ (en se basant toutefois sur une méthode de type NMC<fnref n="1" /> sur le champ précipitation moins évaporation pour estimer l'ordre de grandeur). Ces statistiques ont été calculées pour une journée de ciel clair en été. Comme à Météo-France, les coefficients optimaux $\alpha_i$ et $\beta_i$ sont ensuite réduits par trois fonctions empiriques $F_1$, $F_2$ et $F_3$ dépendant de l'angle zénithal (réduction la nuit et l'hiver), de la transmittivité atmosphérique (réduction en présence de nébulosité) et de l'orographie (réduction en zones montagneuses). Finalement les critères de Météo-France sur le vent, la précipitation et la température permettent d'annuler les coefficients optimaux $\alpha_i$ et $\beta_i$ dans des conditions atmosphériques non informatives sur l'état de la surface. Les corrections pour la température se font de manière complémentaire au contenu en eau (plutôt la nuit et en hiver) en appliquant le coefficient $(1 - F_1)$ :

$$ \mu_1 = 1 \times (1 - F_1)F_3 \quad (3.5) $$

$$ \mu_2 = 0.1 \times (1 - F_1)F_3 \quad (3.6) $$

## 3.2. Le 2D-Var simplifié

Les équations précédentes peuvent s'écrire sous la forme d'une équation d'analyse plus générale. On cherche à analyser un vecteur x connaissant un état a-priori (ou ébauche) $x_b$ (fourni par une prévision à courte échéance) et sa matrice de covariances d'erreurs B ainsi que des observations y₀ avec une matrice de covariances d'erreurs R. On suppose de plus que l'on est capable d'estimer l'équivalent modèle des observations au moyen d'un opérateur $\mathcal{H}$ (soit $y = \mathcal{H}(x)$) ainsi que sa matrice Jacobienne $H$ (définie par $H_{ij} = \partial y_i / \partial x_j$). L'état analysé $x_a$, caractérisé par une variance minimale, est fourni par l'équation :

$$ x_a = x_b + K[y_0 - \mathcal{H}(x_b)] \quad (3.7) $$

où la matrice de gain K est définie par :

$$ K = BH^T(HBH^T + R)^{-1} \quad (3.8) $$

<footnote n="1">¹prévisions décalées valides à une même échéance</footnote>
La matrice de covariances d'erreurs d'analyse A est donnée par :

$$A = (I - KH)B \quad (3.9)$$

Si l'on considère le problème d'analyse des variables dans le sol, le vecteur de contrôle x est de dimension 4 et vaut $(w_g, w_2, T_s, T_2)$, le vecteur d'observation $y_0$ est de dimension 3 pour cette étude avec $(T_{2m}, HU_{2m}, w_g)$. Il faut maintenant introduire la dimension temporelle (essentielle pour le problème qui nous intéresse) qui va nous conduire à définir une approche hybride entre une assimilation variationnelle et un filtre de Kalman étendu. L'équation (3.7) se réfère à un instant $t$. On peut ensuite faire évoluer $x$ par un modèle de prévision $\mathcal{M}$ jusqu'au prochain instant où des observations sont disponibles :

$$x^{t+1} = \mathcal{M}(x^t) \quad (3.10)$$

Il faut également propager la matrice d'analyse A pour obtenir une matrice d'erreur de l'ébauche au nouveau temps $t + 1$:

$$B^{t+1} = MA^T M + Q \quad (3.11)$$

où intervient la matrice Jacobienne $M$ du modèle de prévision (et sa transposée) ainsi qu'une nouvelle matrice représentant les erreurs du modèle $Q$. Cette approche est très coûteuse pour des vecteurs $x$ de grande dimension et nécessite la spécification d'une nouvelle matrice d'erreurs. On lui préfère la linéarisation d'un problème variationnel :

$$J(\mathbf{x}) = \frac{1}{2}(\mathbf{x} - \mathbf{x}_b)^T \mathbf{B}^{-1} (\mathbf{x} - \mathbf{x}_b) + \frac{1}{2}(\mathbf{y}_o - \mathcal{H}(\mathbf{x}))^T \mathbf{R}^{-1} (\mathbf{y}_o - \mathcal{H}(\mathbf{x})) \quad (3.12)$$

qui s'écrit :

$$\mathbf{x}_a^t = \mathbf{x}_b^t + \mathbf{K}_1[y_0^{t+1} - \mathcal{H}(\mathbf{M}(\mathbf{x}_b^t))] \quad (3.13)$$

On cherche un état analysé au début d'une certaine période de temps ($t$) avec des observations disponibles à la fin de la période ($t+1$). On a donc un opérateur d'observation qui combine l'évolution par le modèle et la projection de l'état modèle dans l'espace des observations :

$$\mathcal{H}_1 = \mathcal{H}\mathcal{M} \quad (3.14)$$

La matrice de gain va s'écrire :

$$\mathbf{K}_1 = \mathbf{B}\mathbf{H}_1^T(\mathbf{H}_1\mathbf{B}\mathbf{H}_1^T + \mathbf{R})^{-1} \quad (3.15)$$

Si ensuite on propage les incréments d'analyse par le modèle linéaire tangent $M$ on constate qu'à l'heure d'observation la matrice $\mathbf{B}$ effective est $\mathbf{MBM}^T$ (rappelant l'équation (3.13) sans terme d'erreur modèle).

Ayant obtenu l'incrément $\delta\mathbf{x}^t = \mathbf{x}_a^t - \mathbf{x}_b^t$ en début de fenêtre d'assimilation, l'approche dite "2D-Var simplifié" l'ajoute à $\mathbf{x}_b^t$ et fait évoluer $\mathbf{x}_a^t$ avec le modèle non-linéaire pour obtenir l'état analysé à la fin de la fenêtre :

$$\mathbf{x}_a^{t+1} = \mathcal{M}(\mathbf{x}_b^t + \delta\mathbf{x}^t) \quad (3.16)$$

Pour une période de temps relativement courte (ou pour une variable évoluant lentement dans le temps), on peut faire l'approximation :

$$\mathbf{x}_a^{t+1} \approx \mathcal{M}(\mathbf{x}_b^t) + \delta\mathbf{x}^t \quad (3.17)$$
Soit:

$$x_a^{t+1} = x_b^{t+1} + K_1[y_0^{t+1} - \mathcal{H}(M(x_b^t))] \quad (3.18)$$

Cette approche sera par la suite appelée “filtre de Kalman étendu simplifié”. Cette formulation permet d’avoir une équivalence stricte avec les formulations de l’interpolation optimale présentées précédemment qui ne combinent pas plusieurs instants. Les coefficients de la matrice de gain $K_1$ peuvent ainsi être considérés comme des coefficients optimaux dynamiques qui vont dépendre de la matrice Jacobienne :

$$H_1 = \frac{\partial y_0^{t+1}}{\partial x^t} \quad (3.19)$$

Le fait d'avoir inclu le modèle dans ce Jacobien permet de trouver des sensibilités qui n'existent pas lorsqu'on considère l'état du modèle au même instant que les observations, par exemple $\partial T_{2m}/\partial w_2$ est nul à un instant donné. On verra par la suite que cette méthode permet effectivement d'extraire de l'information sur l'eau du sol profond à partir d'observations sur le sol superficiel. Nous reviendrons sur cet aspect quand il s'agira de combiner plusieurs types d'observations disponibles à des fréquences temporelles différentes.

L'approche précédente peut se généraliser à des observations disponibles à plusieurs instants sur une fenêtre temporelle donnée:

$$x_a^t = x_b^t + \sum_i K_i [y_0^i - \mathcal{H}(M_i(x_b^t))] \quad (3.20)$$

$$\mathcal{H}_i = \mathcal{H} M_i \quad H_i = H M_i \quad (3.21)$$

Les matrices de gain $K_i$ valent :

$$K_i = B H_i^T (H_i B H_i^T + R_i)^{-1} \quad (3.22)$$

où $M_i$ permet l'évolution de l'état $x$ depuis le début de la fenêtre d'assimilation jusqu'à l'instant $i$ où le vecteur d'observations $y_0^i$ est disponible. Cette approche ressemble à un filtre de Kalman étendu mais dans lequel on ne cycle pas la matrice $B$ d'une fenêtre d'assimilation à la suivante et pour lequel l'erreur de modèle est négligée. Elle est équivalente à l'approche variationnelle dans la limite de linéarité de l'opérateur d'observation $\mathcal{H}_i$.

Dans tout ce qui vient d'être décrit l'aspect spatial n'intervient pas. On considère en effet que les observations $y_0$ sont disponibles au même endroit que le vecteur $x$ que l'on cherche à analyser, et qu'elles sont les seules à influencer cet état. On a donc un ensemble de problèmes unidimensionnels indépendants rendant la résolution numérique de l'équation d'analyse particulièrement facile. D'un point de vue pratique, on peut donc calculer la matrice Jacobienne $H$ en différences finies par un faible nombre d'intégrations perturbées, et son stockage s'effectue dans un tableau de faible dimension (ce qui ne serait pas le cas en supposant que de nombreuses observations réparties spatialement puissent affecter un point de grille donné). On peut donc écrire :

$$H_{ij} = \frac{y_i(\mathbf{x} + \delta x_j) - y_i(\mathbf{x})}{\delta x_j} \quad (3.23)$$

La valeur de $\delta x_j$ doit être assez petite pour approximer correctement la dérivée mais pas trop cependant pour éviter les erreurs d'arrondis numériques. En mode forcé une plage assez large est permise ce qui n'est pas le cas lorsque l'on perturbe la surface couplée avec un modèle atmosphérique (voir
Balsamo et al. (2004) pour une discussion relative à ce point). Les contenus en eau volumique dans le sol ont été normalisés en considérant le SWI (soil wetness index) défini comme :

$$SWI = \frac{w - w_{wilt}}{w_{fc} - w_{wilt}} \qquad (3.24)$$

où $w_{fc}$ est le contenu en eau à la capacité au champ et $w_{wilt}$ le contenu en eau au point de flétrissement. Ainsi l'analyse de l'eau du sol, et en particulier les erreurs de l'ébauche et la taille des perturbations, est définie pour cette variable. Dans la mesure où il s'agit d'une transformation linéaire elle est sans influence sur le comportement de l'analyse. On a choisi $\delta T = 0.001$ et $\delta SWI = 0.0001$ pour estimer les Jacobiens de l'opérateur d'observation. Nous verrons dans des tests de sensibilité sur la taille de ces perturbations que des effets non-linéaires peuvent rendre l'estimation des Jacobiens sensible à ce choix. J'explique aussi par la suite la raison de ce choix.

La spatialisation des observations intervient dans l'analyse dans le sol indirectement dans la mesure où, dans la pratique opérationnelle, les vecteurs d'innovations $[y_0 - H(x)]$ utilisés dans l'équation d'analyse sont en fait des incrément d'analyse $[x_a - x_b]$ issus d'une interpolation op- timale des températures et humidités à 2 mètres. On peut donc écrire l'analyse des variables dans le sol en deux étapes :

* Analyse des paramètres à deux mètres :

$$z_a = z_b + K_z [y_0 - \hat{z}_b] \qquad (3.25)$$

où $\hat{z}_b$ est l'équivalent modèle au point d'observation (interpolation spatiale)

* Analyse des variables dans le sol :

$$x_a = x_b + K [z_a - z_b] \qquad (3.26)$$

soit :

$$x_a = x_b + KK_z [y_0 - \hat{z}_b] \qquad (3.27)$$

Cette dernière équation exprime les incrément dans le sol comme une combinaison linéaire des inno- vations pour les paramètres à deux mètres, par contre l'optimalité de la matrice de gain ainsi formée $KK_z$ n'est pas assurée. Des tests seront effectués dans une étude ultérieure pour comparer cette ap- proche en deux étapes avec une analyse directe les variables dans le sol à partir des observations. En effet, dans un avenir proche nous disposerons de données satellitaires irrégulièrement réparties dans le temps et dans l'espace. Va donc se poser rapidement la question de savoir si elle peuvent être as- similées directement ou si elles doivent être d'abord interpolées spatialement avant d'être présentées à l'analyse dans le sol.

## 3.3. Le filtre de Kalman d'ensemble

On résout la même équation d'analyse que pour le filtre de Kalman étendu mais les matrices $BH^T$ et $HBH^T$ sont approximées à partir d'un ensemble de prévisions $(x_f)_i$ avec $i$ entre 1 et $N$ (Evensen, 2003). La valeur de $N$ est en général trop petite pour échantillonner complètement la distribution des erreurs. On écrit :

$$BH^T \approx (x_f - \overline{x_f})(H(x_f) - \overline{H(x_f)})^T \qquad (3.28)$$
$$ \mathbf{HBH}^T \approx \overline{(\mathcal{H}(\mathbf{x}_f) - \overline{\mathcal{H}(\mathbf{x}_f)})(\mathcal{H}(\mathbf{x}_f) - \overline{\mathcal{H}(\mathbf{x}_f)})^T} \quad (3.29) $$

avec les opérateurs de moyenne suivants :

$$ \bar{\mathbf{x}} = \frac{1}{N} \sum_{i=1}^{N} \mathbf{x}_i \quad (3.30) $$

$$ \overline{xy^T} = \frac{1}{N-1} \sum_{i=1}^{N} \mathbf{x}_i \mathbf{y}_i^T \quad (3.31) $$

Un avantage est donc de n'avoir jamais à expliciter les matrices **B**, **H** et **H**^T nécessaires au calcul de la matrice de gain **K**. Les observations peuvent être traitées par lots évitant ainsi l'inversion de trop grosses matrices, puisque HBH^T a la dimension du nombre d'observations. La dispersion de l'ensemble est obtenue en perturbant les observations de manière aléatoire pour chaque membre suivant la matrice d'erreurs **R** (il existe des variantes qui n'ont pas besoin de perturber les observations). L'utilisation d'un échantillon de taille réduite pose divers problèmes dont le plus gênant est la sous-estimation de la matrice d'erreur d'analyse qui conduit rapidement à un effondrement du filtre (les analyses deviennent de plus en plus précises conduisant à un poids de plus en plus faible des observations au fur et à mesure du cyclage). On remédie à ce problème en augmentant empiriquement la dispersion de l'ensemble après chaque analyse et en introduisant un terme d'erreur de modèle dans les équations pronostiques.

D'un point de vue pratique, j'ai mis en place un filtre dans lequel chaque analyse $\mathbf{x}_a^i$ est obtenue par :

$$ \mathbf{x}_a^i = \mathbf{x}_b^i + \mathbf{K}[\mathbf{y}_0 - \mathcal{H}(\mathbf{x}_b^i) + \varepsilon_i] \quad (3.32) $$

où $\varepsilon$ est un bruit aléatoire gaussien de moyenne nulle et de matrice de covariances **R**. La matrice de gain est calculée en utilisant les approximés de **BH**^T et **HBH**^T fournis par les Equations (3.28) et (3.29). Les équations pronostiques pour l'eau du sol $w$ sont modifiées par ajout d'un terme d'erreur de modèle $\varphi$ défini par :

$$ \varphi^{t+1} = \nu\varphi^t + \varepsilon_w \sqrt{1 - \nu^2} \quad (3.33) $$

$$ w^{t+1} = \mathcal{F}(w^t) + \varphi^{t+1} \Delta t \quad (3.34) $$

avec $\nu = 1/(1 + \Delta t/\tau)$. Il s'agit d'un modèle autorégressif d'ordre un avec une constante de corrélation temporelle $\tau$ fixée à 3 jours et une amplitude pour le bruit gaussien $\varepsilon_w$ fixée à 0.001 m³/m³/jour.

Pour éviter une sous estimation de la matrice **B**, après chaque cycle d'analyse on applique un facteur d'inflation $\alpha$ conduisant à redéfinir les membres analysés comme :

$$ \mathbf{x}_a = \bar{\mathbf{x}}_a + \alpha(\mathbf{x}_a - \bar{\mathbf{x}}_a) \quad (3.35) $$

Dans la pratique une valeur de 1.015 a été retenue pour un cycle de 6 heures et de (1.01)¹² pour un cycle de 3 jours.

Un ensemble de 100 membres a été considéré. Pour des valeurs inférieures, le résultat final est très dépendant du nombre de membres, par contre, pour des valeurs supérieures, l'analyse se comporte de la même façon. En augmentant le nombre de membres à des valeurs très élevées (1000) il est toujours nécessaire d'avoir un terme d'erreur de modèle et une inflation des covariances d'analyse, ce qui semble contradictoire avec la théorie qui dit qu'une partie de cette sous-estimation des covariances vient
de la taille réduite de l'échantillon. De même, la technique de Houtekamer et Mitchell (1998) qui consiste à séparer l'ensemble en deux pour calculer la matrice de gain d'une part et effectuer l'analyse d'autre part, n'a pas donné de résultats satisfaisants (i.e. meilleurs qu'avec un seul ensemble). Finallement, l'introduction de perturbations aléatoires (bruit blanc) sur les termes de forçage comme proposé par Reichle et al. (2002) n'a pas non plus amélioré le comportement du filtre. Cette approche n'a pas été explorée dans le détail étant donné le nombre important de spécifications arbitraires des variances pour chacun des paramètres du forçage.

Pour l'ensemble des techniques nécessitant l'inversion de la matrice $(HBH^T + R)$ j'ai utilisé une décomposition de type Cholesky pour une matrice symétrique décrite par Press et al. (1992).
# Chapitre 4

# Présentation des résultats

Le Tableau 4.1 présente les paramètres d'entrée utilisés dans les expériences d'assimilation qui sont décrites brièvement dans le Tableau 4.2.

## 4.1. Assimilation de $T_{2m}$ et $HU_{2m}$

J'ai travaillé ici dans le cadre d'observations simulées permettant de définir une vérité vers laquelle on souhaite converger. J'ai considéré comme référence (notée REF) un sol initialement saturé en eau et intégré ISBA durant un mois. Les températures et humidités relatives à 2m toutes les six heures de cette intégration sont utilisées comme "pseudo-observations" pour les expériences d'assimilation. Elles sont diagnostiquées au moyen d'une interpolation verticale dans la couche limite de surface (entre le sol et le niveau situé à 50 m) suivant la formulation de Geleyn (1988) modifiée pour distinguer les rugosités dynamique et thermique (Mascart et al., 1995). Une seconde intégration (dite "open loop" et notée OL) consiste à initialiser le sol au point de flétrissement et à laisser évoluer le schéma ISBA sans faire de corrections aux variables dans le sol. Dans les expériences d'assimilation, (notées ASSIM sur les figures) on va chercher à se rapprocher de la simulation de référence (REF) et à s'écarter du contrôle (OL). Dans les expériences d'assimilation et dans la simulation OL, le flux de précipitation du forçage de référence est divisé par deux afin d'évaluer la capacité des analyses à réduire un bias du modèle. Les températures dans le sol valent initialement 295 K pour toutes les intégrations qui commencent le 01/07/2002 a 00 UTC (correspondant à 18 h locale la veille).

### a. Les méthodes de type OI

La Figure 4.1 montre l'évolution temporelle des deux réservoirs $w_g$ et $w_2$ durant le mois de Juillet 2002 pour les trois intégrations : référence (REF), assimilation (ASSIM) et contrôle (OL). L'assimilation est ici l'interpolation optimale Météo-France (notée OLMF). La référence met en évidence l'assèchement progressif du sol depuis la saturation vers une valeur proche de la capacité au champ. Le sol superficiel se réhumidifie plus rapidement que $w_2$ lors des épisodes pluvieux visibles sur la Figure 5.1 du forçage. Après le jour 15, le sol profond est en équilibre entre les apports pluvieux et les pertes par évapotranspiration. L'évolution avec assimilation se situe bien entre les courbes du contrôle et de la référence. Les corrections au réservoir superficiel toutes les 6 heures sont rapidement perdues par transfert vers la couche profonde tant que $w_2$ est très sec (10 premiers jours). Les corrections après le jour 15 sont plus faibles car le réservoir $w_g$ s'est humidifié par les précipitations et le rappel vers le
<table><caption>Tableau 4.1: Paramètres d'entrée pour les expériences d'assimilation (en particulier les valeurs des erreurs de l'ébauche (exposant b) et des observations (exposant o))</caption><thead><tr><td>Paramètre (unité)</td><td>Valeur</td></tr></thead><tbody><tr><td>w<sub>sat</sub> (m<sup>3</sup>/m<sup>3</sup>)</td><td>0.440</td></tr><tr><td>w<sub>fc</sub> (m<sup>3</sup>/m<sup>3</sup>)</td><td>0.302</td></tr><tr><td>w<sub>wilt</sub> (m<sup>3</sup>/m<sup>3</sup>)</td><td>0.213</td></tr><tr><td>Δt (s)</td><td>900</td></tr><tr><td>σ<sup>o</sup><sub>T2m</sub> (K)</td><td>1</td></tr><tr><td>σ<sup>o</sup><sub>HU2m</sub> (%)</td><td>10</td></tr><tr><td>σ<sup>b</sup><sub>Ts</sub> (K)</td><td>1</td></tr><tr><td>σ<sup>b</sup><sub>T2</sub> (K)</td><td>1</td></tr><tr><td>σ<sup>b</sup><sub>wg</sub> (% SWI)</td><td>10</td></tr><tr><td>σ<sup>b</sup><sub>w2</sub> (% SWI)</td><td>10</td></tr></tbody></table>

<table><caption>Tableau 4.2: Récapitulatif de l'ensemble des expériences présentées et de leurs principales caractéristiques</caption><thead><tr><td>Expérience</td><td>Description</td><td>Observations</td></tr></thead><tbody><tr><td>REF</td><td>Expérience de référence (vérité)</td><td></td></tr><tr><td>OL</td><td>Expérience de contrôle (open loop)</td><td>aucune</td></tr><tr><td>OL_MF 2c</td><td>Interpolation optimale formulation Météo-France</td><td>T<sub>2m</sub>, HU<sub>2m</sub></td></tr><tr><td>OLEC 2c</td><td>Interpolation optimale formulation CEPMMT</td><td>T<sub>2m</sub>, HU<sub>2m</sub></td></tr><tr><td>2DVAR 2c</td><td>Méthode variationnelle 2D simplifiée</td><td>T<sub>2m</sub>, HU<sub>2m</sub></td></tr><tr><td>EKF 2c</td><td>Filtre de Kalman étendu simplifié</td><td>T<sub>2m</sub>, HU<sub>2m</sub></td></tr><tr><td>ENKF 2c</td><td>Filtre de Kalman d'ensemble</td><td>T<sub>2m</sub>, HU<sub>2m</sub></td></tr><tr><td>2DVAR 3d</td><td>Méthode variationnelle 2D simplifiée (3 jours)</td><td>w<sub>g</sub></td></tr><tr><td>2DVAR 9d</td><td>Méthode variationnelle 2D simplifiée (9 jours)</td><td>w<sub>g</sub></td></tr><tr><td>EKF wg</td><td>Filtre de Kalman étendu simplifié</td><td>w<sub>g</sub></td></tr><tr><td>EKF wg-2m</td><td>Filtre de Kalman étendu simplifié</td><td>w<sub>g</sub>, T<sub>2m</sub>, HU<sub>2m</sub></td></tr><tr><td>ENKF 2d</td><td>Filtre de Kalman d'ensemble</td><td>w<sub>g</sub></td></tr></tbody></table>
réservoir profond qui a été rempli par l'analyse. Le contenu $w_2$ du réservoir profond augmente graduellement au cours de la période d'environ 30 mm au cours des dix premiers jours par l'assimilation (marches d'escalier sur l'évolution de la courbe) et de 10 mm par la suite surtout par les apports pluvieux. Les corrections faibles dans la seconde moitié de Juillet s'expliquent par une diminution de la sensibilité des paramètres à 2 m lorsque l'eau du sol se rapproche de la saturation (cela sera mis en évidence par l'examen des coefficients optimaux dynamiques) et par le fait que les importantes précipitations réduisent les jours où les erreurs à 2 m sont informatives pour les corrections dans le sol.

La Figure 4.2 présente les mêmes courbes mais pour la méthode d'interpolation optimale du CEP-MMT (notée OI_EC). Le réservoir profond se comporte de manière similaire mais les corrections journalières sont plus importantes. La correction au bout de 10 jours est de 50 mm. On peut aussi noter que les corrections sont plus importantes entre les jours 1 à 5, qu'entre les jours 6 à 10. Après cette date pour les raisons évoquées précédemment, $w_2$ reste à une valeur stationnaire proche de 0.27 m$^3$/m$^3$. L'évolution de $w_g$ est assez différente de celle de la Figure 4.1 jusqu'au jour 15 essentiellement parce que les corrections apportées par l'analyse à ce réservoir sont très faibles. Par contre, après cette date ce sont les précipitations (même diminuées d'un facteur deux) qui pilotent les variations simulées.

La Figure 4.3 correspond à l'utilisation d'une assimilation 2D-Var simplifiée (notée 2DVAR) sur une fenêtre de 6 heures. L'évolution de $w_g$ est très proche de celle obtenue avec OI_EC, montrant que les fortes corrections à ce réservoir avec OI_MF ne sont pas justifiées par l'existence de sensibilités $\partial T_{2m}/\partial w_g$ et $\partial HU_{2m}/\partial w_g$. Dans la mesure où ce réservoir, de part sa faible dimension ($d_1=1$ cm), dépend surtout des apports par la pluie et de la valeur de $w_2$, ces corrections non justifiées ne sont pas dommageables. Le réservoir profond se remplit plus rapidement et de façon plus importante qu'avec les méthodes de type OI statiques. Dès le jour 10, la valeur de $w_2$ atteint 0.28 m$^3$/m$^3$. Le report des incréments du 2DVAR en fin de période plutôt qu'en début (Equations 3.16 et 3.17) pour produire un formalisme proche d'un EKF simplifié montre un comportement tout à fait similaire de la méthode d'analyse (Figure 4.4) dans la mesure où la matrice B n'est pas cyclée (évitant ainsi de spécifier la matrice Q.)

La comparaison des coefficients optimaux $\alpha_i$ et $\beta_i$ entre les trois méthodes permet de mieux comprendre les différences de comportement qui viennent d'être décrites (pour le 2DVAR et l'EKF il s'agit des éléments de la matrice de gain K). Le coefficient $\alpha_1$ montre des valeurs négatives très faibles avec le 2DVAR (sauf initialement par un effet de non-linéarité examiné plus loin) de l'ordre de 0.005 K$^{-1}$ (Figure 4.5). Les valeurs sont cependant plus fortes durant les 10 premiers jours et beaucoup plus faibles après le jour 20. Les valeurs nulles lors des épisodes pluvieux (jours 9, 18-19, 22, 29) sont des caractéristiques qui se retrouvent également sur tous les coefficients dynamiques. Par contre la modulation par le cycle diurne n'apparaît pas clairement. Les coefficients analytiques sont environ trois fois plus forts pour OI_EC et 10 fois plus forts pour OI_MF, avec une forte modulation par le cycle diurne (seules les valeurs à 12 et 18 heures locales sont non nulles). On note, en accord avec les coefficients dynamiques, des valeurs nulles lors des épisodes pluvieux.

Le coefficient $\alpha_2$ est positif avec des valeurs plus fortes pour le 2DVAR durant les 15 premiers jours que par la suite (Figure 4.6). Les valeurs non nulles avec OI_EC de l'ordre de 0.15 sont assez proches des valeurs les plus élevées des coefficients dynamiques obtenus durant la période de ciel clair (jours 10-17) après le fort épisode pluvieux du 9 Juillet. L'OI_MF fournit des coefficients 10 fois plus forts qui expliquent les fortes corrections (non justifiées physiquement) appliquées à cette variable durant l'assimilation. Pour ce coefficient on note, avec les formulations analytiques, la modulation
<figure><img src="image_5.png" /></figure>

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

Figure 4.1: Evolution de l'eau du sol superficiel (haut) et du sol profond (bas) au cours d'une as- similation basée sur l'interpolation optimale de Météo-France. Les observations de température et humidité relative à 2 mètres sont assimilées toutes les 6 heures. La courbe en tiretés correspond à la simulation de référence (REF), la courbe en pointillés à une intégration perturbée (conditions initiales et forçage de précipitations) sans assimilation (OL), et la courbe en traits pleins correspond à une assimilation cyclée toutes les 6 heures.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4.2: Identique à la Figure 4.1 mais pour une assimilation par interpolation optimale suivant la formulation du CEPMMT.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 4.3: Identique à la Figure 4.1 mais pour une assimilation basée sur un 2D-Var simplifié avec une fenêtre temporelle de 6 heures.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 4.4: Identique à la Figure 4.1 mais pour une assimilation basée sur un filtre de Kalman étendu simplifié avec une fenêtre temporelle de 6 heures.</figcaption></figure>
par le cycle diurne et les valeurs nulles en présence de pluie.

Pour les coefficients $\beta_1$ et $\beta_2$, les valeurs sont très proches entre les deux formulations analytiques expliquant le comportement similaire pour le réservoir profond entre OI_MF et OI_EC (Figures 4.7 et 4.8). Par contre l'estimation dynamique est plus élevée en début de période là où l'état hydrique est le plus éloigné de la référence et plus faibles en fin de mois (lorsque $w_2$ se rapproche de la capacité au champ). Les différences sont surtout sensibles sur l'analyse durant la première quinzaine lorsque les innovations sont les plus importantes.

En ce qui concerne les coefficients pour la température (Figures 4.9 et 4.10), on note que la valeur de $\mu_1=1$ est trop forte, la méthode dynamique fournissant une valeur proche de 0.25. Par contre la modulation par le cycle diurne utilisée au CEPMMT se justifie par l'estimation dynamique (valeurs non nulles en période nocturne). Pour le coefficient $\mu_2$, l'ordre de grandeur des approches empiriques est en accord avec la valeur dynamique de 0.1 avec un avantage pour la formulation OI_EC qui prend en compte le cycle diurne. On note avec la formulation dynamique une réduction de ce coefficient quand il pleut. Elle pourrait être introduite dans les formulations analytiques.

En conclusion, cette comparaison a montré que pour le réservoir profond $w_2$, les trois méthodes produisent des coefficients optimaux très voisins, expliquant les fortes similarités pour la convergence de cette variable en cours d'analyse. L'estimation dynamique permet de prendre en compte les non-linéarités du système pour conduire à des corrections plus importantes lorsque l'état hydrique est plus loin de la référence. Les coefficients pour $w_g$ sont clairement surestimés avec la formulation OI_MF conduisant à de fortes corrections pour ce réservoir, mais elles ne sont pas dommageables étant donnée sa faible épaisseur.

L'effet des non-linéarités est examiné sur la Figure 4.11 où les coefficients optimaux dynamiques pour $w_g$ et $w_2$ ont été recalculés en modifiant la taille de la perturbation initiale $\delta SWI$ pour l'estimation du Jacobien de l'opérateur d'observation en différences finies (Equation 3.23). En considérant $\delta SWI = -0.01$ on retrouve sur des valeurs très proches de celles présentées précédemment, à l'exception du coefficient $\alpha_1$ où la forte valeur positive initiale de la Figure 4.5 a disparue.

## b. Le filtre de Kalman d'ensemble

Les résultats pour le filtre EnKF sont présentés sur la Figure 4.12. On note un comportement similaire à la méthode EKF pour les 5 premiers jours en ce qui concerne le réservoir profond où les corrections conduisent à une valeur de l'ordre de $0.25 \text{ m}^3/\text{m}^3$. Les corrections durant les 5 jours suivants sont plus faibles. Par la suite, comme avec les autres méthodes, le sol profond est très peu corrigé. Le comportement entre les jours 5 et 10 s'explique en regardant l'évolution de $\sigma^b$ (écart-type d'erreur de l'ébauche) au cours du temps pour $w_2$ (Figure 4.13). Cette erreur décroit de la valeur initiale fixée par la matrice B de $0.09 \text{ m}^3/\text{m}^3$ (0.1 en SWI), à une valeur voisine de $0.06 \text{ m}^3/\text{m}^3$ conduisant à une confiance plus grande dans le terme d'ébauche. Durant les épisodes pluvieux et peu après, l'eau du sol profond contient peu d'information sur les variables atmosphériques et les faibles valeurs de la matrice de gain conduisent à l'augmentation de l'erreur de l'ébauche grâce au facteur d'inflation $\alpha$. Après un mois, la valeur finale est proche de la valeur initiale imposée ($0.011 \text{ m}^3/\text{m}^3$). En choisissant une valeur de $\alpha$ plus forte, des corrections plus importantes seraient apportées en début d'assimilation mais elles auraient conduit à des valeurs trop fortes de $\sigma^b_{w_2}$ en fin de période où l'analyse est moins active principalement car les innovations sont plus faibles. La valeur de $\alpha$ choisie est en accord avec les applications atmosphériques (Evensen, 2003) mais est plus faible que le choix de Muñoz-Sabater et al. (2007) avec une valeur de 1.35 avec un cycle de 3 jours (soit 1.025 toutes les 6 heures). L'examen
<figure><img src="image_13.png" /></figure>

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

<figure><img src="image_15.png" /><figcaption>Figure 4.5: Coefficient optimal $\alpha_1$ reliant l'increment de $w_g$ à l'innovation de $T_{2m}$ avec la formulation analytique de Météo-France (haut), du CEPMMT (milieu) et dynamique [éléments de la matrice de gain] (bas) en Juillet 2002 sur le site choisi (prairie US).</figcaption></figure>
<figure><img src="image_16.png" /></figure>

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

<figure><img src="image_18.png" /><figcaption>Figure 4.6: Comme la Figure 4.5 mais pour le coefficient optimal $\alpha_2$ reliant l'increment de $w_g$ à l'innovation de $HU_{2m}$</figcaption></figure>
<figure><img src="image_19.png" /></figure>

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

<figure><img src="image_21.png" /><figcaption>Figure 4.7: Comme la Figure 4.5 mais pour le coefficient optimal $\beta_1$ reliant l'increment de $w_2$ à l'innovation de $T_{2m}$</figcaption></figure>
<figure><img src="image_22.png" /></figure>

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

<figure><img src="image_24.png" /><figcaption>Figure 4.8: Comme la Figure 4.5 mais pour le coefficient optimal $\beta_2$ reliant l'increment de $w_2$ à l'innovation de $HU_{2m}$</figcaption></figure>
<figure><img src="image_25.png" /></figure>

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

<figure><img src="image_27.png" /><figcaption>Figure 4.9: Comme la Figure 4.5 mais pour le coefficient optimal $\mu_1$ reliant l’increment de $T_s$ à l’innovation de $T_{2m}$</figcaption></figure>
<figure><img src="image_28.png" /></figure>

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

<figure><img src="image_30.png" /><figcaption>Figure 4.10: Comme la Figure 4.5 mais pour le coefficient optimal $\mu_2$ reliant l'increment de $T_2$ à l'innovation de $T_{2m}$</figcaption></figure>
<figure><img src="image_31.png" /></figure>

Figure 4.11: Coefficient optimaux dynamiques $\alpha_i$ et $\beta_i$ obtenus en augmentant d'un facteur 100 la taille de la perturbation et en inversant son signe pour la période Juillet 2002 sur le site choisi (prairie US).
de l'écart-type d'erreur de l'ébauche pour $w_g$ montre que les valeurs prévues sont toujours inférieures à la valeur initiale imposée. Les valeurs sont fortes en situations de ciel clair (jours 5-9 et 11-17 notamment) pour chuter dramatiquement lors d'épisodes pluvieux. Les fluctuations importantes de ces erreurs en fonction de l'état hydrique de sol superficiel semblent un atout indéniable de leur estimation dynamique par le filtre. Toutefois, la valeur de $\sigma^b$ pour $w_g$ joue très peu dans la correction du réservoir $w_2$. C'est pourquoi avec une valeur constante, qui est une très mauvaise approximation pour cette variable, les résultats de l'analyse sont tout à fait satisfaisants. De plus, Zhou et al. (2006) ont montré que dans des conditions très sèches ou très humides les distributions d'erreurs étaient pas gaussiennes (fortes asymétries).

Pour les 4 méthodes d'assimilation, les réductions d'erreurs sur les paramètres à 2 mètres et sur le flux dévaporation sont présentées sur la Figure 4.14 comme le rapport $r$ de l'erreur quadra-tique moyenne journalière de chacune des expériences d'assimilation sur la même erreur issue de l'expérience de contrôle (OL):

$$r = \frac{\sum_{i} (x_{ASSIM}^{i} - x_{REF}^{i})^2}{\sum_{i} (x_{OL}^{i} - x_{REF}^{i})^2} \quad (4.1)$$

avec $x = T_{2m}$, $HU_{2m}$ et $LE$. Une valeur de $r$ égale à un indique que l'assimilation n'améliore pas par rapport au contrôle et une valeur de zéro que l'assimilation s'est parfaitement ajustée à la référence. Pour l'ensemble des assimilations on a une décroissance rapide durant les 10 premiers jours. Le filtre EnKF donne les meilleurs résultats en début de mois, alors que c'est le 2DVAR qui fournit les meilleures réductions d'erreurs sur toute la période. On note au cours des 8 premiers jours la supériorité de l'OI_MF pour réduire les erreurs sur $T_{2m}$ et aussi $HU_{2m}$ par rapport à l'OI_EC grâce aux importantes corrections appliquées à la température de surface même si l'erreur vient de l'eau du sol. Cette meilleure prévision des variables à 2 mètres n'a pas d'impact sur la valeur de $r$ pour le flux d'évaporation, les méthodes statiques gardant des erreurs plus élevées que les approches dynamiques. On note une certaine dégradation des variables à 2 mètres avec l'OI_MF juste après les épisodes pluvieux.

## 4.2. Assimilation de $w_g$

Dans un premier temps, j'ai effectué une assimilation de type "2D-Var simplifié" sur une fenêtre temporelle de 9 jours en considérant des observations une fois tous les trois jours (correspondant au temps de revisite des futurs instruments embarqués sur des satellites en orbite polaire). L'assimilation s'est poursuivie durant le mois d'Août afin de compléter la dernière période de 9 jours. Les résultats ne sont toutefois présentés que pour Juillet (convergence obtenue après un mois). En ce qui concerne la précision de la mesure pour $w_g$, j'ai considéré une valeur de 0.1 en SWI ce qui pour le sol choisi conduit à une incertitude en eau volumique de 0.009 m³/m³ (il faudrait prendre une valeur de 0.4 en SWI pour atteindre la précision attendue des données SMOS). L'annexe 3 présente une étude de sensibilité de l'analyse à la spécification de cette erreur. Comme précédemment, les expériences perturbées partent de conditions initiales au point de flétrissement pour les deux réservoirs et on introduit une erreur sur le forçage des précipitations (réduction d'un facteur deux par rapport à l'expérience de référence).

La Figure 4.15 montre la capacité de l'assimilation à converger en deux cycles (18 jours) vers une valeur proche de la simulation de référence pour les deux réservoirs. Des corrections non-négligeables sont appliquées par la suite montrant que l'assimilation de $w_g$ conduit à des sensibilités non nulles
<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /><figcaption>Figure 4.12: Identique à la Figure 4.1 mais pour une assimilation basée sur un filtre de Kalman d'ensemble avec une fenêtre temporelle de 6 heures.</figcaption></figure>
<figure><img src="image_34.png" /></figure>

Figure 4.13: Ecart-type d’erreur de prévision à 6 heures (ébauche) pour $w_g$ et $w_2$ estimé par un filtre de Kalman d’ensemble assimilation des observations $T_{2m}$ et $HU_{2m}$ toutes les 6 heures en Juillet 2002.

même lorsque le sol est proche de la capacité au champ car ce sont les processus de diffusion qui les définissent plutôt que le lien avec l’évapotranspiration. Ce résultat montre tout l’intérêt d’assimiler des observations plus directement (et linéairement) reliées à la variable que l’on cherche à initialiser.

L’intérêt d’utiliser des fenêtres d’assimilation longues pour restituer l’eau du sol profond à partir de mesures de l’eau du sol superficiel a été mis en évidence par Calvet et al. (1998) et Calvet et Noilhan (2000). Une différence importante par rapport à la méthode proposée ici est de ne pas utiliser d’information a-priori (ébauche fournie par une prévision à courte échéance). Le cyclage proposé ici (l’analyse au temps $t$ sert de condition initiale pour définir l’ébauche au temps $t + 1$) et utilisé dans la pratique de l’analyse atmosphérique permet de prendre en compte dans le terme d’ébauche les observations passées. Il faut également noter qu’une analyse sans terme d’ébauche évite la spécification des erreurs d’observations puisque les observations ne sont pas combinées à une autre source d’information.

J’ai affectué une assimilation de $w_g$ avec un cycle de 3 jours (une seule observation disponible à la fin de la période, ce qui est une configuration optimale, comme précédemment pour l’extraction de l’information contenue dans la dimension temporelle). Les résultats montrent que la convergence est plus rapide car les corrections au réservoir profond sont plus fréquentes (Figure 4.16). Il y aurait un avantage à considérer une fenêtre longue si le Jacobien après 9 jours était très différent de celui obtenu après 3 cycles de 3 jours. Dans l’atmosphère, l’existence de modes instables conduit à des Jacobiens “evolués” significativement différents de ceux estimés à un instant donné (structures baroclines, amplitudes plus importantes). La dimension temporelle n’est pas si intéressante pour le sol où, avec le temps, les trajectoires ont tendance à converger plutôt qu’à diverger.

S’il est nécessaire de considérer la dimension temporelle pour obtenir un Jacobien $\partial w_g^{t+\Delta t} / \partial w_2^t$ non nul (i.e. $\Delta t > 0$), l’échelle de temps, fournie dans les équations de ISBA par la constante de
<figure><img src="image_35.png" /><figcaption>Figure 4.14: Réduction de l’erreur quadratique moyenne des assimilations par rapport à un contrôle et à une référence pour la température à 2 mètres (haut), l’humidité spécifique à 2 mètres (milieu) et le flux d’évaporation (bas). Les observations considérées pour l’assimilation sont $T_{2m}$ et $HU_{2m}$.</figcaption></figure>
<figure><img src="image_36.png" /></figure>

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

Figure 4.15: Evolution de l'eau du sol superficiel (haut) et du sol profond (bas) au cours d'une as- simulation basée sur un 2D-Var simplifié avec une fenêtre de 9 jours et une observation de $w_g$ tous les trois jours.
<figure><img src="image_38.png" /></figure>

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

Figure 4.16: Evolution de l'eau du sol superficiel (haut) et du sol profond (bas) au cours d'une as- simulation basée sur un 2D-Var simplifié avec une fenêtre de 3 jours et une observation de $w_g$ tous les trois jours.
Jacobians of observation operator

<figure><img src="image_40.png" /><figcaption>Figure 4.17: Evolution du Jacobien $\partial w_g^{t+\Delta t} / \partial w_2^t$ entre le 1er Juillet 2002 et le 5 Août 2002 avec $\Delta t=3$ jours pour le "2DVAR 3d" et $\Delta t=3, 6$ et $9$ jours pour le "2DVAR 9d".</figcaption></figure>

rappel $\tau/C_2$, doit plutôt être l'ordre de deux jours. Notamment si une pluie apparaît durant l'intervalle $\Delta t$, l'eau du sol superficiel va complètement perdre la mémoire de l'eau du sol profond conduisant à des Jacobiens très faibles. C'est bien ce que l'on observe sur la Figure 4.17 qui compare les Jacobiens obtenus avec un cycle de 3 jours et un cycle de 9 jours. Les Jacobiens calculés sur une période $\Delta t$ de 3 jours sont presque systématiquement plus élevés. C'est particulièrement le cas au jour 15 où les conditions initiales pour un $\Delta t$ de 3 jours se situent dans un période de ciel clair (valeur de 0.92) alors que pour un $\Delta t$ de 6 jours les conditions initiales correspondent au fort épisode pluvieux du jour 9 conduisant à une valeur de 0.05 pour ce même Jacobien.

Finalement, l'intérêt d'une technique d'assimilation sur une fenêtre d'une dizaine de jours telle que proposée par Calvet et Noilhan (2000) ne se justifie pas lorsqu'on effectue une analyse qui combine les observations présentes avec une ébauche ayant déjà "vu" les observations passées. Autant pour des assimilations hydrologiques ce type de fenêtres longues pourrait être envisagé autant pour des applications météorologiques il semble difficile d'attendre dix jours pour pouvoir fournir un meilleur état hydrique des surfaces continentales, sachant qu'au cours de cette période d'autres observations informatives sur la surface seront disponibles.

Les résultats avec un EnKF pour assimiler $w_g$ montrent sa capacité à corriger rapidement à la fois le sol profond et le sol superficiel (Figure 4.18). La convergence est excellente après 10 jours. L'analyse de $w_g$ est en effet bien meilleure avec le filtre EnKF qu'avec les autres méthodes aussi bien en début qu'en fin de période. Ce résultat s'explique par l'introduction de corrélations fortes (> 0.90) dans la matrice B entre les erreurs du réservoir superficiel et celles du réservoir profond, conduisant à corriger les deux réservoirs de manière cohérente. On constate ici le bénéfice du cyclage de la matrice
B. Sur les deux mois d'assimilation, les erreurs de l'ébauche (Figure 4.19) évoluent de manière similaire à ce qui a été constaté avec ce filtre en assimilant les observations de $T_{2m}$ et $HU_{2m}$. On note des valeurs fortes de $\sigma^b$ pour $w_g$ en situations de ciel clair et des valeurs qui diminuent pour $w_2$ suite aux corrections de l'analyse. Lors d'événements pluvieux, l'erreur de prévision pour $w_g$ prend des valeurs faibles (car quelque soient les conditions initiales, le sol superficiel va être très humide) et une augmentation de $\sigma^b$ pour $w_2$ car l'analyse ne peut pas corriger ce réservoir et seul le facteur d'inflation agit pour augmenter cet écart-type. Pour $w_2$, l'erreur de prévision fluctue entre 0.009 m³/m³ (valeur initiale) et 0.005 m³/m³ autour d'une valeur moyenne située de l'ordre de 0.007 m³/m³.

### 4.3. Assimilation combinée de $T_{2m}$, $HU_{2m}$ et $w_g$

L'approche séquentielle de type EKF a également été testée pour préparer une assimilation combinant différents types d'observations. Dans cette méthode, je considère une fenêtre d'assimilation de 6-h au cours de laquelle je calcule le Jacobien suivant :

$$H_i = \frac{\partial w_g^{i+1}}{\partial \mathbf{x}^i} \qquad (4.2)$$

comme les observations de $w_g$ ne sont pas disponibles toutes les 6 heures, j'utilise le modèle linéaire-tangent $M_i$ pour propager dans le temps ce Jacobien :

$$M_i = \frac{\partial \mathbf{x}^{i+1}}{\partial \mathbf{x}^i} \qquad (4.3)$$

Ainsi au moment où une observation est disponible on peut exprimer le Jacobien par rapport aux conditions initiales au début d'une fenêtre longue (3 jours), en fonction des Jacobiens calculés sur les fenêtres courtes (6 heures), soit :

$$H = \frac{\partial w_g^n}{\partial \mathbf{x}^0} = \left[ \frac{\partial w_g^n}{\partial \mathbf{x}^{n-1}} \right] \prod_{i=1}^{n-1} M_{i-1} \qquad (4.4)$$

Tous les trois jours, l'opérateur $M_i$ est réinitialisé à la matrice identité. Le produit des opérateurs linéaires tangents oblige à avoir une estimation précise des dérivées en différences finies. J'ai choisi la taille des perturbations initiales assez faible (voir section 2) de façon à retrouver le même résultat que pour une intégration longue. Le résultat de cette assimilation est très voisin de celui du 2DVAR avec comme seule différence l'application des corrections à la fin de la période de 3 jours plutôt qu'au début (Figure 4.20). Cela donne l'apparence d'un moins bon résultat durant les 15 premiers jours de Juillet mais les valeurs d'eau du sol sont identiques à la fin de chaque période après analyse. De plus le bon résultat en début de période avec le 2DVAR ne peut s'obtenir qu'au bout de trois jours lorsqu'on dispose d'une observation. Par contre, l'analyse EKF permet d'assimiler d'autres observations sur la période de 3 jours et donc de corriger la variable de contrôle plus tôt qu'avec le 2DVAR. C'est cet aspect que nous examinons maintenant, à savoir la combinaison des observations SYNOP disponibles toutes les 6 heures (et même plus fréquemment comme par exemple avec le réseau RADOME) avec les observations peu fréquentes fournies par les instruments embarqués sur satellites.

J'ai considéré un cycle court de 6 heures où je calcule à la fois les Jacobiens pour les observations à 2 mètres et pour $w_g$. Comme l'eau du sol profond est la variable qui nous intéresse principalement et
qu'elle évolue lentement dans le temps, l'incrément au début de la fenêtre est appliqué au moment où l'observation de $w_g$ est disponible. On combine ici l'approche séquentielle du EKF avec l'utilisation d'un Jacobien "évolué" du 2D-Var simplifié. Le comportement de cette assimilation hybride est très similaire à l'analyse EKF sur 3 jours mais permet d'avoir des analyses plus fréquemment (toutes les 6 heures) compatibles avec les cycles d'assimilation pour l'atmosphère en utilisant les observations synoptiques quand elles sont disponibles, et donc finalement de permettre au système de converger plus rapidement vers l'eau du sol de référence (Figure 4.21). On constate que l'eau du sol superficiel a elle aussi été très rapidement corrigée entre les jours 5 et 15 ce que ne réussissaient à faire les analyses 2DVAR et EKF. La correction maximale est apportée au jour 3 avec un incrément de 50 mm (soit une correction de 13 mm/jour), ce qui, tout en étant élevé, reste acceptable compte-tenu de la forte valeur de l'innovation à cet instant (réservoirs proches de la saturation dans la référence). Comme pour les assimilations précédentes, le comportement du sol superficiel après le jour 15 est le résultat du forçage hydrique (sous-estimé par construction).

Finalement les réductions d'erreurs pour les paramètres atmosphériques montrent que l'assimilation de $w_g$ permet d'améliorer la prévision à court terme de ces paramètres de manière identique à l'assimilation des variables à 2 mètres sauf pour le 2DVAR avec une période de 9 jours qui présente des performances nettement inférieures (Figure 4.22). L'amélioration initiale est ainsi progressive-ment perdue au cours de l'intégration de 9 jours. Le gain sur les réductions d'erreurs de l'assimilation conjointe des deux types d'observations apparaît clairement. On note avec l'approche EKF qu'il faut attendre 3 jours au début du mois pour avoir une première réduction d'erreurs.
<figure><img src="image_41.png" /></figure>

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

Figure 4.18: Evolution de l'eau du sol superficiel (haut) et du sol profond (bas) au cours d'une assimilation basée sur un filtre de Kalman d'ensemble (ENKF) avec une fenêtre de 3 jours et une observation de $w_g$ tous les trois jours.
<figure><img src="image_43.png" /><figcaption>Figure 4.19: Ecart-types d’erreur de prévision pour $w_g$ et $w_2$ estimés par un filtre de Kalman d’ensemble sur la période Juillet-Août 2002 assimilant une observation de $w_g$ tous les 3 jours.</figcaption></figure>
<figure><img src="image_44.png" /></figure>

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

Figure 4.20: Evolution de l'eau du sol superficiel (haut) et du sol profond (bas) au cours d'une as- simulation basée sur un filtre de Kalman étendu (EKF) simplifié avec une fenêtre de 3 jours et une observation de $w_g$ tous les trois jours.
<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /><figcaption>Figure 4.21: Comme la Figure 4.20 mais en assimilant en plus les observations à 2 mètres toutes les 6 heures</figcaption></figure>
<figure><img src="image_48.png" /><figcaption>Figure 4.22: Réduction de l’erreur quadratique moyenne des assimilations par rapport à un contrôle et à une référence pour la température à 2 mètres (haut), l’humidité spécifique à 2 mètres (milieu) et le flux d’évaporation (bas). Les observations considérées pour l’assimilation sont $w_g$ (sauf pour une version du EKF où les valeurs de $T_{2m}$ et $HU_{2m}$ ont aussi été assimilées)</figcaption></figure>
# Chapitre 5

# Discussion et conclusions

J'ai présenté dans cette note une évaluation à l'échelle locale de plusieurs méthodes d'assimilation pour initialiser les variables pronostiques du schéma de surface ISBA. Il s'agit d'expériences qui, tout en utilisant un forçage atmosphérique réaliste (ERA40), ne considèrent que des observations "simulées" extraites d'une intégration de référence.

Un premier aspect était d'examiner la pertinence d'une assimilation découplée de l'atmosphère dans la mesure où il est envisagé de développer un système d'assimilation de surface multi-modèles et multi-observations dans le module de surface externalisé SURFEX. Cet aspect est surtout important pour l'assimilation d'observations de surface à 2 mètres dans la mesure où la sensibilité par rapport au premier niveau atmosphérique est nulle en mode forcé. Pour cela j'ai considéré un niveau de forçage situé à 50 m au dessus du sol. Dans ce contexte j'ai comparé l'analyse dans le sol utilisée à Météo-France et au CEPMMT basée sur une interpolation optimale "locale" avec des méthodes qui calculent dynamiquement les coefficients optimaux (EKF, 2DVAR, EnKF). J'ai montré, sur la période considérée (Juillet 2002, prairie aux USA), que toutes les méthodes donnent des résultats très semblables en terme de convergence vers la référence pour l'eau du sol profond. La comparaison entre les coefficients dynamiques et ceux issus de formulations analytiques a montré que les ordres de grandeurs sont similaires sauf pour l'eau du sol superficiel pour lequel la formulation de Météo-France est trois fois plus forte. Les corrections de la température de surface sont dix fois plus fortes avec les OI statiques. La modulation par le cycle diurne utilisée au CEPMMT est confirmée par les coefficients dynamiques. Toutefois l'application systématique des corrections de la $T_s$ avec le système OI_MF est bénéfique pour réduire les erreurs de $T_{2m}$ et $HU_{2m}$ à court-terme. Dans une évaluation à venir je montrerai que le système d'analyse de Météo-France se compare moins favorablement à celui du CEPMMT à cause de fortes corrections dans les zones montagneuses et d'un terme de relaxation assez fort vers la climatologie Global Soil Wetness Project (GSWP) à un degré de résolution. Pour l'eau du sol, la modulation empirique réduisant les coefficients en conditions pluvieuses et nocturnes se retrouve très bien sur les coefficients dynamiques. Les coefficients dynamiques sont cependant plus élévés en début de période (là où l'état du sol est le plus loin de la référence) qu'en fin de période. Les conclusions de Balsamo et al. (2007), indiquant des coefficients dynamiques en mode couplé et découplé comparables, sont ainsi confirmées et nous reforcent dans l'idée de pouvoir assimiler efficacement les observations à 2 mètres avec le système SURFEX. Certaines non-linéarités ont été mises en évidence pour le coefficient $\alpha_1$ mais leur impact n'apparaît pas dommageable au fonctionnement de l'assimilation. L'utilisation d'un filtre de Kalman d'ensemble à permis d'identifier certaines difficultés liées à la spécification des erreurs de modèle nécessaires pour
augmenter la dispersion de l'ensemble. De même, l'ajout d'un terme d'inflation empirique toutes les 6 heures s'est avéré nécessaire pour maintenir les variances de l'ébauche à un niveau suffisamment élevé pour que les observations aient un impact dans l'assimilation.

Une des limitations de l'approche de type "OI statique" est de ne pas pouvoir prendre en compte facilement de nouveaux types d'observations. En utilisant les approches dynamiques (EKF, 2DVAR, EnKF), j'ai examiné l'assimilation de l'eau du sol superficiel qui pourrait être fournie par des instruments embarqués sur satellites dans un futur proche. Cette étude se rapproche du travail de Muñoz-Sabater et al. (2007) qui ont utilisé des données réelles (campagne de mesures SMOSREX) sur une période de 4 ans. Dans un premier temps, un "2DVAR simplifié" sur 9 jours avec une assimilation tous les trois jours a été examiné, et la convergence de la méthode vers l'eau du sol de référence démontrée. Toutefois cette convergence n'est possible que si l'erreur d'observation pour $w_g$ est suffisamment faible (voir la discussion en Annexe 3). Une assimilation sur une fenêtre de 3 jours converge mieux. L'analyse des Jacobiens $\partial w_g / \partial w_2$ a montré que les valeurs à 6 et 9 jours sont en général plus faibles que la valeur à 3 jours à cause de la perte de sensibilité aux conditions initiales lors de pluies. L'intérêt d'une fenêtre d'assimilation longue, qui se justifiait dans les études de Calvet et al. (1998) et Calvet et Noilhan (2000) à cause de l'absence de terme d'ébauche, a plutôt un effet négatif dans un cadre d'assimilations plus réalistes (cyclage des analyses).

La convergence de l'eau du sol profond s'est avérée de meilleure qualité avec un EnKF grâce aux corrélations implicites dans la matrice B entre les erreurs sur $w_2$ et celles sur $w_g$. Dans la mesure où la matrice B n'a pas été cyclée dans les approches EKF et 2DVAR que j'ai mises en œuvre il faudra examiner cet aspect à l'avenir. De même, on constate avec le filtre EnKF une très forte modulation des erreurs de l'ébauche pour $w_g$ avec le forçage de pluie (très forte réduction quand le sol superficiel est très humide). Toutefois la variance de $w_g$ joue assez peu sur la correction du réservoir profond.

La possibilité de combiner l'assimilation des observations à 2 mètres disponibles toutes les 6 heures et des humidités superficielles disponibles tous les 3 jours a été étudiée avec un filtre EKF dans lequel le Jacobien estimé toutes les 6 heures est propagé dans le temps avec le modèle linéaire-tangent. Cet aspect n'avait jamais été examiné par le passé mais requiert une attention particulière pour les applications opérationnelles dans un futur proche. Cela permet d'avoir des analyses toutes les 6 heures tout en bénéficiant de la dimension temporelle permettant de relier les variables $w_g$ et $w_2$. Cette fenêtre de 3 jours pourrait être réduite et/ou on pourrait utiliser un 2DVAR sur cette fenêtre raccourcie (solution finale probablement identique mais méthodologie plus rigoureuse puisqu'avec le filtre EKF on applique les incrément au "mauvais instant").

L'approche "EKF simplifié" va être développée dans SURFEX et dans un premier temps utilisée pour assimiler les $T_{2m}$ et $HU_{2m}$ dans un cadre identique à ce qui est fait opérationnellement avec le modèle ARPEGE. Les aspects relatifs à l'analyse en hiver devront être examinés attentivement, notamment les non-linéarités liées au gel/dégel de l'eau dans le sol. L'étude de l'approche EnKF devra être poursuivie pour améliorer sa définition (spécification des diverses erreurs) et en analyser plus en détail les résultats. Certains problèmes comme la persistence de la divergence du filtre avec un grand nombre de membres et l'absence de sensibilité de sa réponse à des perturbations du forçage devront être compris, en collaborant avec des spécialistes de cette technique. Le problème de la spatialisation des observations avant introduction dans l'analyse du sol sera abordé dans un cadre idéalisé afin d'examiner la future méthodologie à mettre en œuvre pour l'utilisation des observations satellitaires. Finalement, la deuxième partie de ce rapport va présenter une évaluation de l'analyse dans le sol opérationnelle à Météo-France sur la zone Europe pour l'été 2006 où la méthode "EKF simplifié" va être initialement testée.
Références

Balsamo, G., F. Bouyssel, and J. Noilhan, 2004: A simplified bi-dimensional variational analysis of soil moisture from screen-level observations in a mesoscale numerical weather prediction model. *Quart. J. Roy. Meteor. Soc.*, **130**, 895-915.

Balsamo, G., J.-F. Mahfouf, S. Bélair, and G. Deblonde, 2006: A global root-zone soil moisture analysis using simulated L-band brightness temperature in preparation for the HYDROS satellite mission. *J. Hydrometeor.*, **7**, 1126-1146.

Balsamo, G., J.-F. Mahfouf, S. Bélair, and G. Deblonde, 2007: A land data assimilation system for soil moisture and temperature : an information content study. *J. Hydrometeor*. (accepté pour publication)

Bouttier, F., J.-F. Mahfouf, and J. Noilhan, 1993: Sequential assimilation of soil moisture from atmospheric low-level parameters. Part I: Sensitivity and calibration studies. *J. Appl. Meteor.*, **32**, 1335-1351.

Bouyssel, F., V. Cassé, and J. Pailleux, 1999: Variational surface analysis from screen-level atmospheric parameters. *Tellus*, **51A**, 453-468.

Calvet, J.-C., J. Noilhan, and P. Bessemoulin, 1998 : Retrieving the root-zone soil moisture from surface soil moisture or temperature estimates: a feasibility study based on field measurements. *J. Appl. Meteor.*, **37**, 371-386.

Calvet, J.-C., and J. Noilhan, 2000: From near-surface to root-zone soil moisture using year-round data. *J. Hydrometeor.*, **1**, 393-411.

Crow, W.T., and E.F. Wood, 2003 : The assimilation of remotely sensed brightness temperature imagery into a land surface model using ensemble Kalman filtering: A case study based on ESTAR measurements during SGP97. *Adv. Water Resour.*, **26**, 137-149.

Douville, H., P. Viterbo, J.-F. Mahfouf, and A.C.M. Beljaars, 2000: Evaluation of optimal interpolation and nudging techniques for soil moisture analysis using FIFE data. *Mon. Wea. Rev.*, **128**, 1733-1756.

Evensen, G., 2003 : The ensemble Kalman Filter: Theoretical formulation and practical implementation. *Ocean Dyn.*, **53**, 343-367.

Geleyn, J.-F., 1988: Interpolation of wind, temperature, humidity values from model levels to the height of measurement. *Tellus*, **40A**, 347-351.
Germain, M.J., 1990: Validation d'une paramétrisation des échanges de surface sur différents types de végétation. Note de travail de l'ENM, 140 pp.

Giard, D., and E. Bazile, 2000: Implementation of a new assimilation scheme for soil and surface variables in a global NWP model. *Mon. Wea. Rev.*, **128**, 997-1015.

Habets, F., A Boone, J.-L.Champeaux, P. Etchevers, L. Franchisteguy, E. Leblois, E. Ledoux, P. Le Moigne, E. Martin, S. Morel, J. Noilhan, P. Quintana-Seguì, F. Rousset-Regimbeau, and P. Viennot, 2007: The SAFRAN-ISBA-MODCOU hydrometeorological model applied over France. *J. Geophys. Res.* (soumis).

Hess, R., 2001: Assimilation of screen-level observations by variational soil moisture analysis. *Meteor. Atmos. Phy.*, **77**, 145-154.

Houtekamer, P. and H.L. Mitchell, 1998: Data assimilation using an Ensemble Kalman filter technique. *Mon. Wea. Rev.*, **126**, 796-811.

Mahfouf, J.-F., 1991: Analysis of soil moisture from near-surface parameters: A feasibility study. *J. Appl. Meteor.*, **30**, 506-526.

Mahfouf, J.-F., and J. Noilhan, 1991: Comparative study of various formulations of evaporation from bare soil using in-situ data. *J. Appl. Meteor.*, **30**, 1354-1365.

Mahfouf, J.-F., A.O. Manzi, J. Noilhan, H. Giordani, and M. Déqué, 1995: The land surface scheme ISBA within the Météo-France climate model ARPEGE. Part I: Implementation and preliminary results. *J. Climate*, **8**, 2039-2057.

Mahfouf, J.-F., H. Douville, P. Viterbo, A. Beljaars, and S. Saarinen, 2000: A revised land surface analysis scheme in the Integrated Forecasting System. *ECMWF Newsletter*, **88**, 8-13.

Mascart, P., J. Noilhan, and H. Giordani, 1995: A modified parameterization of the surface layer flux-profile relationships using different roughness length values for heat and momentum. *Bound. Layer. Meteor.*, **72**, 331-344.

Mitchell, K.E. et al-auteurs, 2004: The multi-institution North American Land Data Assimilation System (NLDAS): Utilizing multiple GCIP products and partners in a continental distributed hydrological modeling system. *J. Geophys. Res.*, **109**, D07S90, doi:10.1029/2003/2003JD003823.

Muñoz-Sabater, J., L. Jarlan, J.-C. Calvet, F. Bouyssel, and P. De Rosnay, 2007: From near-surface to root-zone soil moisture using different assimilation techniques. *J. Hydrometeor.*, **8**, 194-206.

Noilhan, J., and J.-F. Mahfouf, 1996: The ISBA land surface parameterization scheme. *Global Planet. Change*, **13**, 145-159.

Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, 1992: Numerical recipes in FORTRAN. The art of scientific computing. Cambridge University Press, 963 pp.

Reichle, R.H., J.P. Walker, R.D. Koster, and P.R. Houser, 2002: Extended versus ensemble Kalman filtering for land data assimilation. *J. Hydrometeor.*, **3**, 728-740.
Rhodin, A., F., Kucharski, U. Callies, and D.P. Eppel, 1999: Variational soil moisture analysis from screen-level atmospheric parameters: Application to a short-range weather forecast model. *Quart. J. Roy. Meteor. Soc.*, **125**, 2427-2448.

Seuffert, G., H. Wilker, P. Viterbo, M. Drusch, and J.-F. Mahfouf, 2004: The usage of screen-level parameters and microwave brightness temperature for soil moisture analysis. *J. Hydrometeor.*, **5**, 516-531.

Staley, D.O., and G.M. Jurica, 1972: Effective atmospheric emissivity under clear skies. *J. Appl. Meteor.*, **11**, 349-356.

Uppala, S., P.W. Kallberg, A.J. Simmons, U. Andrae, V. da Costa Bechtold, M. Fiorino, J.K. Gibson, J. Haseler, A. Hernandez, G.A. Kelly, X. Li, K. Onogi, S. Saarinen, N. Sokka, R.P. Allan, E. Andersson, K. Arpe, M.A. Balmaseda, A.C.M. Beljaars, L. van de Berg, J. Bidlot, N. Bormann, S. Caires, F. Chevallier, A. Dethof, M. Dragosavac, M. Fisher, M. Fuentes, S. Hageman, E. Holm, B.J. Hoskins, L. Isaksen, P.A.E.M. Janssen, R. Jenne, A.P. McNally, J.-F. Mahfouf, J.-J. Morcrette, N.A. Rayner, R.W. Saunders, P. Simon, A. Sterl, K.E. Trenberth, A. Untch, D. Vasiljevic, P. Viterbo and J. Woollen, 2005: The ERA-40 reanalysis, *Quarterly Journal of the Royal Meteorological Society*, **131**, 2961-3012.

Zhou, Y., D. McLaughlin, and D. Entekhabi, 2006: Assessing the performance of the Ensemble Kalman Filter for land surface data assimilation. *Mon. Wea. Rev.*, **134**, 2128-2142.
# Annexe 1: Constitution du forçage atmosphérique

J'ai utilisé les données de la réanalyse de 40 ans du Centre Européen (ERA40: Uppala et al., 2005) disponibles toutes les 6 heures sur une grille latitude/longitude de 2.5° pour les périodes de Juillet et Août 2002. On dispose des paramètres suivants :

* Composante zonale du vent à 10 m : $U_{10m}$

* Composante méridienne du vent à 10 m : $V_{10m}$

* Température de l'air à 2 m : $T_{2m}$

* Température du point de rosée à 2 m : $T_d$

* Rayonnement solaire descendant en surface accumulé sur 6 heures : $\sum R_s$

* Rayonnement thermique descendant en surface accumulé sur 6 heures : $\sum R_l$

* Précipitation accumulée sur 6 heures : $\sum PR$

* Pression réduite au niveau de la mer : $P_{mer}$

* Géopotentiel de surface : $gZ_{orog}$

* Longueurs de rugosité dynamique et thermique : $Z_0, Z_{0H}$

A partir de cet ensemble de paramètres, j'ai constitué les champs suivants à une fréquence horaire (pour être ensuite interpolées linéairement tous les pas de temps du modèle) :

* Composante zonale du vent à 50 m : $U_{50m}$

* Composante méridienne du vent à 50 m : $U_{50m}$

* Température de l'air à 50 m : $T_{50m}$

* Humidité spécifique de l'air à 50 m : $q_{50m}$

* Flux instantané de rayonnement solaire descendant en surface : $R_s$

* Flux instantané de rayonnement thermique descendant en surface : $R_l$

* Flux de précipitation instantané : $PR$
• Pression de surface : $P_s$

On calcule dans un premier temps la température $T_s$ du niveau de la mer en utilisant un gradient vertical standard $\Gamma$=-6.5K/km

$$T_s = T_{2m} + \Gamma(Z_{2m} + Z_{orog}) \quad (5.1)$$

où $Z_{orog}$ est l'altitude du relief. Ensuite la relation hydrostatique et l'équation des gaz parfaits permet le calcul de la pression de surface $P_s$ :

$$P_s = P_{mer} \left[ \frac{T_{2m}}{T_s} \right]^{R_d \Gamma / g} \quad (5.2)$$

Le vent au niveau 50 m est interpolé à partir du vent à 10 m en utilisant un profil logarithmique (couche neutre) :

$$U_{50m} = U_{10m} \frac{\log \left[ \frac{Z_{50m} + Z_0}{Z_0} \right]}{\log \left[ \frac{Z_{10m} + Z_0}{Z_0} \right]} \quad (5.3)$$

La température au niveau 50 m est obtenue par simple extrapolation avec un gradient standard $\Gamma$ :

$$T_{50m} = T_{2m} - \Gamma(Z_{50m} - Z_{2m}) \quad (5.4)$$

L'humidité spécifique au niveau 50 m est obtenue en supposant l'humidité relative constante sur la couche 0-50 m:

$$q_{50m} = \frac{e_{sat}(T_d)}{e_{sat}(T_{2m})} \times q_{sat}(P_s - 10Z_{50m}, T_{50m}) \quad (5.5)$$

Les précipitations instantanées sont obtenues en divisant le cumul sur 6 heures par 21600 s. Pour le flux de rayonnement solaire l'échantillonage toutes les 6 heures ne permet pas de décrire correctement le cycle diurne, j'ai donc considéré une formule empirique permettant d'avoir la dépendance explicite avec l'angle zénithal $\mu$ :

$$R_s = S_0 \cos \mu T^{(1/\cos \mu)} \quad (5.6)$$

où $S_0 = 1370 \text{ W/m}^2$ est la constante solaire et $T = 0.7$ la transmittivité de l'atmosphère. Si la précipitation en surface est non nulle, cette valeur est multipliée par 0.4. Le rayonnement thermique descendant est également exprimé de manière empirique par la formule de Staley et Jurica (1972):

$$R_l = 0.67 [e_{sat}(P_s, T_{50m})]^{0.08} \sigma T_{50m}^4 \quad (5.7)$$

où $\sigma$ est la constante de Stefan-Boltzmann.

L'échantillonage temporel pour les variables atmosphériques ($U_{50m}$, $V_{50m}$, $T_{50m}$, $q_{50m}$, $P_s$) a été accru au pas de temps horaire en utilisant une interpolation de valeurs toutes les 6 heures au moyen de splines cubiques. Un lissage additionnel des flux radiatifs toutes les heures $R(i)$ a finalement été appliqué avec deux passes du filtre suivant :

$$\hat{R}(i) = 0.25R(i-1) + 0.5R(i) + 0.25R(i+1) \quad (5.8)$$

Ce dernier lissage est plutôt d'ordre esthétique son influence sur les résultats du schéma ISBA étant faible. La Figure 5.1 présente le forçage résultant au point choisi pour les tests d'assimilation (Juillet 2002). On constate que le site va présenter un forçage radiatif important en dehors de périodes pluvieuses permettant un couplage fort entre la surface et la basse atmosphère. Le flux thermique est
assez constant autour de 400 W/m² alors que le flux solaire est caractérisé par le cycle diurne et des maxima plus faibles en situations perturbées (par définition). Les précipitations sont réparties assez irrégulièrement au cours du mois avec un cumul mensuel de 100 mm. On notera un fort épisode au jour 9 ainsi qu'une période de 6 jours avec des précipitations importantes entre les jours 17 et 23. La période de ciel clair entre les jours 9 et 17 est caractérisée par une advection d'air sec et un module du vent élevé favorisant l'évaporation. Le module du vent reste assez faible au voisinage de 4 m/s sauf vers le jour 10 et en fin de période (jours 25 à 30). Le cycle diurne très régulier pour la température est la conséquence de l'interpolation par splines cubiques effectuée sur les valeurs horaires par la forte modulation du cycle diurne. L'effet de lissage temporel est moins évident sur l'humidité et le vent.
<figure><img src="image_49.png" /><figcaption>Figure 5.1: Série temporelle du forçage atmosphérique utilisé pour les expériences d'assimilation à partir des données de ERA40 extraites au point (40N, 275E) en Juillet 2002. Ce forçage horaire a été constitué à partir de données toutes les 6 heures suivant une procédure décrite dans le texte.</figcaption></figure>
# Annexe 2: Comparaison avec la version ISBA-2L de SURFEX

La version de ISBA développée pour cette étude est dérivée de la description faite dans Noilhan et Mahfouf (1996) à l'exception du flux d'évaporation sur sol nu pour lequel la formulation en résistance de sol utilisée au CEPMMT, a été considérée :

$$E_g = (1 - veg) \rho \frac{q_{sat}(T_s) - q_a}{R_{soil} + R_a} \qquad (5.9)$$

avec

$$R_{soil} = R_{min} \left[ \frac{w_{fc} - w_{wilt}}{w_2 - w_{wilt}} \right] \qquad (5.10)$$

et $R_{min}=50$ s/m. Mahfouf et Noilhan (1991) ont montré que les formulations en résistance et en humidité relative étaient très proches en régime diurne. La formulation en humidité relative nécessite un traitement particulier de la rosée sur sols très secs. La formulation en résistance est aussi moins sensible à la spécification de la rugosité. D'autres différences mineures avec la version ISBA-2L de SURFEX concernent la formulation de la tension de vapeur saturante en fonction de la température et d'autres aspects que je n'ai pas identifié (comme le détail des fonctions de stabilité dans la couche limite de surface). L'avantage de cette version est de gérer de manière plus flexible les entrées/sorties qui intervenissent plus fréquemment dans des cycles d'assimilation que dans de simples intégrations. Actuellement SURFEX ne lit et n'écrit que des fichiers NetCDF dans sa version "offline" nécessitant des étapes supplémentaires d'encodage/décodage. Le code est informatiquement plus efficace, ce qui pour des intégrations nombreuses requises par le filtre de Kalman d'ensemble, est un avantage certain. En effet un mois d'intégration en un point sur un ordinateur PC Pentium IV avec un pas de temps de 15 min nécessite 100 s avec SURFEX et 2 s avec la version simplifiée. Pour le forçage utilisé dans les expériences, je présente dans la Figure 5.2 l'évolution des quatre variables pronostiques du schéma en Juillet 2002 échantillonnées toutes les heures. La version simplifiée se compare bien avec SURFEX. On note deux différences principales : l'amplitude du cycle diurne pour $T_s$ est plus forte avec ISBA et pour $w_s$ la déplétion significative entre les jours 10 et 15 (période de ciel clair suite à un fort épisode pluvieux) est plus prononcée avec SURFEX. En ce qui concerne le bilan radiatif présenté sur la Figure 5.3, j'ai sélectionné une période de 10 jours (10-19 Juillet 2002) de façon à mieux mettre en évidence le cycle diurne. Les flux de rayonnement net et de chaleur latente se comparent très bien et ce sont les plus importants. Le maximum du flux de chaleur sensible est plus élevé avec ISBA d'environ 20 W/m² ce qui est cohérent avec les valeurs plus élevées de $T_s$ en milieu de journée. Le flux de chaleur dans le sol est celui qui présente le plus de différences entre les deux versions ce qui n'est pas étonnant dans la mesure où il est obtenu comme résidu des autres flux et reflète ainsi l'ensemble des différences. L'amplitude du maximum est plus élevée avec ISBA et apparaît en début
de matinée lorsque les flux turbulents restent faibles mais que le forçage radiatif augmente (on doit
donc passer du régime stable au régime instable plus tardivement). En conclusion, le schéma ISBA
simplifié reproduit de manière réaliste l'évolution des variables dans le sol (que l'on va chercher à
initialiser) et des composantes du bilan d'énergie.
<figure><img src="image_50.png" /><figcaption>Figure 5.2: Evolution des variables prognostiques du schéma ISBA ($T_s$, $T_2$, $w_g$, $w_2$) en Juillet 2002 forcé par les analyses ERA40 décrites dans l'Annexe 1. "ISBA" correspond à la version utilisée dans cette note et "SURFEX" à celle disponible dans le module externalisé (ISBA-2L).</figcaption></figure>
<figure><img src="image_51.png" /><figcaption>Figure 5.3: Evolution des différents termes du bilan d'énergie en surface simulés par le schéma ISBA entre le 10 et le 19 Juillet 2002 forcé par les analyses ERA40 décrites dans l'Annexe 1. "ISBA" correspond à la version utilisée dans cette note et "SURFEX" à celle disponible dans le module externalisé (ISBA-2L).</figcaption></figure>
# Annexe 3: Sensibilité à la spécification de l'erreur d'observation pour $w_g$

Je discute ici brièvement du problème de l'erreur associée à la mesure de l'eau du sol superficiel $w_g$ pour analyser $w_2$. L'équation d'analyse se réduit à la formulation scalaire suivante :

$$w_2^a = w_2^b + \left\{ \frac{\kappa[\sigma_{w2}^b]^2}{\kappa^2[\sigma_{w2}^b]^2 + [\sigma_{wg}^o]^2} \right\} (w_g^o - w_g^b) \quad (5.11)$$

$$\Delta w_2 = \chi \Delta w_g \quad (5.12)$$

La matrice de gain (réduite à un scalaire) $\chi$ est fonction du Jacobien $\kappa = \partial w_g / \partial w_2$ et du rapport des erreurs d'observation pour le contenu en eau du sol superficiel $\sigma_{wg}^o$ sur l'erreur de l'ébauche pour le contenu en eau dans la zone racinaire $\sigma_{w2}^b$. La Figure 4.17 montre que $\kappa$ varie entre 0.2 dans des situations peu informatives à des valeurs maximales de 1.0 (fin d'une période de ciel clair après une pluie). Cette valeur maximale se compare bien avec les estimés de Calvet et Noilhan (2000) obtenus dans des conditions similaires mais avec des perturbations d'amplitude finie ($\Delta w_g / \Delta w_2$). Les variations du coefficient $\chi$ en fonction de $\kappa$ sont présentées pour quatre valeurs du rapport $\sigma_{wg}^o / \sigma_{w2}^b$ (entre 1 et 4) sur la Figure 5.4. Une innovation de 0.1 m³/m³ correspondant à toute la plage de la réserve utile peut être considérée comme une valeur maximale (e.g. pluie sur sol sec). Pour un rapport de un entre les erreurs et une valeur $\kappa = 1$ on trouve une correction dans le sol profond de 0.05 m³/m³, soit 50 mm pour 1 mètre de sol (soit 13 mm/jour). Ce sont les valeurs maximales des corrections effectuées pour les assimilations présentées dans la Section 4b. Par contre en utilisant la précision SMOS, le rapport $\sigma_{wg}^o / \sigma_{w2}^b$ vaudrait 4, conduisant à une matrice de gain 10 fois plus faible, soit une correction journalière de 1.3 mm, conduisant à une convergence extrêmement lente du système. D'une part on est pénalisé par le fait que les corrections ne sont disponibles que tous les 3 jours, et d'autre part dès qu'une pluie apparaît la valeur de $\chi$ chute à cause de la diminution de la sensibilité exprimée par le coefficient $\kappa$. Dans l'étude de Muñoz-Sabater et al. (2007) le rapport $\sigma_{wg}^o / \sigma_{w2}^b$ vaut 6 conduisant à une correction de 0.9 mm par jour pour un sol d'un mètre d'épaisseur, et pose la question de l'influence des observations de $w_g$ dans leurs assimilations. L'étude de Balsamo et al. (2006) montre une correction mensuelle globale de 0.03 m³/m³ en assimilant des données de température de brillance en bande L ($T_b$) simulées avec une précision de 3 K, soit une correction de 1.0 mm/jour en accord avec notre estimation en utilisant la précision SMOS.<fnref n="1" /> Le poids de ces

<footnote n="1">¹Une étude de sensibilité de Balsamo et al. (2006) montre que $T_b$ varie d'environ 15 K pour une variation de SWI entre 0 et 1 (Fig. 1). On peut donc considérer (dans un régime linéaire) qu'une erreur en $T_b$ de 3K équivaut à une erreur en SWI de 0.2. Par ailleurs, ils ont imposé une erreur de l'ébauche de 0.05 SWI, conduisant à un rapport $\sigma_{wg}^o / \sigma_{w2}^b$ de 4.</footnote>
<figure><img src="image_52.png" /></figure>

Figure 5.4: Coefficient de la matrice de gain χ associé à l'incrément dans le réservoir profond $w_2$ en fonction de l'innovation pour l'eau du réservoir superficiel $w_g$. Sa dépendance est exprimée en fonction de la valeur du Jacobien $\partial w_g / \partial w_2$ de l'opérateur d'observation et du rapport des erreurs $\sigma_{wg}^o / \sigma_{w2}^b$. La courbe supérieure est correspond à un rapport de 1 et les courbes suivantes à des rapports de 2, 3 et 4 respectivement.

données sera faible dans des régions où les observations à 2 mètres seront disponibles comme cela a également été confirmé par Balsamo et al. (2007).

En considérant une erreur de l'ébauche deux fois plus importante $\sigma_{w2}^b = 0.2 \, SWI$, et une valeur de l'erreur d'observation $\sigma_{wg}^o = 0.4 \, SWI$ conforme à la précision attendue du produit SMOS, on constate, avec un 2DVAR simplifié sur 3 jours, qu'il faut 2 mois (Juillet et Août 2002) pour corriger de façon significative le contenu en eau du réservoir profond (Figure 5.5). Cette valeur forte de $\sigma_{w2}^b$ ne semble pas confirmée par le filtre EnKF qui l'estime dynamiquement et conduirait aussi à donner un poids trop fort aux observations à 2 mètres par rapport à ce qui a été présenté dans cette note.

On sait que la précision SMOS va dépendre du type de surface (végétation, présence d'eau libre dans le pixel, sols gelés, neige, ...) et donc cette valeur de 0.04 m³/m³ est plutôt générique. L'assimilation de l'eau du sol superficiel issue d'une inversion de températures de brillance va donc nécessiter une connaissance précise de l'erreur associée (spatialement et temporellement). Des efforts sont également nécessaires du côté de la modélisation pour améliorer la spécification des erreurs de l'ébauche. Ces deux conditions sont indispensables pour pouvoir assimiler avec succès les produits SMOS dans les modèles de PN du temps
<figure><img src="image_53.png" /><figcaption>Figure 5.5: Evolution de l'eau du sol superficiel (haut) et du sol profond (bas) au cours d'une as- simulation basée sur une méthode "2D-Var simplifiée" avec une fenêtre de 3 jours. Une observation de contenu en eau superficiel $w_g$ est assimilée tous les 3 jours. L'erreur d'observation est fixée à $0.4 \times \text{SWI}$ (soit $0.036 \text{ m}^3/\text{m}^3$) et les erreurs de l'ébauche valent $0.2 \times \text{SWI}$ pour chacun de réservoirs.</figcaption></figure>