THEORIE MATHEMATIQUE DU REBOND POST-GLACIAIRE

Isostasie. Cas général et formalisme.











Calcul de déformations de la Terre dues aux effets de charge. Cas général et formalisme.

La prédiction de déformations de la Terre solide associées à une charge posée à la surface est un problème classique de géophysique (voir par exemple Farrell, 1972. Comme, depuis les années 1980, et surtout 1990, le développement des techniques de positionnement spatiales telles que SLR, VLBI et GPS a permis des mesures de positions et vitesses de points situés à la surface de la Terre avec une précision de plus en plus grande (autour de 1 mm pour le VLBI), les prédictions de réponse de la surface terrestre à des effets de charge doivent également se faire de plus en plus précises. Il importe en effet de pouvoir incorporer des modèles de mouvements verticaux dans les inversions de positions et vitesses des sites GPS ou VLBI, et ce avec une précision suffisante. Cela concerne aussi bien les réponses à la pression atmosphérique (vanDam, 1987) à une surcharge océanique (Scherneck, 1991), ou aux effets visco-élastiques de la déglaciation du Pleistocène (Tushingham, 1991). Dans la suite, nous présentons de façon progressive le formalisme utilisé dans la plupart des modélisations récentes de rebond post-glaciaire (fonctions de Green et décomposition en harmoniques sphériques). On considérera différents cas de comportement du manteau supérieur, selon une rhéologie élastique, et visco-élastique respectivement.


L'approche des fonctions de Green.

D'après Farrell, 1972, cette approche consiste dans un premier temps à calculer les fonctions de Green pour les déplacements radial et tangentiel associés au potentiel de marée et à une charge sur une Terre à symétrie sphérique ayant une rhéologie pûrement élastique. On sait depuis longtemps que le comportement du manteau est assez éloigné de celui d'un solide élastique (viscosité comprise entre 1020 et 1023 Pas approximativement), mais pour un effet de charge ayant une constante de temps suffisamment brève (de la semaine à quelques dizaines d'années), l'approximation élastique peut suffire à modéliser la réponse en surface avec une exactitude satisfaisante.
Dans tous les calculs de réponse terrestre à une sollicitation externe, on doit en fait tenir compte de deux effets agissant sur la surface et le potentiel de gravitation : celui dû au potentiel de marée et celui dû à la charge en elle-même. On commence par l'exposé de cet effet de marée terrestre, qui sera ensuite incorporé aux calculs globaux.
La réponse du modèle de Terre à une charge massique arbitraire posée à la surface peut ensuite être calculée par le biais d'un produit de convolution en surface avec la fonction de Green appropriée. En pratique, la fonction de Green est discrétisée sur l'ensemble de la surface de la Terre, sous forme d'éléments finis, la taille de la surface élémentaire étant conditionnée par la résolution de la charge appliquée.
Le formalisme spectral se fonde sur les nombres de Love Love, 1909. Ces nombres gouvernent la réponse d'un modèle planétaire à un forçage de potentiel gravitationnel ou de charge en surface, respectivement.


Réponse au potentiel de marée.

Les nombres de Love HL, LL, et KLreprésentent le coefficient d'harmoniques sphériques de degré L des déplacements radial et tangentiel, et de la perturbation du potentiel de gravitation dûs à la redistribution de masse (marées terrestres) provoquée par le potentiel de marée. Si l'on désigne par W(, ,t). Ce potentiel de marée terrestre au point de la surface terrestre de latitude , de longitude et à l'instant t, la réponse de la Terre à ce potentiel s'écrira :

U =
YLm(, ) correspond aux harmoniques sphériques normalisées sur la surface de la Terre, et où le symbole désigne l'opérateur gradient à deux dimensions, soit :
=


Réponse à une charge. Cas pûrement élastique.

Dans ce cas, les nombres de Love s'écrivent HL, LL et KL, où fait référence à la rhéologie élastique adoptée, et L au degré de la décomposition en harmoniques sphériques, et désignent, respectivement, les coefficients d'harmoniques sphériques des déplacements radial et tangentiel, et de la perturbation gravitationnelle provoqués par la charge.

Les fonctions de Green correspondant à cet effet de charge s'écrivent :

G =

P =
respectivement, où a désigne le rayon moyen de la Terre, M sa masse, est un vecteur unitaire tangent à la surface, selon le grand cercle qui passe par le point de charge (', ') vers le point d'observation (, ), et désigne la distance angulaire entre les deux points.
Notons que peut être donnée par :
= cos cos' + sin sin' cos(-')


Enfin, PL représente la polynôme de Legendre de degré L.

Dans le cas particulier élastique, si l'on désigne par G() une fonction de Green, la réponse de la Terre à une masse surfacique appliquée au point (', ') et au temps t sera donnée par :
R =

G =
d désigne la surface élémentaire autour du point ( ', ').

La réponse R(,,t) sera donc une quantité scalaire pour les fonctions de Green radiale et de potentiel, et un vecteur pour la réponse tangentielle correspondant à la fonction.
En appliquant directement les équations précédentes, et à condition de discrétiser suffisamment la charge, puis de procéder à l'intégration par éléments finis, on obtient la réponse de la Terre (e.g. (Jentzsch, 1985).

Réponse à une charge. Cas visco-élastique.

Dans le cas de déformations de la surface associées à l'effet de charge du rebond post-glaciaire, l'approximation élastique est insuffisante. Du fait des constantes de temps, longues, associées au comportement d'une calotte de glace, les nombres de Love correspondants deviennent dépendant du temps. D'après le principe de correspondance (Biot, 1954), on peut les écrire sous la forme :



Le formalisme explicité par exemple par (Peltier, 1985) permet de calculer les amplitudes de mode normal et les temps de relaxation correspondant aux pulsations propres.
On écrit ensuite de la même façon que dans le cas pûrement élastique les fonctions de Green associées :


L'intégration de la réponse sur l'ensemble de la surface terrestre se fait de la même façon que dans l'équation ?, à ceci près que la fonction de Green dépend du temps et qu'il faut distinguer le temps t' auquel s'applique la charge surfacique du temps t d'observation. L'intégration devient un double produit de convolution sur la surface de la Terre et le temps, soit :


La charge est discrétisée, en général sous forme de disques d'épaisseur constante, avant intégration. Cette méthode, incluant la prise en compte de la redistribution océanique, ainsi que les conséquences de la déformation crustale, a été utilisée à maintes reprises (par exemple (Wu et Peltier, 1983; James et Morgan 1990, ou Spada et al., 1992)


Réponse de la Terre à une charge arbitraire : une approche spectrale.

Cette méthode, reprenant le formalisme de base exposé précédemment, a été exposée par Mitrovica et al., 1994a. Elle s'applique dans le cas d'un modèle de Terre à symétrie sphérique, et consiste à décomposer la charge en surface selon les harmoniques sphériques, plutôt que de la discrétiser, soit :



Si l'on reprend les fonctions de Green de l'équation ? et la réponse terrestre de l'équation ?, on va obtenir par exemple pour la réponse en déplacement tangentiel :


Comme et que l'opérateur ne s'applique qu'aux quantités , et comme , on peut écrire :


Ceci compte tenu de l'orthonormalisation des harmoniques sphériques. De la même façon, on aura :


pour le déplacement radial à la surface de la Terre, et :


pour l'effet sur le potentiel de gravité.

L'énorme avantage de cette décomposition de la charge en harmoniques sphériques est que cela permet d'adapter la résolution spatiale du calcul par simple troncature du degré maximal de décomposition. Cette stratégie se révèle particulièrement payante lorsque la résolution de la charge varie spatialement, ou encore lorsqu'on reprend un calcul déjà effectué, avec un modèle de charge amélioré, dont la résolution s'est affinée.


Application pratique. Cas d'une rhéologie élastique.

Dans la pratique, la charge peut être constituée d'une charge atmosphérique, océanique, ou d'une surcharge glaciaire, ou encore une combinaison de plusieurs de ces éléments.
En reprenant la décomposition de la charge, la réponse locale de la Terre, prenant en compte aussi bien l'effet de marée terrestre que l'effet de charge, s'écrit :



et


où le vecteur est défini par :


On a considéré jusqu' à présent que les nombres de Love, qu'ils s'appliquent au potentiel de marée ou à l'effet de charge, sont indépendants de la fréquence. On sait que ce n'est pas le cas, et qu'il faudrait écrire :


Il faut ensuite introduire les expressions explicites du potentiel de marée et de l'effet de charge.
Le potentiel de marée est exprimé sous forme d'une somme de termes sinusoïdaux, soit pour le terme de degré L et d'ordre m :


où , et représentent la fréquence, la phase et l'amplitude du nème composants sur un total de nm. De la même façon :


et on a au total les expressions :


et


Il faut remarquer que ce formalisme n'est pas adapté à tous les cas de charge. Par exemple, pour la réponse à la fonte actuelle des calottes et des glaciers, qui a des constantes de temps très longues, on peut considérer que les nombres de Love sont indépendants de la fréquence.

Résolution des équations d'équilibre du niveau des mers.

Les équations ?, ? et fournissent le système permettant d'obtenir la réponse à l'effet de charge compte tenu des modifications du potentiel de perturbation gravitationnel. La perturbation du potentiel de gravitation total peut se décomposer en harmoniques sphériques :



où chaque terme inclut la réponse au potentiel de marée, la charge de l'océan et la charge autre que celle de l'océan (par exemple l'allègement de la charge glaciaire).
On définit la fonction océan, , comme égale à l'unité aux endroits recouverts par la mer, nulle partout ailleurs (soit les composantes de sa décomposition en harmoniques sphériques). On désigne par ? l'épaisseur de l'eau dans l'océam , ? les composantes associées, et ?, où w est la densité de l'eau de mer.
L'équation qui gouverne la hauteur d'eau dans les océans s'écrit au total :


Physiquement, elle reflète le fait que la charge océanique perturbe le champ de gravité, qui a son tour vient perturber l'équilibre océanique. L'expression pour la quantité ? peut être obtenue en imposant une contrainte de conservation de la masse sur l'ensemble de l'océan.
En intégrant ? la surface ?, en utilisant le fait que , on obtient :


et


est obtenu par conservation de masse, en fonction des autres effets de charge.

Application à la fonte glaciaire actuelle.

On considère la redistribution océanique résultant de la fonte des glaciers et calottes. On a alors :



et


i est la densité de la glace et où . est la composante de degré L, et d'ordre m de l'épaisseur de la glace. Pour la redistribution océanique, une première approximation est donnée par la redistribution eustatique :


qui se place dans le cas d'une redistribution de l'eau de fonte d'une manière indépendante de la géographie.
Dans le cas de cette variation, séculaire, on calcule la réponse globalement plutôt que pour une fréquence spécifique.

Application dans le cas d'une rhéologie visqueuse.

On se place à nouveau dans le cas d'une Terre sphérique avec une rhéologie mantellique visco-élastique. Les équations ?, ?, et ?, mènent à :



et


avec


Le premier terme du membre de droite de représente la composante élastique de la réponse, le terme suivant la composante non-élastique.
Si on applique cette méthode à la déglaciation du Pléistocène (fonte des grandes calottes de l'hémisphère Nord et d'une partie de l'Antarctique), on représente généralement la décroissance de ces calottes au cours du temps comme une succession de paliers. Cette modélisation n'entraîne pas de perte de généralité dans la mesure où les intervalles de temps peuvent être aussi petits qu'on le désire. En appliquant la même discrétisation à la charge en surface, on peut écrire :


H représente la fonction de Heaviside et N le nombre total d'incréments. On a alors :


En réinjectant ces quantités dans les équations ? et ? on obtient :


et


Les sommes sur n incluent tous les éléments de charge dont le début se situe avant le temps d'observation t. N(t) représente donc le nombre d'incréments sur un total de N, qui satisfont ?.
De façon générale, on peut séparer les effets de charges provenant de la redistribution océanique, des autres effets, dits résiduels. On peut donc écrire :


et


où les Ylm sont les coefficients de la décomposition en harmoniques sphériques de la hauteur d'eau dans l'océan à léquilibre, et les ? les coefficients associés pour le changement incrémental.

Application à la déglaciation Pléistocène.

On utilise les expressions précédentes de ? et ? avec la décomposition de la charge exprimée ci-dessus. Dans ce cas particulier, le terme résiduel concerne l'épaisseur de glace, soit :



et


Ylm représentent les coefficients d'harmoniques sphériques de la charge glaciaire, et ? les changements incrémentaux associés pour le pas de déglaciation.

Calcul à partir du mouvement du Pôle ou de l'accélération de gravité.

Le développement de modèles de rebond à partir d'une Terre présentant une structure visco-élastique stratifiée a permis de montrer que, contrairement à l'idée longtemps admise, la déglaciation continue à produire des effets sur la dérive du pôle longtemps après que la variation de masses ait cessé. La quantité , où est la vitesse angulaire actuelle de rotation de la Terre, et l'accélération angulaire selon l'axe de rotation, est directement proportionnelle à la dérivée temporelle du terme axial de degré 2 du développement en harmoniques sphériques du champ de gravité. Ce terme est évalué par les mesures de géodésie spatiale (tirs SLR sur les satellites de mesure du champ de gravité Starlette et Lageos) à des valeurs comprises entre -2,5 +/- 0,3 10-11 /an et -3,5 +/-0,3 10-11 /an (Peltier et Jiang, 1996; Rubincam, 1984).
Les variations de la rotation de la Terre provoquées par la déglaciation affectent les composantes du tenseur d'inertie (?) par le biais des effets de charge. La conservation du moment d'inertie de la Terre nécessite un ajustement de la vitesse de rotation (de composantes ?. Les deux quantités sont reliées entre elles par les équations d'Euler, qui traduisent le principe de conservation. Dans un repère fixe lié à la Terre, on peut écrire (Peltier et Jiang, 1996) :



où ? est le tenseur de Levi-Cevita. Si on suppose que ? reste toujours très proche de sa valeur initiale ?, l'équation? peut être linéarisée en introduisant le terme de perturbation :


A, B, et C sont les trois moments d'inertie principaux, les perturbations adimensionnées des composantes de la vitesse angulaire ?, et les perturbations d'inertie. En substituant ? dans ? on obtient :


avec les fonctions d'excitation données par :


La première étape dans le calcul de l'effet de charge sur ? consiste à évaluer les expression de ?, ? et ?. Pour un modèle visco-élastique, la perturbation d'inertie est la somme d'un terme dit ;SPMlt;;SPMlt; de Terre rigide;SPMgt;;SPMgt; résultant du seul effet de la charge de surface (indépendant de la rhéologie) et un terme provenant de la déformation terrestre due à l'effet de charge. La part rigide inclut la contribution de la glace et de l'océan. Si on considère que la charge en surface est sur une couche très fine par rapport au rayon de la Terre, on peut approcher le terme rigide de la perturbation d'inertie par :


en intégrant sur toute la surface de la Terre, dont a est le rayon moyen en harmoniques sphériques on obtient pour le terme ? :


ce qui signifie que les coefficients de degré 2 de la charge de glace et de l'océan suffisent au calcul.
Pour exprimer ? en fonction de ?, il suffit ensuite de résoudre l'équation ? avec la perturbation ? donnée par :


dans laquelle ? est la convolution de ? par le nombre de Love de degré 2. Ce nombre de Love de degré 2 est donné par :


En substituant cette dernière expression dans l'équation ? et en intégrant, on obtient l'expression de ? :


? se déduit de ? par le biais d'une constante multiplicative :


M est la masse de la Terre.
La figure de Peltier et Jiang, 1996, récapitule les valeurs du ? actuelles observées, à partir des analyses des données anciennes d'éclipses pour Stephenson, 1995, à partir d'observations satellitaires pour les autres.
L'expression explicite du mouvement du pôle est plus difficile à obtenir. Il s'agit toujours de résoudre le système des équations en tenant compte de ?. Le moyen de plus simple semble être l'utilisation de la transformée de Laplace, et la résolution dans le domaine de Laplace de ces équations sous la forme de ?, ce qui donne, dans le domaine temporel :


où le symbole ? désigne le produit de convolution. Les fonctions ?, ?, H' et H'' ont pour expression, dans le domaine de Laplace :


avec


et


? est le nombre de Love fluide de degré 2 défini par ?, et ? est le nombre de Love de l'effet de marée de degré 2. Il s'agit ensuite d'obtenir des expressions explicites pour les fonctions ?, ?, H' et H'' dans le domaine de Laplace, puis d'obtenir leur valeur dans le domaine temporel par une transformée inverse.






Tous les calculs exposés ici sont extraits de la thèse de M.N. Bouin, effectuée au laboratoire de géologie de l'ENS, sous la direction conjointe de Ch. Vigny (ENS) et C. Boucher (IGN), et soutenue le 14 Octobre 1999, devant un Jury composé de MM F. Barlier (président), K. Lambeck (rapporteur), M. Kasser (rapporteur), F. Remy (examinateur), B. Ambrosius (examinateur), C. Boucher, et Ch. Vigny.