# Initialisation de l'eau du sol dans les modèles de prévision. Tests préliminaires sur sol nu

Jean-François Mahfouf

► To cite this version:

Jean-Francois Mahfouf. Initialisation de l'eau du sol dans les modèles de prévision. Tests préliminaires sur sol nu. Note de travail du CNRM N°6, Centre National de Recherches Météorologiques. 1989. hal-04478334

HAL Id: hal-04478334

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

Submitted on 26 Feb 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
MINISTERE DE L'EQUIPEMENT, DU LOGEMENT,
DES TRANSPORTS ET DE LA MER

DIRECTION DE LA METEOROLOGIE NATIONALE

ETABLISSEMENT D'ETUDES ET DE RECHERCHES METEOROLOGIQUES

NOTE DE TRAVAIL

DU

CENTRE NATIONAL DE RECHERCHES METEOROLOGIQUES

N°6

INITIALISATION DE L'EAU DU SOL DANS
LES MODELES DE PREVISION
TESTS PRELIMINAIRES SUR SOL NU

Par

Jean-François MAHFOUF
# SOMMAIRE

1-Introduction 1

2-Les équations de l'analyse optimale

A-Position du problème 3

B-Solution générale 3

C-Remarques 5

3-Présentation des simulations de référence 7

4-Assimilation variationnelle de l'eau du sol

A-Présentation générale de l'assimilation quadridimensionnelle 9

B-Assimilation de l'état hydrique des sols 10

C-Résultats 12

5-Assimilation séquentielle de l'eau du sol

A-Méthode de Geleyn 15

B-Méthode de Geleyn modifiée 17

C-Interpolation optimale 18

6-Conclusion 22

Annexe 23

Références 24
# 1-Introduction

La nécessité pour les modèles atmosphériques de représenter de façon réaliste les transferts d'eau et d'énergie entre le sol et l'atmosphère n'est maintenant plus à prouver tant les études de sensibilité sur le sujet sont nombreuses (McCumber et Pielke, 1981; Mahfouf, 1986; Jacquemin, 1987; Lanicci et al., 1987; Benjamin et Carlson, 1986). C'est ainsi que depuis quelques années, on assiste à une recrudescence dans la littérature de nouveaux schémas d'interface sol-atmosphère (Wilson et al., 1987; Sellers et Dorman, 1987; Segal et al., 1988; Abramopoulos et al., 1988; Pinty et al., 1988; Noilhan et Planton, 1989).

Si l'on veut que ces nouveaux schémas puissent un jour remplir le rôle pour lequel ils ont été développés, c'est à dire améliorer la prévision du temps, il faudra être capable de répondre à la question suivante: "Comment initialiser les variables pronostiques des schémas de surface dans les modèles numériques ?". Le but de ce rapport est d'effectuer quelques investigations préliminaires dans ce sens.

La température et l'humidité des sols ne sont pas des quantités mesurées en routine par le réseau météorologique. Il faut donc accéder à ces quantités de manière indirecte. Quelques essais ont déjà été tentés en utilisant la température radiative de surface télédétectée (Wetzel et al., 1984; Taconet et al., 1986). L'approche retenue ici est d'essayer de tirer parti de l'information météorologique du réseau de surface pour déduire les valeurs les plus vraisemblables de l'eau et la température dans le sol. Elle s'appuiera sur le jeu de données de la campagne HAPEX-MOBILHY, nous permettant de disposer à la fois de mesures au voisinage de la surface et dans le sol.

Nous verrons tout d'abord que la solution de ce problème s'inscrit dans le cadre général de l'analyse optimale dont les équations ont été résumées de manière élégante par Lorenc (1986). A partir de là, disposant du schéma de surface de Noilhan et Planton (1989) et d'une version 1D du modèle de prévision opérationnel PERIDOT, nous explorerons deux types d'initialisation pour l'eau du sol. Le premier type se place dans la perspective des schémas d'analyse de la prochaine génération, utilisant pleinement, au travers du formalisme variationnel, la structure temporelle des observations et la dynamique du modèle numérique. Le second type, qui est de nature à pouvoir être introduit facilement dans le schéma d'analyse opérationnel actuel, consiste en une suite de
prévisions et de corrections statiques par une méthode de type interpolation optimale.

Les tests présentés seront relatifs à la journée du 19 mai 1986 à Bats. Il s'agit d'une journée de ciel clair sur sol nu. Notre objectif sera, à l'aide des méthodes précitées, de retrouver le contenu en eau du sol issu d'une expérience de référence, correspondant à la valeur observée. Dans un premier temps, il sera fait abstraction des erreurs d'observation.
## 2- Les équations de l'analyse optimale

On résume ici une présentation faite par Lorenc (1986,1988) qui fournit la solution générale du problème de l'analyse optimale.

### A-Position du problème

Le problème de l'analyse en prévision du temps consiste à déterminer les meilleures conditions initiales d'un modèle numérique c'est à dire celles conduisant à la meilleure prévision.

On se définit un état du modèle par un vecteur $\vec{x}$ de dimension $N_x$ et un ensemble de données d'observation par un vecteur $\vec{y}$ de dimension $N_y$. On a besoin de connaître la relation entre l'état du modèle $\vec{x}$ et les observations $\vec{y}$. On la représente sous la forme d'un opérateur explicite $K_n$ non linéaire tel que:

$$ \vec{y} = K_n(\vec{x}) \quad (1) $$

soit la meilleure estimation de $\vec{y}$ pour un $\vec{x}$ donné.

On peut alors redéfinir le problème de l'analyse comme trouver la meilleure valeur de $\vec{x}$ inver-
sant (1) étant donnés un ensemble d'observations $\vec{y}_0$, la connaissance d'erreurs d'observations et
d'une information préalable.

### B-Solution générale

La solution du problème précédent repose sur des arguments bayesiens (probabilités condition-
nelles). En effet, la connaissance de la densité de probabilité:

$$ P_a(\vec{x}) = P(\vec{x} = \vec{x}_t | \vec{y} = \vec{y}_0) $$

donnant la probabilité qu'a un état $\vec{x}$ de l'atmosphère d'être l'état vrai étant donné un ensemble
d'observations $\vec{y}_0$, contient toute l'information nécessaire à notre problème d'analyse.

Etant donnée cette densité de probabilité, la meilleure estimation $\vec{x}_a$ de $\vec{x}$ sera prise comme le
vecteur le plus vraisemblable dans l'espace des représentations possibles de l'atmosphère, soit

$$ \vec{x} = \vec{x}_a \text{ tel que } P_a(\vec{x}) \text{ soit maximum} $$

En utilisant le théorème de Bayes:

$$ P_a(\vec{x}) \propto P(\vec{y}_0 | \vec{x} = \vec{x}_t) P(\vec{x} = \vec{x}_t) \quad (2) $$
$P(\vec{x} = \vec{x}_t)$ englobe toute la connaissance sur l'état $\vec{x}$ avant que ne soit effectuée l'observation $\vec{y}_0$. Cette connaissance est disponible sous la forme d'une estimation préalable (prévision ou climatologie) de la valeur la plus vraisemblable de $\vec{x}$: $\vec{x}_b$. On écrit alors:

$$P(\vec{x} = \vec{x}_t) = P(\vec{x} | \vec{x}_b)$$

Du fait d'erreurs de mesure, l'observation $\vec{y}_0$ est différente de la vraie valeur $\vec{y}_t$. Cela conduit à introduire la probabilité conditionnelle $P(\vec{y}_0 | \vec{y} = \vec{y}_t)$.

Il faut aussi prendre en compte les erreurs de représentativité liées à l'opérateur $K_n$. Si $\vec{x}_t$ est la vraie valeur de $\vec{x}$, $K_n(\vec{x}_t)$ peut être différent de $\vec{y}_t$. On introduit donc la probabilité conditionnelle $P(\vec{y} | \vec{x} = \vec{x}_t)$.

La quantité $P(\vec{y}_0 | \vec{x} = \vec{x}_t)$ englobe à la fois les erreurs de représentativité et les erreurs instrumentales. Elle s'écrit:

$$P(\vec{y}_0 | \vec{x} = \vec{x}_t) = \int P(\vec{y}_0 | \vec{y} = \vec{y}_t) P(\vec{y} | \vec{x} = \vec{x}_t) d\vec{y}$$

L'équation (2) prend alors la forme:

$$P_a(\vec{x}) = P(\vec{y}_0 | \vec{x} = \vec{x}_t) P(\vec{x} | \vec{x}_b) / \int P(\vec{y}_0 | \vec{x}_t = \vec{x}_t) P(\vec{x}_t | \vec{x}_b) d\vec{x}_t$$

En définissant la meilleure analyse $\vec{x}_a$ comme l'estimation du maximum de vraisemblance, elle correspond à la valeur de $\vec{x}$ qui maximise la densité de probabilité $P_a(\vec{x})$. Ce qui est équivalent à minimiser la fonction de coût $J$ suivante:

$$J(\vec{x}) = - \ln\{P(\vec{y}_0 | \vec{x} = \vec{x}_t)\} - \ln\{P(\vec{x} | \vec{x}_b)\}$$

Considérons le cas particulier de fonctions de densité gaussiennes, c'est à dire:

$$P(\vec{y} | \vec{x} = \vec{x}_t) \propto \exp \left[ -\frac{1}{2} \{\vec{y} - K_n(\vec{x})\}^* F^{-1} \{\vec{y} - K_n(\vec{x})\} \right]$$

$$P(\vec{y} | \vec{y} = \vec{y}_t) \propto \exp \left[ -\frac{1}{2} \{\vec{y} - \vec{y}_0\}^* O^{-1} \{\vec{y} - \vec{y}_0\} \right]$$

$$P(\vec{x} | \vec{x}_b) \propto \exp \left[ -\frac{1}{2} \{\vec{x} - \vec{x}_b\}^* B^{-1} \{\vec{x} - \vec{x}_b\} \right]$$
$$P(\bar{y}_0 | \bar{x} = \bar{x}_t) \propto \exp \left[ -\frac{1}{2} \{\bar{y}_0 - K_n(\bar{x})\}^* (O + F)^{-1} \{\bar{y}_0 - K_n(\bar{x})\} \right]$$

d'où

$$J(\bar{x}) = \frac{1}{2} [(\bar{y}_0 - K_n(\bar{x}))^* (O + F)^{-1} (\bar{y}_0 - K_n(\bar{x})) + (\bar{x} - \bar{x}_b)^* B^{-1} (\bar{x} - \bar{x}_b)]$$

$B$, $O$ et $F$ sont les matrices de covariances associées à chacune des densités de probabilité. $(\bar{x})^*$ correspond à la transposée du vecteur $\bar{x}$.

Le problème d'analyse se présente donc comme la minimisation d'une fonctionnelle. Le premier terme représente l'ajustement aux observations et le second l'ajustement à l'estimation préalable (guess). Lorsque l'opérateur de liaison entre $\bar{x}$ et $\bar{y}$ est linéaire, l'estimation $\bar{x}_a$ minimisant $J(\bar{x})$ peut être obtenue explicitement (Lorenc, 1986):

$$\bar{x}_a = \bar{x}_b + BK^* (KBK^* + O + F)^{-1} \{\bar{y}_0 - K_n(\bar{x}_b)\} \quad (3)$$

Une des caractéristiques du problème d'analyse en prévision numérique est d'être sous déterminé ($N_x > N_y$). Cette indétermination est levée en introduisant de l'information supplémentaire visant à augmenter virtuellement $N_y$. Ainsi, $P(\bar{x} | \bar{x}_b)$ constitue une forme d'information, elle prend le nom de contrainte faible dans la terminologie de Sasaki (1970). Un autre type d'information peut être la vérification de certaines lois (hydrostatisme, géostrophisme). On parle alors de contrainte forte.

### C-Remarques

Ce formalisme peut s'appliquer au problème de l'initialisation des variables pronostiques d'un schéma de surface ($\bar{x}$) à partir d'observations atmosphériques au voisinage du sol ($\bar{y}_0$). On fait donc l'hypothèse implicite qu'il existe une corrélation assez forte entre ces quantités pour que notre problème puisse avoir une solution. Les équations du schéma de surface de Noilhan et Planton (1989) nous donnent une relation entre la température, l'humidité spécifique et le vent au premier niveau de calcul du modèle et les températures et contenus en eau des deux réservoirs du sol. Toute- fois, cette relation n'est ni linéaire, ni diagnostique, ni biunivoque. Les deux premiers inconvénients peuvent être résolus dans le cadre de l'assimilation variationnelle. Par contre l'ambiguïté venant des relations entre l'humidité superficielle (ou la résistance de surface) et l'état hydrique du sol a peu
de chances d'être résolue facilement. En effet, pour des contenus en eau compris entre la capacité au champ et la valeur à saturation (plage qui peut être très grande comme pour le sable, voir le tableau 1), l'évaporation s'effectuera toujours au taux potentiel. Nous verrons les conséquences de cela dans la présentation des résultats.
### 3-Présentation des simulations de référence

Nous présentons dans cette partie un ensemble de six simulations réalisées au site de Bats (SAMER 7) pour la journée du 19 mai 1986. Il s'agit d'un sol nu limoneux et très humide (la phase de désaturation des sols est à peine amorcée). Ces simulations ont été effectuées en utilisant une version unidimensionnelle du modèle PERIDOT avec le schéma de surface de Noilhan et Planton (1989). Les termes de forçage de grande échelle (advections et vitesse verticale) proviennent d'analyses réalisées toutes les 6h avec le modèle PERIDOT opérationnel pour l'ensemble de la Période d'Observations Intensives de la campagne HAPEX. La simulation 1 correspond aux conditions initiales fournies dans le tableau 2 issues des données de la campagne HAPEX (notamment pour l'eau du sol). Les autres simulations n'ont pour différences que le contenu en eau initial des deux réservoirs (Tableau 3). On considère la simulation 1 comme représentant la réalité et on se propose, partant des conditions initiales des simulations 2 à 6, de retrouver le contenu en eau du sol de cette simulation. Ainsi, pour chacune des méthodes qui vont être introduites on se place dans des conditions de fonctionnement "optimal" dans le sens où toute erreur sur les paramètres atmosphériques pour les simulations 2 à 6 vient, par construction, d'une erreur sur l'eau du sol.

La figure 1 rassemble l'évolution diurne du flux de chaleur latente pour les six simulations, ainsi que la mesure fournie par la station SAMER. La figure 2 représente l'évolution du réservoir superficiel pour chacune de ces simulations. L'examen global de la figure 1 montre que les conditions initiales pour l'eau du sol correspondent à une gamme de réponse pour l'évaporation assez étendue. Dans la simulation 1, le flux de chaleur latente est en bon accord avec la mesure sauf dans les phases nocturnes où le flux de rosée n'est pas reproduit par le modèle. Compte-tenu des conditions de surface, l'énergie radiative incidente est convertie principalement en flux de conduction et en flux d'évaporation. La simulation 2 conduit à une évolution du flux de chaleur latente très proche de la simulation de référence. L'évaporation à l'interface est directement pilotée par le contenu en eau du réservoir superficiel. Comme dans la simulation 1 ce réservoir était très proche de la capacité au champ (voir le tableau 1 qui regroupe les valeurs critiques de l'eau du sol pour l'ensemble des classes texturales), le fait d'y rajouter de l'eau modifie peu l'humidité relative toujours proche de 100 %. Dans la simulation 3, le réservoir superficiel est proche de la valeur au point de flétrissement,
toutefois, comme le réservoir profond est encore très humide, il peut réalimenter la surface évitant
ainsi un assèchement trop important, notamment en début de simulation jusqu'au lever du soleil
(Figure 2c). Bien que l'évaporation soit réduite par rapport à la simulation 1, elle a des valeurs très
importantes jusqu'à 12h. Le réservoir superficiel est plus humide dans la simulation 4 que dans la
simulation de référence. Les valeurs plus fortes du flux en milieu de journée s'expliquent par une
augmentation de l'humidité spécifique saturante en surface due à un échauffement plus important
de la surface (influence de l'eau du sol sur l'inertie thermique discutée plus loin). Le manque d'eau
dans le réservoir profond se fait sentir à partir de 14h, le réservoir superficiel s'étant fortement
vidé sans pouvoir être réalimenté par la couche profonde. Cela est nettement visible sur l'évolution
diurne de W, avec une décroissance régulière depuis le début de la simulation jusqu'à 17h (Figure
2d). C'est également ce qui se passe pour les simulations 5 et 6, où la diminution de l'évaporation
superficielle s'opère d'autant plus tôt que le réservoir superficiel contient peu d'eau.
# 4-Assimilation variationnelle de l'eau du sol

## A-Présentation générale de l'assimilation quadridimensionnelle

Dans l'assimilation 4D, les équations pronostiques du modèle à initialiser sont supposées être suffisamment exactes pour les imposer comme une contrainte forte au vecteur 4D: $\bar{x}$. On augmente ainsi $N_y$ en utilisant les observations non plus à un instant donné, mais sur une certaine période de temps. Le Dimet et Talagrand (1986) ont proposé de résoudre ce problème par réduction de la variable de contrôle. En effet, si certains éléments ($\bar{v}$) de $\bar{x}$ peuvent être calculés à partir d'autres ($\bar{w}$):

$$ \bar{x} = (\bar{w}^*, \bar{v}^*) $$

$$ \bar{v} = G(\bar{w}) $$

alors, le problème d'analyse sous contraintes pour $\bar{x}$ se transforme en un problème sans contrainte pour $\bar{w}$ (dont la résolution pratique est beaucoup plus simple).

Dans un problème d'assimilation, l'opérateur $G$ représente le modèle de prévision et $\bar{w}$ les conditions initiales. Ainsi, le problème initialement quadridimensionnel se réduit à la détermination du champ tridimensionnel des conditions initiales.

Il s'agit alors de minimiser:

$$ J(\bar{w}) = - \ln[P(\bar{y}_0 | \bar{w} = \bar{w}_t)] - \ln[P(\bar{w} | \bar{w}_b)] $$

Dans la pratique, cette quantité est minimisée par une méthode itérative nécessitant au minimum le calcul du gradient $J'$ de $J$ (voir Navon et Legler (1987) qui font une revue des méthodes numériques couramment utilisées en météorologie). Posons:

$$ d\bar{y} = K_n\{G(\bar{w})\} - \bar{y}_0 $$

$$ J_{of}(d\bar{y}) = - \ln[P(\bar{y}_0 | \bar{w} = \bar{w}_t)] $$

$$ J_b(\bar{w} - \bar{w}_b) = - \ln[P(\bar{w} | \bar{w}_b)] $$

On écrit alors la fonctionnelle et son gradient:

$$ J(\bar{w}) = J_{of}(d\bar{y}) + J_b(\bar{w} - \bar{w}_b) $$
$$J'(\bar{w}) = G^* K^* J'(d\bar{y}) + J_b'(\bar{w} - \bar{w}_b)$$

On constate que le calcul du gradient de la fonctionnelle fait apparaître les adjoints du modèle linéaire tangent ($G^*$) et de l'opérateur linéaire tangent de $K_n$ ($K^*$).

Le Dimet et Talagrand (1986) ont montré qu'une estimation de $J'$ peut s'obtenir par une intégration directe du modèle, suivie d'une intégration à rebours du modèle linéaire tangent adjoint. Le coût de cette méthode est donc d'environ deux intégrations du modèle de prévision sur la période d'assimilation par itération de l'algorithme de descente.

## B-Assimilation de l'état hydrique des sols

Ce qui vient d'être présenté sur l'assimilation peut s'appliquer à l'estimation de l'état initial des variables du sol.

Dans un premier temps, on peut privilégier les échanges verticaux par rapport aux échanges horizontaux, ce qui permet de se ramener à un problème unidimensionnel. On formule alors notre problème comme suit:

Connaissant des observations atmosphériques $\bar{y}_0$ au voisinage du sol ($HU_{2m}$) sur un certain intervalle de temps $[-T, 0]$, on cherche l'état initial $\bar{w}$ à $t = -T$ pour les variables pronostiques du schéma de surface ($W_s, W_p$) qui conduit à une évolution du modèle passant au plus près des observations tout en satisfaisant la contrainte forte des équations pronostiques.

Etant donné que le nombre de variables de contrôle est faible et que le modèle de prévision réduit à sa version 1D est très peu coûteux en temps-calcul, on peut, au moins pour une première étude, se dispenser de l'écriture de l'adjoint du modèle de sol. En effet, il est possible d'obtenir le gradient de la fonctionnelle à minimiser en perturbant l'état initial au moyen de $M$ intégrations du modèle, où $M$ est le nombre de variables de contrôle. Cette méthode, déjà utilisée par Hoffman (1986), est ici reprise.

Considérons un ensemble de $N$ observations $\bar{y}$ sur une période d'assimilation $[-T, 0]$. Soit $K_n(\bar{x})$ l'estimation de $\bar{y}$ faite par le modèle. On peut écrire cette quantité comme une fonction de l'état initial $\bar{w}$ à $t = -T$ soit: $\bar{y}_n(\bar{w})$. Hoffman (1986) a défini la valeur optimale de $\bar{w}$ comme étant
celle qui minimise l'erreur quadratique totale:

$$S(\vec{w}) = \sum_{n} [f_{n}(\vec{w})]^{2}$$

où

$$f_n(\vec{w}) = \alpha_n [\vec{y}_n(\vec{w}) - \vec{y}]$$

est le résidu pondéré de la n-ième observation.

La fonctionnelle $S(\vec{w})$ apparaît comme un cas particulier de la fonctionnelle $J(\vec{w})$ dans laquelle on a négligé l'ajustement à une estimation préalable, et où l'on a considéré que la matrice de covariance des erreurs d'observations était diagonale. C'est une fonctionnelle de cette forme qui a été définie par Courtier (1987) dans l'étude de l'assimilation du géopotentiel à 500 hPa.

Les poids $\alpha_n$ doivent inclure les erreurs de mesure et dépendre de l'âge des données pour rendre compte des imperfections du modèle.

L'utilisation d'une fonctionnelle quadratique $S(\vec{w})$ permet de mettre en oeuvre un algorithme d'optimisation très simple utilisant la méthode de Gauss-Newton.

Soit $\vec{w}^{(k)}$ une estimation de l'état initial à la k-ième itération du processus d'optimisation. On obtient $\vec{w}^{(k+1)}$ en linéarisant $f_n$ au voisinage de $\vec{w}^{(k)}$. On considère une approximation de la méthode de Newton. La méthode de Newton cherche $\vec{w}$ tel que le gradient de $S$ s'annule, soit:

$$\frac{1}{2} \frac{\partial S}{\partial w_i} = \sum_n f_n(\vec{w}) \frac{\partial f_n(\vec{w})}{\partial w_i} = 0$$

En développant $f_n$ et $\partial f_n / \partial w_i$ en série de Taylor, et en négligeant les termes d'ordre élevés, on obtient:

$$\sum_n \sum_j \left[ \frac{\partial f_n}{\partial w_i} \frac{\partial f_n}{\partial w_j} + f_n \frac{\partial^2 f_n}{\partial w_i \partial w_j} \right] P_j = - \sum_n f_n \frac{\partial f_n}{\partial w_i}$$

où

$$P_j = w_j^{(k+1)} - w_j^{(k)}$$

$f_n$ et ses dérivées sont évalués au point $\vec{w}^{(k)}$.

La méthode de Gauss-Newton consiste à négliger les dérivées secondes de $f_n$ (ce qui peut être une mauvaise approximation en cas de problème fortement non linéaire). On a alors à résoudre le
système:

$$J.\vec{P} = -\vec{f}$$

où $\vec{f}$ est le vecteur de résidus et J la matrice jacobienne $J_{ni} = \partial f_n / \partial w_i$.

Ce système est surdéterminé dès que le nombre d'observations est supérieur aux nombre de variables de contrôle. On résout ce système en décomposant J sous la forme QR où R est une matrice triangulaire supérieure et Q une matrice orthogonale. La méthode de factorisation utilise les matrices de Householder (Ciarlet, 1982) conduisant au nouveau système:

$$R.\vec{P} = -Q^T \vec{f}$$

que l'on peut facilement résoudre par une méthode de remontée.

Dans la pratique, la matrice jacobienne est approximée en différences finies décentrées:

$$J_{ni}(\vec{w}) = \frac{f_n(\vec{w} + h_i \vec{u}_i) - f_n(\vec{w})}{h_i}$$

où $\vec{u}_i$ est le n-ième vecteur unitaire et $h_i$ un pas approprié.

### C- Résultats

Nous considérons ici que la variable de contrôle est le couple $(W_s, W_p)$ et que la fonctionnelle à minimiser représente l'ajustement à l'humidité à 2m issue de la simulation 1 échantillonnée toutes les heures.

Nous venons de voir que les éléments de la matrice jacobienne $J_{ni}(\vec{w})$ se calculaient en différences finies. Ainsi:

$$\frac{\partial f_n}{\partial W_s} = \frac{f_n(W_s + \delta W_s, W_p) - f_n(W_s, W_p)}{\delta W_s}$$

On écrit l'incrément $\delta W_s$ sous la forme $\varepsilon W_s$ où $\varepsilon$ est un nombre arbitrairement petit. Une condition nécessaire pour que la méthode décrite précédemment soit utilisable est que les quantités $\partial f_n / \partial W_s$, aient un sens, à savoir que:

$$\lim_{\epsilon \to 0} \frac{f_n(W_s(1+\epsilon), W_p) - f_n(W_s, W_p)}{\epsilon W_s}$$

existe.
Partant des conditions initiales de la simulation 3 ($W_s = 1.7\,\text{mm}, W_p = 320\,\text{mm}$) nous avons calculé les éléments de la matrice J lors de la première itération du processus de minimisation en faisant varier $\log \varepsilon$ depuis -1 jusqu'à -10. La figure 3 représente l'évolution de ces éléments. Lorsque $\varepsilon$ est trop grand, le calcul des dérivées est trop imprécis. De même, lorsque $\varepsilon$ devient trop faible, des problèmes liés à la précision des calculs apparaissent. Dans une gamme comprise entre $10^{-4}$ et $10^{-7}$, le calcul est stable assurant la différentiabilité des $f_n$. Les résultats présentés par la suite ont été obtenus avec $\varepsilon = 10^{-6}$.

On peut noter les valeurs importantes des dérivées par rapport à $W_s$, alors que les dérivées par rapport à $W_p$ sont beaucoup plus faibles. L'algorithme va donc modifier principalement l'eau du réservoir superficiel, ce qui est conforme à ce que l'on attend de lui pour retrouver l'eau de la simulation 1.

Pour cette initialisation de l'eau du sol, nous avons reporté la valeur du contenu en eau du sol à chaque itération du processus de descente (Figure 4). On constate une convergence très rapide de la méthode. Après 7 itérations on retrouve précisément les valeurs initiales de la simulation 1. L'évolution de la fonctionnelle montre qu'en 4 itérations la valeur initiale a été réduite d'un facteur supérieur à $3.10^4$ (Figure 5). Elle tend très rapidement vers zéro qui est par construction son minimum.

Le tableau 4 regroupe l'ensemble des résultats de la méthode variationnelle pour les initialisations de l'eau du sol correspondant aux simulations 2 à 6.

Partant de $W_s = 3\,\text{mm}$ et $W_p = 320\,\text{mm}$, on a initialement une valeur de la fonctionnelle $S$ très faible (correspondant à un écart-type d'erreur d'humidité relative $\sigma_{HU} = \sqrt{S/N}$ de 0.2%). On est dans une situation où la réponse atmosphérique est très voisine de la situation de référence (Figures 1a et 1b). Il semble que l'algorithme n'est pas capable de lever cette ambiguïté. Dans ce cas, la meilleure estimation de l'eau du sol correspond aux valeurs initiales introduites dans le processus de minimisation. Il ne semble pas gênant de conserver ces valeurs puisque ce qui nous intéresse en définitive c'est une bonne prévision des paramètres atmosphériques près de la surface. Toutefois, on peut penser qu'en prenant une période d'assimilation plus longue on puisse arriver à lever l'ambiguïté. Pour la simulation 3, comme nous l'avons vu précédemment, le processus est très
efficace car alors les évolutions diurnes de l'humidité relative sont très contrastées. Dans le cas de la simulation 4, on réduit fortement la fonction-coût en deux itérations (rapport 2000) en ajoutant de l'eau au réservoir profond. Toutefois, l'algorithme perd ensuite son efficacité puisque on se trouve alors dans une zone de l'espace des paramètres proche de celle de la simulation 2. Les trop faibles gradients d'humidité atmosphérique par rapport à l'humidité du sol font diverger la méthode à l'itération suivante. Les mêmes conclusions sont applicables pour la simulation 5. La simulation 6, correspondant à un état initial très loin de l'état cherché conduit à une convergence rapide du fait que l'on se situe dans une zone de l'espace des paramètres où les gradients sont élevés.

De cette première étude, il ressort que l'assimilation de l'eau du sol à partir de l'observation de l'humidité à 2m sur une période de 24h donne ici des résultats très encourageants. Un résultat nouveau par rapport aux assimilations "dynamiques" vient de ce que l'on n'a pas de relation biunivoque entre la quantité que l'on cherche à estimer et l'observation. La conséquence est que l'algorithme de descente peut conduire à des valeurs de l'eau du sol sensiblement différentes des valeurs "vraies" mais qui correspondent à un accord avec les observations tout à fait satisfaisant (comme on peut le voir sur les valeurs de $σ_{HU}$ final reportées dans le tableau 4).

Signalons pour terminer un problème lié au fait que les variables de contrôle de notre problème sont des quantités bornées ($0 \le W \le W_{sat}$). L'algorithme peut générer des valeurs en dehors de ce domaine de définition. Dans de telles situations nous avons ramené les points extérieurs dans le domaine de façon arbitraire. Il serait préférable de rajouter une pénalité à la fonction-coût comme cela a été fait par Lorenc (1988) et suggéré par Courtier (1987). Toutefois, la fonctionnelle à minimiser n'est alors plus quadratique, ce qui oblige à utiliser un algorithme de descente plus complexe que celui retenu pour la présente étude en faisant appel à des bibliothèques informatiques (IMSL, NAG).
# 5-Assimilation séquentielle de l'eau du sol

L'assimilation séquentielle (ou cycle d'assimilation) est constituée par l'ensemble d'opérations suivantes:

a) Prévision à 6h d'échéance

b) Analyse à l'échéance 6h consistant en une combinaison linéaire optimale de la prévision et des observations disponibles

c) Retour en (a) avec comme conditions initiales l'analyse obtenue en (b)

## A-Méthode de Geleyn

La méthode des incrémentes proposée par Geleyn et utilisée opérationnellement est décrite en annexe.

C'est l'adaptation de cette méthode au nouveau schéma de surface et réduite à l'initialisation de l'eau du sol que nous décrivons maintenant.

On se donne l'incrément d'humidité relative au premier niveau de calcul du modèle entre l'analyse (a) et la prévision (p), soit:

$$ \Delta HU = HU^a - HUP^p $$

On fait l'hypothèse que cette erreur de prévision dans l'atmosphère provient d'une erreur de prévision de l'eau du sol par l'intermédiaire de l'humidité relative superficielle. Ainsi, si $M^p$ est l'humidité relative superficielle prévue, la valeur analysée s'obtient en reportant l'incrément atmosphérique en surface soit:

$$ M^a - M^p = HU^a - HUP^p $$

$$ M^a = M^p + \Delta HU $$

Comme l'humidité relative en surface est reliée à l'eau du sol superficiel par la relation suivante:

$$ M = \frac{1}{2} \left[ 1 - \cos \left( \frac{W_s \pi}{W_{fc}} \right) \right] $$

on peut alors fournir une nouvelle estimation de l'eau du réservoir superficiel:

$$ W_s^a = \frac{W_{fc}}{\pi} \arccos(1 - 2M^a) \quad \text{si} \quad 0 \le M^a \le 1 $$
Si l'estimation conduit à une valeur de $M^a$ supérieure à un, on majore $W_s^a$ par la valeur correspondant à la capacité au champ ($W_{fc}$). On applique une correction identique sur le réservoir profond soit:

$$W_p^a = W_p^p + (W_s^a - W_p^p)$$

Notons que cette adaptation de la méthode de Geleyn au nouveau schéma de surface n'est possible que sur sol nu. En effet, on montre facilement que sur une végétation partiellement couvrante l'humidité de surface s'exprime en fonction de l'eau des deux réservoirs.

Nous avons testé cette méthode sur un cycle d'assimilation de 24h. Les analyses atmosphériques introduites toutes les 6h correspondent aux valeurs calculées de l'humidité relative au premier niveau du modèle dans la simulation 1. L'efficacité de la méthode va se juger à son aptitude à reproduire les contenus en eau du sol obtenus après 24h de simulation dans la simulation 1. Ces valeurs sont reportées pour l'ensemble des simulations dans le tableau 5. Lorsque l'eau du réservoir profond est initialement correcte, la méthode corrige dans le bon sens l'eau du sol superficiel. On peut toutefois émettre une petite réserve venant du fait que la valeur à retrouver est proche de la capacité au champ (2.3mm pour 2.4mm). Lorsque l'humidité relative analysée vaut 1, l'eau du sol peut être comprise entre 2.4mm et 4.5mm et rien ne nous permet de choisir une valeur plutôt qu'une autre. En choisissant d'imposer la valeur à la capacité au champ on biaise quelque peu la méthode en notre faveur. Si le contenu en eau superficiel était voisin de la saturation, la correction imposerait la valeur à la capacité au champ. Toutefois, on peut penser que la prévision qui en résulterait conduirait à une atmosphère trop sèche d'où une réhumidification du sol au prochain cycle. Par contre, lorsque le sol profond est trop sec, les quantités d'eau réinjectées à chaque cycle sont trop négligeables pour qu'une correction sensible puisse être notée en 24h. Ce qui pouvait marcher avec le schéma opérationnel où le rapport entre les épaisseurs des deux réservoirs n'est que de 5, se révèle totalement inadapté avec le schéma de Noilhan et Planton où l'on a un rapport 100 entre ces quantités. En conclusion de ce test, on retrouve avec la méthode de Geleyn, l'ambiguïté sur l'eau du sol superficiel venant de la relation entre $M$ et $W_s$ déjà évoquée dans l'assimilation variationnelle. Il semble de plus que cette méthode ne soit pas capable de fournir une information utile sur le
contenu en eau du sol profond. Etant donnée la disparité entre les épaisseurs de deux réservoirs, il faudrait des temps d'assimilation extrêmement longs pour que les corrections appliquées au réservoir superficiel puissent avoir une influence notable sur l'eau du sol profond.

## B-Méthode de Geleyn modifiée

On se propose ici de savoir quelle est l'exactitude de la relation reliant les incréments atmosphériques aux incréments en surface:

$$M^a - M^p = HU^a - HU^p$$

Pour les simulations 2 et 3, on a comparé l'incrément d'humidité que l'on ajoute effectivement au niveau de la surface ($\Delta HU$) avec celui que l'on devrait ajouter ($\Delta M$) pour retrouver l'eau du sol de la simulation 1. On constate (Fig. 6) que ces incréments ont une évolution diurne assez bien correlée, toutefois, l'information atmosphérique conduit à sous estimer le déficit en humidité relative de la surface en milieu de journée (soit pour les analyses de 12h et 18h).

On peut alors améliorer la formulation initiale en effectuant une régression linéaire entre les incréments atmosphériques et en surface:

$$M^a = M^p + \alpha (HU^a - HU^p)$$

Supposons que $HU^a$ et $M^a$ représentent la valeur vraie pour ces quantités (pas d'erreur d'analyse), en notant les erreurs de prévision $\epsilon M = M^p - M^a$ et $\epsilon HU = HU^p - HU^a$, il vient pour $\alpha$:

$$\alpha = \frac{\overline{\epsilon M \epsilon HU}}{\overline{\epsilon HU^2}}$$

On considère que la valeur vraie est celle issue de la simulation 1. En perturbant cette expérience par modification du contenu en eau initial des sols, on peut calculer les corrélations $\overline{\epsilon M \epsilon HU}$ et $\overline{\epsilon HU^2}$. Soit:

$$\overline{\epsilon M \epsilon HU} = \frac{1}{N} \sum_{i=1}^{N} (M^p - M^a)_i (HU^p - HU^a)_i$$

$$\overline{\epsilon HU^2} = \frac{1}{N} \sum_{i=1}^{N} (HU^p - HU^a)^2$$
On a fait varier la contenu en eau volumique de 0.14 $m^3 m^{-3}$ à 0.32 $m^3 m^{-3}$ par pas de 0.02 pour les deux réservoirs. Cela correspond à un ensemble de 100 simulations perturbées (150s de CRAY 2) permettant de calculer les variances et covariances aux heures d'analyse (6h, 12h, 18h, 24h). Les valeurs du paramètre $\alpha$ reportées dans le tableau 6 montrent que l'on devrait augmenter de façon significative l'incrément atmosphérique au niveau de la surface en milieu de journée. Le résultat des assimilations réalisées avec cette méthode modifiée (Tableau 7) révèlent une meilleure correction du réservoir superficiel. L'eau du sol profond n'est bien sûr que très légèrement modifiée. Pour la simulation 5 nous avons comparé l'évolution de $W_s$ et $W_p$ au cours de l'assimilation avec la méthode de Geleyn (Figure 7a) et avec la méthode de Geleyn modifiée (Figure 7b). C'est principalement à 18h que l'augmentation de l'incrément conduit à une valeur plus proche de la quantité cherchée. On constate dans cette situation particulière (réservoir superficiel correct, réservoir profond trop sec) que la surface s'assèche rapidement entre 2 analyses à cause d'un défaut en profondeur, d'où un manque d'efficacité de la correction appliquée toutes les 6h.

### C-Interpolation optimale

L'analyse optimale consiste à estimer un champ météorologique comme une combinaison linéaire d'autres champs. Cette combinaison est optimale dans le sens où elle minimise la variance de l'erreur d'analyse. Si $x$ est un champ à estimer (prédicand) et $y_i$ un ensemble de prédicteurs, on préfère effectuer l'analyse sur les résidus "observation-guess". On écrit alors:

$$x^a = x^p + \sum_{i=1}^{n} \lambda_i (y_i^o - y_i^p)$$

Cette expression est équivalente à l'équation (3) introduite dans la présentation générale du chapitre 2. Les indices $a, o$ et $p$ représentent respectivement l'analyse, l'observation et la prévision. Minimiser $\overline{\varepsilon x^a}$ conduit à l'estimation des coefficients $\lambda_i$ qui sont solution du système linéaire suivant:

$$\sum_{i=1}^{n} \lambda_i (\overline{\varepsilon y_i^o \varepsilon y_j^o} + \overline{\varepsilon y_i^p \varepsilon y_j^p}) = \overline{\varepsilon x^p \varepsilon y_j^p}$$

Cette présentation nous montre que la méthode de Geleyn est une interpolation optimale dégradée avec comme prédicteur $HU$ et comme prédicand $M$, et un coefficient $\lambda = 1$ non optimal. La méthode modifiée consiste à choisir $\lambda$ de façon optimale. Les conclusions que nous avons obtenues
avec la méthode de Geleyn sont que pour l'eau superficielle on converge assez bien vers la valeur cherchée avec toujours une ambiguïté entre la capacité au champ et la saturation, mais par contre elle ne peut fournir aucune information concernant l'eau du sol profond. Pour pallier ces problèmes on peut essayer de pousser plus loin l'analogie avec l'interpolation optimale et considérer que les prédictands sont les contenus en eau des deux réservoirs $W_s$ et $W_p$ (évitant ainsi l'inversion de la relation $M = M(W_s)$) et comme prédicteurs les deux quantités au voisinage du sol avec qui ils semblent à priori bien corrélés: l'humidité relative $HU$ et la température $T$. On écrit alors:

$$W_s^a = W_s^p + \lambda_s^1(T^o - T^p) + \lambda_s^2(HU^o - HU^p)$$

$$W_p^a = W_p^p + \lambda_p^1(T^o - T^p) + \lambda_p^2(HU^o - HU^p)$$

Les coefficients $\lambda$ sont solution des systèmes:

$$\begin{bmatrix} \frac{\overline{\varepsilon T^p + \varepsilon T^o}}{\overline{\varepsilon T^p \varepsilon H U^p}} & \frac{\overline{\varepsilon T^p \varepsilon H U^p}}{\overline{\varepsilon H U^p + \varepsilon H U^o}} \\ \frac{\overline{\varepsilon T^p + \varepsilon T^o}}{\overline{\varepsilon T^p \varepsilon H U^p}} & \frac{\overline{\varepsilon T^p \varepsilon H U^p}}{\overline{\varepsilon H U^p + \varepsilon H U^o}} \end{bmatrix} \begin{bmatrix} \lambda_s^1 \\ \lambda_s^2 \end{bmatrix} = \begin{bmatrix} \frac{\overline{\varepsilon W_s^p \varepsilon T^p}}{\overline{\varepsilon W_s^p \varepsilon H U^p}} \\ \frac{\overline{\varepsilon W_p^p \varepsilon T^p}}{\overline{\varepsilon W_p^p \varepsilon H U^p}} \end{bmatrix}$$

Les termes $\varepsilon X^p$ et $\varepsilon X^o$ représentent les erreurs du guess et de l'observation respectivement. Dans notre étude nous négligerons les erreurs d'observation ($\overline{\varepsilon X^o} = 0$). La valeur vraie d'un paramètre sera considérée comme celle issue de la simulation 1. En perturbant cette expérience par modification du contenu en eau des sols, on peut calculer les corrélations $\overline{\varepsilon X^p \varepsilon Y^p}$ et $\overline{\varepsilon Y^p}$. Nous avons effectué un ensemble de 100 simulations en initialisant le contenu en eau des sols de la façon suivante:

$$W_s = W_s^* \left( \zeta + \frac{1}{2} \right)$$

$$W_p = W_p^* \left( \eta + \frac{1}{2} \right)$$

$W_s^* = 2.3 \, mm$ et $W_p^* = 320 \, mm$. $\zeta$ et $\eta$ sont deux nombres aléatoires compris entre 0 et 1. Ce choix à été fait pour que $\overline{\varepsilon W_i} = \overline{W_i} - W_i^* = 0$ puisque $\overline{\zeta} = \overline{\eta} = 0.5$.

Les coefficients $\lambda$ issus de cet ensemble de simulations sont rassemblés dans le tableau 8 aux heures d'analyse. Dans le tableau 9, ont été reportées les valeurs des coefficients de corrélation multiple toutes les 6 heures. Ils sont donnés par la formule:

$$r_{ijk} = \sqrt{\frac{\rho_{ij}^2 + \rho_{ik}^2 - 2\rho_{ij}\rho_{ik}\rho_{jk}}{1 - \rho_{jk}^2}}$$
Les quantités $ρ_{ij}$ représentent les coefficients de corrélation d'ordre 0 entre les variables i et j (tableau 9). L'indice i prend les valeurs $W_s$ et $W_p$, alors que les indices j et k prennent les valeurs T et HU.

On constate que l'hypothèse de linéarité est assez bonne puisqu'à l'exception des coefficients de corrélation incluant $W_p$ à 18h on a des valeurs comprises entre 0.80 et 0.99. On peut donc espérer obtenir une information utile par une simple interpolation optimale sur les contenus en eau du sol.

Le tableau 10 rassemble les valeurs des contenus en eau obtenus après 24h d'assimilation par interpolation optimale. L'apport principal par rapport à la méthode de Geleyn est d'être capable de réalimenter le sol profond en eau lorsque celui-ci est initialement trop sec. Dans le cas de la simulation 2, où la réponse atmosphérique est pratiquement la même que pour la simulation de référence, la méthode n'introduit aucune modification au niveau de l'eau du sol. Sur la figure 8 sont rassemblées les évolutions temporelles des contenus en eau pour les 5 assimilations et pour la simulation 1. La correction la plus importante pour le sol profond apparaît à 6h où l'on injecte plus de 100 mm d'eau lorsqu'on initialise avec $W_p = 200\text{mm}$. Les corrections effectuées par la suite sont toujours beaucoup plus faibles. Sur la figure 8a, on constate que seules des corrections infimes ont été appliquées à l'eau du sol. La figure 8b montre que l'on applique une correction importante au réservoir superficiel à 18h. Par contre on enlève 20mm à 12h au réservoir profond alors que la prévision était correcte. Comme il s'agit d'une méthode statistique, on est pas capable d'appliquer la bonne correction systématiquement. Globalement les corrections sont satisfaisantes puisqu'on humidifie la surface sans modifier de façon importante l'eau du sol profond. Sur la figure 8c, un défaut du même type apparaît. A 6h les incréments atmosphériques indiquent qu'il faut injecter de l'eau dans le sol. C'est bien ce que fait l'algorithme, mais il ajoute également de l'eau dans le réservoir superficiel qui était déjà trop humide. Toutefois on constate que la correction importante qu'il fallait effectuer sur le réservoir profond a bien été prise en compte. L'information provient de la liaison entre le coefficient d'inertie thermique $C_G$ et l'eau du sol profond (relation 10 dans Noilhan et Planton, 1989). Ainsi lorsque le réservoir profond est trop sec, la température atmosphérique est plus froide que dans la simulation de référence. On voit ici l'intérêt d'avoir considéré la température comme prédicteur, ainsi que la supériorité du nouveau schéma puisque
dans le schéma opérationnel (annexe) $C_{sol}$ est un coefficient constant. A l'analyse de 12h, on corrige quelque peu la réhumidification trop importante des deux réservoirs faite à 6h. Un comportement identique est observé sur la figure 8d, avec pour seule différence le fait que l'humidification du réservoir superficiel est dans ce cas justifiée. Finalement, en partant d'un sol trop sec dans les deux couches, on corrige très bien l'eau du sol (Figure 8e). On constate qu'il suffit d'une seule analyse (6h) pour retrouver les valeurs de la simulation de référence.

En conclusion, l'interpolation optimale est plus satisfaisante que la méthode de Geleyn dans la mesure où elle permet de corriger l'eau du sol profond quand c'est nécessaire. Il n'y a bien évidemment pas de miracle, cela provient de l'information supplémentaire injectée sous forme des covariances entre les variables observées et estimées, ainsi que de l'utilisation de la température comme prédicteur.

Lorsque les variables atmosphériques sont très différentes de celles de la simulation de référence, comme pour l'assimilation 4D on retrouve facilement les valeurs cherchées. Lorsque l'eau du sol ne modifie pas les variables atmosphériques, on garde les valeurs prévues. Dans les situations intermédiaires, les corrections sont globalement satisfaisantes, avec quelques erreurs liées à la nature statistique de la méthode et à l'hypothèse de linéarité.
## 6-Conclusion

Il s'agissait dans ce travail d'effectuer une étude préliminaire sur la faisabilité de méthodes d'analyse déjà existantes appliquées à l'initialisation de l'eau du sol dans les modèles de prévision numérique.

L'assimilation variationnelle est apparue très satisfaisante dans la mesure où elle prend en compte de façon explicite la liaison temporelle non linéaire $K_n$ entre les observations et les paramètres à analyser. C'est avec cette méthode que l'on a utilisé le maximum d'informations: les équations du modèle de prévision et des observations horaires. Il faudra regarder à l'avenir les modifications apportées par une cadence temporelle moins rapide. Il semble aussi nécessaire d'avoir un algorithme de minimisation plus performant, ne serait ce que pour prendre en compte le caractère borné des variables à estimer. On cherchera également à comprendre les causes de la mauvaise convergence de l'algorithme pour certaines valeurs initiales (mauvais calcul du gradient, existence de plusieurs bassins d'attraction,...).

La méthode des incrémentes proposée par Geleyn corrige raisonnablement l'eau du sol superfi ciel, mais est inefficace pour retrouver l'eau du sol profond.

La méthode d'interpolation optimale, bien qu'elle ne considère que des liaisons linéaires entre les observations et l'état à estimer, a donné des résultats encourageants. Pour utiliser cette méthode de façon plus générale, il faudra modéliser les matrices de variances-covariances en fonction de la couverture végétale, du type de sol et de l'heure locale d'analyse.

Cette première étude va être poursuivie en utilisant pleinement les données de la campagne HAPEX. Nous réaliserons tout d'abord des simulations de référence de longue durée. Elles serviront d'une part à vérifier la gestion du bilan hydrique de surface par le schéma de Noilhan et Planton, et d'autre part à fournir des valeurs d'eau du sol entre les mesures hebdomadaires effectuées par sondages neutroniques. Sur des périodes montrant une bonne cohérence entre l'observation et la simulation dans l'atmosphère et dans le sol, nous testerons les méthodes introduites dans cette note.
# ANNEXE

## Equations pronostiques du schéma de surface opérationnel

$$ \frac{\partial T_s}{\partial t} = C_{sol}H_{at} - \frac{2\pi}{\tau_1}(T_s - T_p) \quad \tau_1 = 24h $$

$$ \frac{\partial T_p}{\partial t} = \frac{2\pi}{\tau_2}(T_s - T_p) \quad \tau_2 = 5\text{jours} $$

$$ \frac{\partial W_s}{\partial t} = P - J_q - \frac{W_s - W_p(W_{smax}/W_{pmax})}{\tau_1} \quad W_{smax} = 20\text{mm} $$

$$ \frac{\partial W_p}{\partial t} = \frac{W_s - W_p(W_{smax}/W_{pmax})}{\tau_1} \quad W_{pmax} = 100\text{mm} $$

avec

$$ H_{at} = R_N - (J_s - (C_{pv} - C_{pa})T J_q + LJ_q) $$

$$ J_s = -\rho C_H V_{NK} (C_p T_{NK} + g Z_{NK} - C_p T_s) $$

$$ J_q = -\rho C_H V_{NK} (q_{NK} - H U_s q_{sat}(T_s)) $$

$$ HU_s = 0.5 \left[ 1 - \cos \left( \frac{\pi W_s}{W_{smax}} \right) \right] \quad (1) $$

## Méthode d'analyse par incréments

On se donne les incréments au premier niveau atmosphérique du modèle (NK):

$$ \Delta T = T_{NK}^{a} - T_{NK}^{p} $$

$$ \Delta HU = HU_{NK}^{a} - HU_{NK}^{p} $$

Les valeurs analysées en surface s'écrivent:

$$ T_s^a = T_s^p + \Delta T $$

$$ T_p^a = T_p^p + \Delta T \left( \frac{\tau_1}{\tau_2} \right) $$

$$ HU_s^a = HU_p^a + \Delta HU $$

On en déduit $W_s^a$ par la relation diagnostique (1), et on impose en outre

$$ W_p^a = W_p^p + \Delta W_s $$

avec

$$ \Delta W_s = W_s^a - W_p^a $$
REFERENCES

Abramopoulos, F., C. Rosenzweig and B. Choudhury, 1988 : Improved ground hydrology calculations for Global Climate Models (GCMs): Soil water movement and evapotranspiration. *J. Climate*, **1**, 921-941

Benjamin, S.G. and T.N. Carlson, 1986 : Some effects of surface heating and topography on the regional severe storm environment. Part I: 3-D simulations. *Mon. Wea. Rev.*, **114**, 307-329

Ciarlet, P.G., 1982 : Introduction à l'analyse numérique matricielle et à l'optimisation. Masson, Paris, 395pp.

Courtier, P., 1987 : Application du contrôle optimal à la prévision numérique en météorologie. Thèse de doctorat de l'Université de Paris VI.

Hoffman, R.N., 1986 : A four-dimensional analysis exactly satisfying equations of motion. *Mon. Wea. Rev*, **114**, 388-397.

Jacquemin, B., 1987 : Test des paramétrisations des processus de surface dans le modèle PERIDOT de mésoéchelle. Note de travail de l'ENM, 216pp.

Lanicci J.A., T.N. Carlson and T.T. Warner, 1987 : Sensitivity of great plains severe-storm environment to soil-moisture distribution. *Mon. Wea. Rev.*, **115**, 2660-2673

Le Dimet, F-X and O. Talagrand, 1986 : Variational algorithms for analysis and assimilation of meteorological observations: theoritical aspects. *Tellus*, **38A**, 97-110.

Lorenc, A.C., 1986 : Analysis methods for numerical weather prediction. *Q. R. J. Meteorol. Soc.*, **112**, 1177-1194.

Lorenc, A.C., 1988 : Optimal nonlinear objective analysis. *Q. R. J. Meteorol. Soc.*, **114**, 205-240.

McCumber, M.C., and R.A. Pielke, 1981 : Simulation of the effects of surface fluxes of heat and moisture in a mesoscale numerical model. Part I: Soil layer. *J. Geophys. Res.*, **86**, 9929-9938.
Mahfouf, J.F., 1986 : Contribution à la définition d'une paramétrisation des transferts entre le sol, la végétation et l'atmosphère. Analyse de sensibilité et insertion dans un modèle mésoéchelle. Thèse de doctorat de l'Université de Clermont II.

Navon I.M and D.M. Legler, 1987 : Conjugate-gradients for large-scale minimisation in meteorology
*Mon. Wea. Rev.*, **115**, 1479-1502

Noilhan J., and S. Planton, 1989 : A simple parameterization of land surface processes for meteorological models. *Mon. Wea. Rev.*, **117**, 536-549

Pinty J.P., P. Mascart, E. Richard and R. Rosset, 1989 : An investigation of mesoscale flows induced by vegetation inhomogeneities using an evapotranspiration model calibrated against HAPEX-MOBILHY data. *J. Appl. Meteor.*, **28**, 976-992

Sasaki, Y., 1970 : Some basic formalisms on numerical variational analysis. *Mon. Wea. Rev.*, **98**, 875-883

Segal M., R. Avissar, M.C. McCumber and R.A. Pielke, 1988 : Evaluation of vegetation effects on the generation and modification of mesoscale circulations. *J. Atmos. Sci.*, **45**, 2268-2292

Sellers P.J. and J.L. Dorman : Testing the Simple Biosphere Model (SiB) using point micrometeorological and biophysical data. *J. Climate Appl. Meteor.*, **26**, 622-651.

Taconet, O. R. Bernard and D. Vidal-Madjar, 1986 : Evapotranspiration over an agricultural region using surface flux/temperature model based on NOAA-AVHRR data. *J. Climate Appl. Meteor.*, **25**, 284-307.

Wetzel, P.J., D. Atlas and R.H. Woodward, 1984 : Determining soil moisture from geosynchronous satellite infrared data : A feasibility study. *J. Climate Appl. Meteor.*, **23**, 375-391.

Wilson M.F., A. Henderson-Sellers, R.E. Dickinson and P.J. Kennedy, 1987 : Sensitivity of the Biosphere-Atmosphere Transfert Scheme (BATS) to the inclusion of variable soil characteristics. *J. Climate Appl. Meteor.*, **26**, 341-362
<table><thead><tr><th>TYPE DE SOL</th><th>η<sub>s</sub> (m<sup>3</sup>/m<sup>3</sup>)</th><th>η<sub>fc</sub> (m<sup>3</sup>/m<sup>3</sup>)</th><th>η<sub>wilt</sub> (m<sup>3</sup>/m<sup>3</sup>)</th></tr></thead><tbody><tr><td>Sable</td><td>0.395</td><td>0.135</td><td>0.068</td></tr><tr><td>Sable limoneux</td><td>0.410</td><td>0.150</td><td>0.075</td></tr><tr><td>Limon sableux</td><td>0.435</td><td>0.195</td><td>0.114</td></tr><tr><td>Limon fin</td><td>0.485</td><td>0.255</td><td>0.179</td></tr><tr><td>Limon</td><td>0.451</td><td>0.240</td><td>0.155</td></tr><tr><td>Limon sablo-argileux</td><td>0.420</td><td>0.255</td><td>0.175</td></tr><tr><td>Limon fin argileux</td><td>0.477</td><td>0.322</td><td>0.218</td></tr><tr><td>Limon argileux</td><td>0.476</td><td>0.325</td><td>0.250</td></tr><tr><td>Argile sableuse</td><td>0.426</td><td>0.310</td><td>0.219</td></tr><tr><td>Argile limoneuse</td><td>0.482</td><td>0.370</td><td>0.283</td></tr><tr><td>Argile</td><td>0.482</td><td>0.367</td><td>0.286</td></tr></tbody></table>

**Tableau 1:** Valeurs critiques des contenus en eau des sols de la classification texturale de Clapp et Hornberger. La capacité au champ ($\eta_{fc}$) correspond à une conductivité hydraulique $K_{\eta}$ de 0.1mm/jour. Le point de flétrissement ($\eta_{wilt}$) est défini par un potentiel hydrique $\psi$ de -15bars.
<table><caption>**Tableau 2:** Conditions initiales de la simulation 1 (prise comme référence)</caption><tr><td>W<sub>s</sub> (mm)</td><td>2.3</td></tr><tr><td>W<sub>p</sub> (mm)</td><td>320.</td></tr><tr><td>T<sub>s</sub> (C)</td><td>18.</td></tr><tr><td>T<sub>p</sub> (C)</td><td>21.</td></tr><tr><td>veg (%)</td><td>0.0</td></tr><tr><td>Z<sub>0</sub> (m)</td><td>0.01</td></tr><tr><td>Type de sol</td><td>Limon</td></tr></table>

<table><caption>**Tableau 3:** Valeurs initiales de l'eau du sol des réservoirs superficiel (W<sub>s</sub>) et profond (W<sub>p</sub>) pour les simulations de 1 à 6</caption><thead><tr><td>Simulation</td><td>W<sub>s</sub> (mm)</td><td>W<sub>p</sub> (mm)</td></tr></thead><tbody><tr><td>1</td><td>2.3</td><td>320.</td></tr><tr><td>2</td><td>3.0</td><td>320.</td></tr><tr><td>3</td><td>1.7</td><td>320.</td></tr><tr><td>4</td><td>3.0</td><td>200.</td></tr><tr><td>5</td><td>2.3</td><td>200.</td></tr><tr><td>6</td><td>1.7</td><td>200.</td></tr></tbody></table>
<table><thead><tr><th></th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Simulation</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>W<sub>s</sub> initial</td><td>3.0</td><td>1.7</td><td>3.0</td><td>2.3</td><td>1.7</td></tr><tr><td>W<sub>p</sub> initial</td><td>320.</td><td>320.</td><td>200.</td><td>200.</td><td>200.</td></tr><tr><td>S initiale</td><td>0.931</td><td>3323.</td><td>5151.</td><td>20922.</td><td>30242.</td></tr><tr><td>σ<sub>HU</sub> initial</td><td>0.196</td><td>11.76</td><td>14.65</td><td>29.52</td><td>35.50</td></tr><tr><td>Nb d'itérations</td><td>0</td><td>4</td><td>2</td><td>1</td><td>4</td></tr><tr><td>W<sub>s</sub> final</td><td>3.0</td><td>2.28</td><td>2.82</td><td>2.93</td><td>2.32</td></tr><tr><td>W<sub>p</sub> final</td><td>320.0</td><td>320.6</td><td>312.5</td><td>321.3</td><td>319.0</td></tr><tr><td>S finale</td><td>0.931</td><td>0.0041</td><td>1.789</td><td>0.816</td><td>0.0134</td></tr><tr><td>σ<sub>HU</sub> final</td><td>0.196</td><td>0.013</td><td>0.273</td><td>0.184</td><td>0.024</td></tr><tr><td>Divergence</td><td>oui</td><td>non</td><td>oui</td><td>oui</td><td>non</td></tr></tbody></table>

**Tableau 4:** Récapitulatif de l'assimilation variationnelle de l'eau du sol. On donne les valeurs initiales de l'algorithme de minimisation, la valeur initiale de la fonctionnelle S à minimiser, ainsi que l'écart-type d'erreur de l'humidité relative à 2m $\sigma_{HU} = \sqrt{S/N}$. Ces quantités sont fournies après n itérations du processus de descente.
<table><thead><tr><td>Simulation</td><td>W<sub>s</sub></td><td>W<sub>p</sub></td></tr></thead><tbody><tr><td>1</td><td>2.36</td><td>314.</td></tr><tr><td>2</td><td>2.40</td><td>315.</td></tr><tr><td>3</td><td>2.40</td><td>315.</td></tr><tr><td>4</td><td>2.10</td><td>199.</td></tr><tr><td>5</td><td>2.11</td><td>200.</td></tr><tr><td>6</td><td>2.20</td><td>201.</td></tr></tbody></table>

**Tableau 5:** Valeurs de l'eau du sol après 24h d'assimilation séquentielle utilisant la méthode de Geleyn

<table><tr><td>α</td><td>1.412</td><td>2.974</td><td>2.085</td><td>1.084</td></tr><tr><td>h</td><td>6</td><td>12</td><td>18</td><td>24</td></tr></table>

**Tableau 6:** Rapport entre l'erreur de prévision de l'humidité relative de surface et l'erreur de prévision de l'humidité relative au premier niveau atmosphérique aux heures synoptiques.

<table><thead><tr><td>Simulation</td><td>W<sub>s</sub></td><td>W<sub>p</sub></td></tr></thead><tbody><tr><td>1</td><td>2.36</td><td>314.</td></tr><tr><td>2</td><td>2.40</td><td>315.</td></tr><tr><td>3</td><td>2.40</td><td>316.</td></tr><tr><td>4</td><td>2.40</td><td>200.</td></tr><tr><td>5</td><td>2.40</td><td>200.</td></tr><tr><td>6</td><td>2.40</td><td>202.</td></tr></tbody></table>

**Tableau 7:** Valeurs de l'eau du sol après 24h d'assimilation séquentielle utilisant la méthode de Geleyn modifiée
<table><caption>**Tableau 8:** Poids de l'interpolation optimale</caption><thead><tr><th>heure</th><th>λ<sup>1</sup><sub>s</sub></th><th>λ<sup>2</sup><sub>s</sub></th><th>λ<sup>1</sup><sub>p</sub></th><th>λ<sup>2</sup><sub>p</sub></th></tr></thead><tbody><tr><td>6</td><td>1.1492</td><td>0.0917</td><td>266.55</td><td>13.62</td></tr><tr><td>12</td><td>-2.661</td><td>-0.2103</td><td>-399.70</td><td>-37.38</td></tr><tr><td>18</td><td>-0.0841</td><td>0.0284</td><td>-86.77</td><td>-9.018</td></tr><tr><td>24</td><td>1.097</td><td>0.1049</td><td>123.53</td><td>19.94</td></tr></tbody></table>

<table><caption>**Tableau 9:** Coefficients de corrélation d'ordre zéro ($ρ_{ij}$) et de corrélation multiple ($r_{ijk}$) entre les variables atmosphériques T et HU et l'eau du sol</caption><thead><tr><th>heure</th><th>ρW<sub>s</sub>,T</th><th>ρW<sub>s</sub>,HU</th><th>ρW<sub>p</sub>,T</th><th>ρW<sub>p</sub>,HU</th><th>ρT,HU</th><th>rW<sub>s</sub>,T,HU</th><th>rW<sub>p</sub>,T,HU</th></tr></thead><tbody><tr><td>6</td><td>0.13</td><td>0.71</td><td>0.63</td><td>0.54</td><td>-0.31</td><td>0.80</td><td>0.99</td></tr><tr><td>12</td><td>-0.80</td><td>0.76</td><td>-0.61</td><td>0.51</td><td>-0.98</td><td>0.82</td><td>0.89</td></tr><tr><td>18</td><td>-0.78</td><td>0.81</td><td>-0.69</td><td>0.67</td><td>-0.99</td><td>0.82</td><td>0.70</td></tr><tr><td>24</td><td>-0.76</td><td>0.86</td><td>-0.64</td><td>0.77</td><td>-0.97</td><td>0.92</td><td>0.92</td></tr></tbody></table>

<table><caption>**Tableau 10:** Valeur de l'eau du sol après 24h d'assimilation séquentielle par interpolation optimale</caption><thead><tr><th>Simulation</th><th>W<sub>s</sub></th><th>W<sub>p</sub></th></tr></thead><tbody><tr><td>1</td><td>2.36</td><td>314.</td></tr><tr><td>2</td><td>2.76</td><td>316.</td></tr><tr><td>3</td><td>2.03</td><td>298.</td></tr><tr><td>4</td><td>2.95</td><td>315.</td></tr><tr><td>5</td><td>2.63</td><td>316.</td></tr><tr><td>6</td><td>2.57</td><td>318.</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Evolution diurne du flux de chaleur latente pour les six simulations de référence (voir Tableau 3) et mesuré par la SAMER 7 le 19 mai 1986 à Bats (trait pointillé).</figcaption></figure>
EAU DU RESERVOIR SUPERFICIEL EN MM

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

EAU DU RESERVOIR SUPERFICIEL EN MM

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

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

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

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

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

**Figure 2:** Evolution diurne de l'eau du réservoir superficiel pour les six simulations de référence.
<figure><img src="image_9.png" /><figcaption>Figure 3: Influence du facteur $\epsilon$ sur le calcul des éléments de la matrice jacobienne J.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 4:** Eau du sol à chaque itération de l'algorithm de minimisation en partant des conditions de la simulation 3.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 5:** Fonction de coût à chaque étape de l'algorithm de minimisation en partant des conditions de la simulation 3.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 6a:** Evolution de l'erreur d'humidité atmosphérique (A) et de l'erreur de l'humidité superficielle (B) pour la simulation 3.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 6b:** Evolution de l'erreur d'humidité atmosphérique (A) et de l'erreur de l'humidité superficielle (B) pour la simulation 2.</figcaption></figure>
<figure><img src="image_15.png" /></figure>

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

**Figure 7a** Evolution de l'eau du sol dans une assimilation séquentielle par la méthode de Geleyn partant des conditions hydriques de la simulation 5. On a reporté en pointillés l'évolution de l'eau du sol pour la simulation de référence (1).
<figure><img src="image_17.png" /></figure>

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

**Figure 7b** Evolution de l'eau du sol dans une assimilation séquentielle par la méthode de Geleyn modifiée partant des conditions hydriques de la simulation 5. On a reporté en pointillés l'évolution de l'eau du sol pour la simulation de référence (1).
<figure><img src="image_19.png" /></figure>

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

**Figure 8a** Evolution de l'eau du sol dans une assimilation séquentielle par interpolation optimale partant des conditions hydriques de la simulation 2. On a reporté en pointillés l'évolution de l'eau du sol pour la simulation de référence (1).
<figure><img src="image_21.png" /></figure>

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

**Figure 8b** Evolution de l'eau du sol dans une assimilation séquentielle par interpolation optimale partant des conditions hydriques de la simulation 3. On a reporté en pointillés l'évolution de l'eau du sol pour la simulation de référence (1).
<figure><img src="image_23.png" /></figure>

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

**Figure 8c** Evolution de l'eau du sol dans une assimilation séquentielle par interpolation optimale partant des conditions hydriques de la simulation 4. On a reporté en pointillés l'évolution de l'eau du sol pour la simulation de référence (1).
<figure><img src="image_25.png" /></figure>

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

**Figure 8d** Evolution de l'eau du sol dans une assimilation séquentielle par interpolation optimale partant des conditions hydriques de la simulation 5. On a reporté en pointillés l'évolution de l'eau du sol pour la simulation de référence (1).
<figure><img src="image_27.png" /></figure>

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

**Figure 8e** Evolution de l'eau du sol dans une assimilation séquentielle par interpolation optimale partant des conditions hydriques de la simulation 6. On a reporté en pointillés l'évolution de l'eau du sol pour la simulation de référence (1).