le massif (contraintes naturelles et perturbées par la mine), on peut difficilement confirmer
avec certitude l'endommagement de la paroi du forage. On peut néanmoins affirmer que cet
endommagement s'il a lieu est certainement faible et doit être localisé à la paroi de la
chambre puisque les contraintes radiales et orthoradiales en coordonnées cylindriques
évoluent en 1/r². Néanmoins, afin d'estimer l'extension d'une zone endommagée par le
forage, deux cas extrêmes sont considérés. Le forage est réalisé instantanement :

(a) sans fluide (pression nulle en paroi). Le calcul fera donc intervenir le critère d'endommagement du sel "sec";

(b) avec un fluide sous pression (dont la valeur serait celle du premier palier). Le calcul utilise donc le critère d'endommagement du sel "humide";

En l'absence de fluide dans le forage, un calcul numérique réalisé à partir d'un modèle de comportement élastoplastique montre que la zone endommagée au voisinage d'un puits ne dépasse guère 50 % du rayon du puits (Thorel, 1994). Compte tenu du rayon de la chambre (7,3 cm), la zone endommagée calculée s'étend sur 4 cm au maximum. La valeur de ce rayon doit être comparée à la distance de pénétration* estimée à 15 cm (pour les essais en saumure). Le facteur 3,75 entre ces deux grandeurs suggère que l'endommagement même s'il existe ne peut expliquer, à lui seul la perméabilité mesurée.

En présence d'un fluide (pression de fluide non nulle en paroi), l'extension de cette
zone endommagée peut être facilement estimée à l'aide des deux hypothèses :

(i) le seuil d'endommagement étant faiblement dépendant de la contrainte moyenne, on suppose qu'il obéit à un critère de Tresca :

$$
F = \operatorname{Sup}_{i,j} |\sigma_i - \sigma_j| - 2C
$$

$\sigma_i$ et $\sigma_j$ : contraintes principales extrema; C : cohésion prise égal à 5 MPa, d'après la figure 13.

(ii) d'une part, le contraste de perméabilité étant très important entre la zone endommagée et la zone saine (cf figure 10 du chapitre I) et d'autre part, l'extension de la zone endommagée étant faible (vérifiable *a posteriori*), les pertes de charge dans cette zone sont négligées.

Le calcul de ce rayon d'endommagement R (formellement un calcul de rayon plastique) conduit à l'expression suivante (cf annexe III) :

$$
\ln \frac{R}{r_0} = - \frac{1}{2C} \left[ C + \sigma_0 + P_1 - b \left( \frac{1 - 2v_0}{1 - v_0} \right) \left( \frac{P_1 - P_0}{2} \right) \right]
$$

r₀ : rayon du forage; P₀ : pression interstitielle initiale; P₁ : pression en paroi; v₀ : coefficient de Poisson drainé.

D'après cette expression, on note que plus la pression appliquée en paroi est
importante, moins la zone endommagée est étendue : la pression en paroi joue un rôle de
soutènement. Avec les valeurs notées dans le tableau 8, une application numérique donne
un rapport R/r₀ de 1,3 pour une pression de fluide en paroi de 2 MPa (premier palier
appliqué, pression en paroi de fluide minimum). Compte tenu du rayon de la chambre (7,3
cm), la zone endommagée calculée s'étend sur 2 cm. Là encore, cette valeur doit être
comparée à la distance de pénétration (15 cm). Le facteur 7,5 entre ces deux grandeurs
souligne une nouvelle fois, le rôle limité de l'endommagement.

<footnote>* Il convient de souligner ici que la "distance de pénétration" utilisée dans ce calcul est un *minorant* du rayon de la zone influencée par l'écoulement, lorsque le milieu poreux est initialement saturé en saumure. Ceci contribue alors à justifier plus facilement l'hypothèse c). Nous y reviendrons plus loin dans le paragraphe 2.1 de ce chapitre relatif à la modélisation des essais.</footnote>
Vu les distances impliquées dans ces calculs (plusieurs centimètres), la zone endommagée ne doit concerner que quelques grains qui auront probablement vu leur joints se dilater, voir se rompre. Le réseau hydraulique correspondant à cette zone serait alors associé à une perméabilité très importante devant celle du milieu sain loin du forage, de sorte que cette perméabilité ne serait impliquée que dans les premiers instants de la partie transitoire de la courbe de débit. Par ailleurs, les masses injectées de saumure mesurées, durant les périodes où la thermorégulation a bien fonctionné et dans la phase du régime permanent ne montre aucune évolution significative de la perméabilité avec la distance à la paroi.

## Endommagement dû au chargement thermique

L'endommagement d'origine thermique de la paroi du puits pourrait également provoquer des changements de la perméabilité du massif. D'une manière générale, il résulte d'une rupture en traction de la paroi consécutivement à un refroidissement du puits. En effet, la contrainte orthoradiale $\sigma_{\theta}$ et la contrainte vertical $\sigma_{z}$ initialement en compression, pourrait se mettre en traction, suite au refroidissement, et dépasser la faible résistance en traction du sel voisine de 2 MPa (égale à 1,65 MPa pour le sel des MDPA cf partie 1.2 Chapitre II).

Dans un premier temps, considérons que le massif est poro-élastique. Les composantes du tenseur des contraintes en paroi du forage s'écrivent alors :

$$ 
\begin{aligned}
\sigma_r &= -P_I \\
\sigma_{\theta} &= 2\sigma_o + P_I - \frac{1-2v_o}{1-v_o}[b(P_I - P_o) + 3\alpha_o K_o (T_I - T_o)] \\
\sigma_z &= \sigma_o - \frac{1-2v_o}{1-v_o}[b(P_I - P_o) + 3\alpha_o K_o (T_I - T_o)]
\end{aligned}
$$

$\alpha_o$ : coefficient de dilatation linéique drainé, $K_o$ : module de compressibilité drainé, $T_I$ : la température à la paroi, $T_o$ : la température initiale du massif. $P_I$ : pression en paroi ; $P_o$ : pression de pores initiale.

Pour ne pas engendrer une rupture en traction, la contrainte orthoradiale et la contrainte verticale doivent satisfaire l'inégalité :

$$ 
\begin{gathered}
\sigma_{\theta} \le 1,65 \text{ MPa et } \sigma_{z} \le 1,65 \text{ MPa} \\
\begin{aligned}
2\sigma_o + P_I - \frac{1-2v_o}{1-v_o}[b(P_I - P_o) + 3\alpha_o K_o (T_I - T_o)] &\le 1,65 \text{ MPa} \\
\sigma_o - \frac{1-2v_o}{1-v_o}[b(P_I - P_o) + 3\alpha_o K_o (T_I - T_o)] &\le 1,65 \text{ MPa}
\end{aligned}
\end{gathered}
$$

La variation de température doit alors satisfaire :

$$ (T_o - T_I) \le - \frac{2\sigma_o - P_I - 1,65}{3\alpha_o K_o} \frac{1 - v_o}{1 - 2v_o} + \frac{b}{3\alpha_o K_o} (P_I - P_o) $$

et

$$ (T_o - T_I) \le - \frac{\sigma_o - 1,65}{3\alpha_o K_o} \frac{1 - v_o}{1 - 2v_o} + \frac{b}{3\alpha_o K_o} (P_I - P_o) $$
Procédons à une application numérique avec les valeurs déjà retenues pour les calculs précédents $\sigma_o=10 \text{ MPa}, v_o=0,3; \alpha_o=10^{-5} \text{ $C^{-1}$}; K_o=22,5 \text{ GPa}$ (calculé à partir des valeurs des paramètres } K_s=25 \text{ GPa et } b=0,1 \text{ donnés par Mc Tigue, 1986); on obtient :

$(T_o-T_1)< 69 \text{ $^\circ C$}, \text{ et } (T_o-T_1)< 30 \text{ $^\circ C$} \text{ avec } P_1=0 \text{ MPa (le cas le plus défavorable).}$

Un tel refroidissement si intense par rapport à l'état initial n'a évidemment pas eu lieu. Néanmoins, compte tenu de la nature viscoplastique du sel, la relaxation de l'état de contrainte conduit à une situation beaucoup plus défavorable vis-à-vis du risque de rupture en traction. Le déviateur limite à la paroi, s'écrit pour une géométrie cylindrique (cf annexe) :

$$ \sigma_r - \sigma_\theta = -\frac{2}{n}(\sigma_o + P_1) $$

La contrainte orthoradiale $\sigma_\theta$ et la contrainte verticale $\sigma_z$ vérifie alors :

$$ \sigma_{\theta} = -P_{1} + \frac{2}{n}(\sigma_{o} + P_{1}) \quad \text{et} \quad \sigma_{z} = \sigma_{o} + 2v_{o}(\frac{1}{n} - 1)(\sigma_{o} + P_{1}) - b(1 - 2v_{o})(P_{1} - P_{o}) $$

n : exposant de contrainte dans la loi de fluage de Norton-Hoff.

Procédons, une nouvelle fois encore, à une application numérique (avec n=4 ; $\sigma_o=-10 \text{ MPa}$). Avec $P_1=0$ (cas le plus défavorable), on obtient $\sigma_\theta=-7 \text{ MPa}$ et $\sigma_z=-6 \text{ MPa}$. Lorsqu'un refroidissement instantané se produit en paroi du forage ($\delta T<0$), une contrainte orthoradiale et une contrainte verticale de traction d'origine thermo-élastique prennent naissance :

$$ \delta\sigma_{\theta}^{th} = \delta\sigma_{z}^{th} = \frac{E_{o}}{1 - v_{o}}\alpha_{o}\delta T > 0 $$

Si la somme de ces contraintes thermiques et des valeurs des contraintes $\sigma_\theta$ et $\sigma_z$ calculées est supérieure à la résistance en traction, alors la rupture par traction a lieu. Ainsi, afin d'éviter la rupture, les inégalités suivantes doivent être satisfaites :

$$ -\frac{E_o}{1 - v_o}\alpha_o \delta T + \sigma_\theta \le 1,65 \text{ MPa} \quad \text{ou} \quad -\frac{E_o}{1 - v_o}\alpha_o \delta T + \sigma_z \le 1,65 \text{ MPa} $$

c'est à dire :

$$ \delta T \geq -\frac{1 - v_o}{\alpha_o E_o} \text{Inf}(1,65 - \sigma_\theta, 1,65 - \sigma_z) $$

Une application numérique donne avec $E_o=20 \text{ GPa}$ (calculée à partir de $v_o$ et de $K_o$), $v_o=0,3$; $\alpha_o=10^{-5} \text{ $C^{-1}$}$:

$$ \delta T \geq 27 \text{ $^\circ C$} $$

Une nouvelle fois encore, à aucun moment de l'expérimentation, la chambre d'essai n'a subi un refroidissement de cette amplitude. Les variations de température dans le forage au niveau de la zone d'essai située à plus de seize mètres de la niche ont toujours été inférieures à 10 °C. Le risque d'endommagement par un tel processus peut donc être écarté.

**d- l'écoulement est radial** (perpendiculairement à l'axe vertical de la chambre). Si on compare la distance de pénétration de la saumure (de l'ordre de quinze centimètres) et la hauteur de la chambre (75 cm), cette hypothèse parait plausible au voisinage de la chambre d'injection.

**e- le Volume Élémentaire Représentatif associé à la perméabilité est suffisamment petit** devant la dimension caractéristique du problème hydraulique étudié (le rayon de la chambre ou mieux le "rayon d'action" de l'écoulement, traduisant l'extension de la
zone influencée par la percolation de saumure) Cette hypothèse doit être vérifiée pour que la notion de perméabilité puisse s'appliquer pour l'interprétation des essais. Une façon simple de justifier cette hypothèse est de comparer *a posteriori* le "rayon d'action" de l'écoulement obtenu grâce à un modèle (cf. paragraphe 1.2) et la taille du VER estimée à plusieurs centimètres vu les caractéristiques géométriques du réseau conducteur (taille moyenne des grains, forme des microcanaux, longueur moyenne, épaisseur moyenne, densité des canaux, .. cf. Lebrun, 1993). Cette dernière estimation de la taille du VER sera confirmée dans le paragraphe 1.2 du chapitre IV.

## 1.2 Validité de la loi de Darcy pour la saumure

A l'aide des hypothèses précédentes, la perméabilité du sel peut être calculée en résolvant l'équation de diffusion de la pression interstitielle P dans un massif poro-viscoplastique :

$$ \frac{k}{\eta} \nabla^2 P = \frac{1}{M} \frac{\partial P}{\partial t} + b \frac{\partial e}{\partial t} + \underbrace{\frac{\partial \phi^{VP}}{\partial t} - b \frac{\partial e^{VP}}{\partial t}}_{\text{contribution de la viscoplaticité}} $$

k : perméabilité intrinsèque (m²); η : viscosité dynamique (Pa·s); ∇ : l'opérateur laplacien; P : pression du fluide (Pa); M : module de Biot (Pa); b : coefficient de Biot; e et eP : déformation volumique totale et viscoplastique respectivement, φVP : porosité viscoplastique.

Dans la partie 1. de ce chapitre, on a pu voir que, si l'écoulement a lieu dans une zone non endommagée mécaniquement (sel sain), l'équation de diffusion précédente pouvait se simplifier et prendre la forme suivante :

$$ \frac{\partial(P - P_o)}{\partial t} = D \nabla^2 (P - P_o) $$

p : pression interstitielle, p₀ : pression interstitielle initiale, D diffusivité définie par :

$$ D = \frac{k}{\eta} \frac{3K_o + 4\mu}{3K + 4\mu} \left( \frac{b - \phi}{K^s} + \frac{\phi}{K^{fl}} \right)^{-1} $$

k : perméabilité intrinsèque, η : viscosité dynamique du fluide (1,7 10³ Pa·s à 25 °C), Kfl : module de compressibilité du liquide (3,7 GPa), Ks : module de compressibilité de la matrice solide (25 GPa), φ : porosité (0,1 %), b : coefficient de Biot pris égal à 0,1 d'après Mc Tigue (1986), K : module de compressibilité non drainé pris égal à 24,9 GPa à partir de b et K₀ : le module de compressibilité drainé (22,5 GPa calculé à partir des paramètres Ks et b).

Le débit volumique à la paroi et par conséquent le volume injecté sont donnés par une résolution numérique de l'équation de diffusion, à l'aide de la méthode des différences finies en géométrie cylindrique. L'historique du chargement de la chambre a été prise en compte : le calcul de la courbe d'injection d'un palier donné prend comme profil initial de pression, celui obtenu à la fin du palier précédent. Le volume injecté théorique est alors comparé aux données obtenues pour différentes valeurs de perméabilité et de pression de pores initiale. Le meilleur ajustement est obtenu pour le couple des valeurs de **perméabilité et de pression de pores initiale de 2.10⁻²¹ m² et 1 MPa respectivement** (figure 14).

Conformément à la première hypothèse utilisée (écoulement isotherme), on rappelle que la comparaison entre mesures et modèle a porté sur les paliers 4 et 6 MPa.
<figure><img src="image_44.png" /></figure>

Figure 14. Comparaison entre les données expérimentales et le modèle tenant compte du régime transitoire pour les valeurs de perméabilité de $2\times10^{-21} m^2$ et de pression de pores initiale de 1 MPa.

La valeur de la perméabilité intrinsèque de $2 \times 10^{-21} m^2$ est du même ordre de grandeur que celle d'un sel non endommagé mesurée sur le site du WIPP (cf. paragraphe 2.2 du chapitre I). Comme également sur le site du WIPP, la valeur de la pression de pores initiale de 1 MPa est plus faible que la pression hydrostatique probable de 6 MPa à la profondeur de la chambre d'essai (536 m). Cette valeur plus faible de 1 MPa, traduit-elle un phénomène de drainage à l'échelle de l'infrastructure minière ? Des mesures complémentaires de perméabilité à des distances variables de la paroi de la galerie apporteraient certainement des éléments de réponse à cette question.

La modélisation transitoire permet également d'estimer le "rayon d'action" hydraulique, c'est à dire l'extension de la zone influencée par l'écoulement. Le "rayon d'action" est défini arbitrairement comme la distance à la paroi de la chambre où la pression atteint la valeur de $P_1+\Delta P_{max}/10=P_1+(P_1-P_0)/10$ (P₁ : pression de la chambre et P₀ : pression de pores initiale). Le calcul donne 2 mètres, soit une distance nettement plus grande que la taille du VER et du rayon de la zone endommagée. Les hypothèses c. et e. sont ainsi vérifiées.

En outre, les résultats de cette modélisation suggèrent que des effets de non linéarités associés à la perméabilité auraient un rôle mineur ou alors se compenseraient. Ces effets sont les suivants :

- le changement de concentration en sel à saturation de la saumure consécutif aux variations de pression, induisant des phénomènes de dissolution-cristallisation,

- la cicatrisation sous contrainte,

- la variation de la contrainte effective,

Tous ces effets n'agissent pas dans le même sens. L'augmentation de la concentration à saturation suite à l'augmentation de pression du fluide (lors du changement de palier) aurait plutôt tendance à faire augmenter la perméabilité, par dissolution de la paroi des pores, alors que la cicatrisation sous contrainte et l'augmentation de la contrainte effective tendraient plutôt à diminuer la taille des canaux de l'écoulement. Examinons à la lumière de l'état actuel des connaissances l'effet que pourrait avoir chacune de ces non linéarités et proposons quelques hypothèses qui expliqueraient leur faible influence sur les résultats de l'essai.

## 1.3. Effets des non linéarités

### a- Effet de la contrainte moyenne effective

Comme nous l'avons rappelé dans la partie 3.3 du chapitre I, plus la contrainte moyenne effective appliquée est grande, plus la perméabilité du milieu poreux est faible.
Rencontré sur beaucoup de géomatériaux, cet effet a été mis en évidence in situ sur l'argile de Boom (Volkaert et al., 1995). Regardons l'effet possible de la variation de contrainte effective consécutive aux changement de paliers de pression sur la perméabilité.

Considérant le sel comme un milieu élasto-poreux et assimilant la chambre d'essai à un cylindre, la contrainte effective moyenne est exprimée en paroi (là où elle est maximum) par :

$$ \frac{1}{3} \mathrm{Tr} \sigma^{\mathrm{eff}} = \sigma_o + \frac{1}{3} \frac{1+v_o}{1-v_o} b P_1 + \frac{2}{3} \frac{1-2v_o}{1-v_o} b P_o $$

r : rayon, a : rayon du forage, $\sigma_o$ : contrainte lithostatique, $v_o$ : coefficient de Poisson drainé, $P_o$ : pression de pores initiale, $P_1$ : pression de saumure à la paroi, b : coefficient de Biot.

<figure><img src="image_45.png" /><figcaption>Figure 15. Évolution du rapport contrainte effective moyenne (sigmaeff) sur contrainte totale (sigma0) en fonction du rapport de la pression appliquée en paroi (P1) sur la contrainte totale (sigma0), pour différentes valeurs de coefficient de Biot, b = 0,1-0,5-1.</figcaption></figure>

On a représenté sur la figure 15, l'évolution du rapport contrainte effective moyenne ($\frac{1}{3} \mathrm{Tr} \sigma^{\mathrm{eff}}$) sur contrainte totale ($\sigma_o$) en fonction du rapport de la pression appliquée (P1) en paroi sur la contrainte totale ($\sigma_o$), pour différentes valeurs de coefficient de Biot. Au cours des essais, la pression en paroi a pu varier dans un rapport de 1 à 3. Il est facile de voir à l'aide de la figure 15, que l'effet d'une telle variation sur la contrainte effective est faible, compte tenu de la valeur probable du coefficient de Biot.

## b- Effet des phénomènes chimiques

Les effets chimiques considérés ici sont les changements de concentration en sel à saturation, et le phénomène de cicatrisation sous contrainte.

Concernant le deuxième phénomène, on rappelle que Peach (1991) a pu mettre en évidence une diminution d'un ordre de grandeur de la perméabilité d'échantillons préalablement mouillés à la saumure et soumis à un chargement hydrostatique durant deux semaines.

Les deux échantillons de sel de la mine de Asse étudiés par Peach avaient des perméabilités initiales de $3,3 \times 10^{-17} \text{ m}^2$ et $4,6 \times 10^{-18} \text{ m}^2$. Malgré ces valeurs de perméabilité très faibles, elles sont néanmoins respectivement 4 et 3 ordres de grandeur au-dessus de la valeur vraisemblable de la perméabilité in situ mesurée sur le site des MDPA. Ainsi, les échantillons ont été préalablement déformés axialement provoquant une dilatation volumique d'environ 1,5 %. Compte tenu de la grande pureté du sel de Asse (teneur en halite supérieure à 98 % de halite), la perméabilité qui est mise en jeu dans les essais de Peach est donc clairement associée à la microfissuration de l'échantillon : les discontinuités
qui seraient cicatrisées sont bien les microfissures artificiellement provoquées. La durée de
ses essais étant trop courte, il n'a pas été possible de mettre en évidence une éventuelle
restauration de la perméabilité d'un sel sain. Par conséquent, malgré la tendance observée,
il n'a pas été prouvé expérimentalement que le processus de cicatrisation sous contrainte
était susceptible d'annuler complètement la perméabilité du sel.

Or, les essais en saumure *in situ* réalisés aux MDPA mettent en jeu une
perméabilité qui n'est certainement pas de même nature que celle des essais de Peach
(microfissuration intense des carottes). En outre, d'une durée totale de 80 jours, ils n'ont
pas permis de mettre en évidence une baisse significative de la perméabilité mesurée. Cet
effet de cicatrisation, s'il existe pour un matériau salifère sain, doit avoir une constante de
temps très largement supérieure à celle de l'expérimentation *in situ* réalisée.

Quant aux conséquences induites par les changements de concentration en sel à
saturation, il semble que ceux-ci ne soient pas non plus susceptibles de modifier l'ordre de
grandeur de la perméabilité. Cette question sera examinée plus loin dans ce mémoire au
chapitre V.

1.4 Récapitulation des essais à la saumure

Nous avons pu montrer que les volumes de saumure injectés étaient relativement
bien reproduits par un modèle basé sur la loi de Darcy avec une perméabilité de 2.10<sup>-21</sup> m<sup>2</sup>
et une pression de pores initiale de 1 MPa. L'ordre de grandeur de la perméabilité d'un sel
non endommagé et de la pression de pores initiale, à plusieurs rayons de la galerie est donc
retrouvé (cf. Stormont *et al.*, 1991). Ce résultat expérimental permet également de
souligner la grande précision du système expérimental de perméamètrie mis en œuvre *in
situ*.

2. ESSAIS AU GAZ

2.1. Déplacement miscible de gaz dissous dans la phase liquide

En premier lieu, le déplacement de gaz dissous en solution est examiné. Dans ce
mode de transport, on ne peut identifier d'interfaces entre les phases fluides.

Ainsi, dans le cas du système gaz/saumure, le gaz se dissous d'abord dans la
saumure. Puis, celui-ci est transporté en solution par différents processus tels que la
convection (due au mouvement de la phase liquide), la diffusion moléculaire (provoquée
par l'agitation moléculaire) et la dispersion cinématique (phénomène associé à
l'hétérogénéité des vitesses microscopiques) (Marsily, 1986). On négligera l'adsoprtion des
éléments en solution sur la phase solide, faute de connaissance de ce phénomène dans le
contexte de l'essai dans le sel.

a- Dissolution du gaz dans la phase liquide

La quantité de gaz dissous dans la saumure peut être estimée à l'aide de la loi de
Henry (1803) :

"lorsqu'un gaz, à une température et sous pression données, est en équilibre avec
un liquide dans lequel il est soluble (sans réaction chimique mutuelle), la masse de gaz
dissoute dans un volume donné de liquide est proportionnelle à la pression " (Lévy, 1988).
Cette loi a été validée à l'aide des concepts de la thermodynamique classique pour un gaz parfait et une solution idéale. Elle s'écrit :

$$P = K X \frac{Ml}{Mg}$$

P : pression du gaz, K : constante de Henry à une température donnée, X : fraction massique de gaz dans le liquide, Ml, Mg : masses moléculaires du liquide et du gaz respectivement.

Selon cette loi, la masse dissoute de gaz (azote) est très faible : à 15 °C et à la pression atmosphérique, seulement 21 milligrammes peuvent se dissoudre dans un litre d'eau. Selon toujours la même loi, aux pression de 4 et 6 MPa, les quantités de 0,84 grammes et de 1,26 grammes sont dissoutes dans un litre d'eau. Comme le liquide mis en jeu est de la saumure (la constante de Henry pour la saumure étant inconnue), les valeurs réelles sont probablement du même ordre voire plus faibles car la dissolution diminue avec la température. On rappelle que la chambre est thermorégulée à environ 27 °C.

## b- Transport en solution du gaz dissous

Lorsque le gaz dissous se trouve en phase liquide, il peut migrer de deux manières dans le réseau poreux : par convection et par diffusion. On néglige d'emblée la dispersion cinématique vu la faible perméabilité du sel.

Par conséquent, l'équation de conservation de masse de gaz dissous en solution dans le milieu poreux s'écrit :

$$\text{div}(\underline{\underline{\underline{\text{gradC}}}} - \underset{\text{Convection}}{\overline{\overline{\overline{CV}}}}) = \frac{\partial(\phi C)}{\partial t}$$

C: concentration en gaz dissous (g/l), D: tenseur de diffusion (m²/s), V: vitesse de Darcy, (m/s), φ: porosité.

La fuite de gaz dissous dans la formation est un processus en série (dissolution puis transport) : elle est donc contrôlée par le processus le plus lent. Selon les résultats des expériences de Volckaert et al. (1995), quelques heures seulement suffisent pour dissoudre de l'hydrogène dans un matériau poreux. Par conséquent, le temps caractéristique de dissolution est probablement petit comparé aux constantes de temps de la convection et de la diffusion. En considérant un écoulement radial, ces constantes de temps et la masse de gaz dissoute injectée durant un palier, pour chacun des modes de transport, peuvent être estimées. A titre d'illustration, une application pour le palier de 6 MPa, d'une durée de 16 jours est proposée dans le tableau 9.

<table><caption>Tableau 9 : Calcul des temps caractéristiques des phénomènes relatifs à deux modes de transport en solution: la convection et la diffusion et des masses injectées calculées pour le palier de 6 MPa.</caption><thead><tr><th></th><th>Convection</th><th>Diffusion</th></tr></thead><tbody><tr><th>Temps caractéristique τ</th><td>τ<sub>c</sub> = <sup>a</sup>/<sub>V</sub> = <sup>η a<sup>2</sup> ln(b/a) φ</sup>/<sub>k(P<sub>l</sub> - P<sub>0</sub>)</sub><br>≈ 1000 ans</td><td>τ<sub>d</sub> = <sup>φ a<sup>2</sup></sup>/<sub>D</sub> = <sup>a<sup>2</sup></sup>/<sub>d<sub>o</sub></sub><br>≈ 30 jours</td></tr><tr><th>Masse injectée pour le palier de 6 MPa</th><td>M = <sup>C<sub>l</sub> k (P<sub>l</sub> - P<sub>0</sub>)</sup>/<sub>η ln(b/a)</sub> 2πh t<br>&lt; 1 mg</td><td>M = S<sub>o</sub> ∫<sub>0</sub><sup>t</sup> -D <sup>∂C</sup>/<sub>∂r</sub> |<sub>a</sub> dτ<br>numériquement &lt; 1 mg</td></tr></tbody></table>

Données : P₁=6 MPa; k= 2 10⁻²¹ m²; d₀: le coefficient de diffusion moléculaire de l'azote dans l'eau égal à 1,65 10⁻⁹ m²/s (Jost, 1960); C₁=1,26 g/l en r=a à la paroi du forage, b: rayon extérieur suffisamment loin du forage (b=100a), h: hauteur du découvert.
Un palier de pression de 6 MPa donne des valeurs de la masse injectée théorique M inférieures à 1 mg sur une durée 16 jours, correspondant à la durée du palier. Or, il a été injecté expérimentalement plusieurs dizaines de grammes (figure 7). La comparaison entre les valeurs calculées et celles mesurées suggère qu'un déplacement miscible de gaz dissous est un phénomène négligeable.

Seule une petite quantité de gaz se dissous et diffuse dans la saumure. Par conséquent, la migration de gaz mise en jeu au cours de l'expérimentation doit être considérée comme un écoulement multiphasique de fluides immiscibles.

## 2.2. Écoulement de fluides immiscibles. Rôle de la pression capillaire initiale

Supposons qu'il existe plusieurs phases fluides (non miscibles) séparées par des interfaces. Dans un milieu poreux constitué de pores pouvant atteindre une taille nanométrique tel que le sel, ces interfaces sont soumises à des forces de capillarité à l'origine d'une différence de pression. Cette différence de pression appelée la pression capillaire a été estimée précédemment dans la partie 1.3 du chapitre I : elle est comprise dans un intervalle d'une fraction de MPa à plusieurs centaines de MPa. Son influence est certaine comme peuvent l'attester les données de l'essai au gaz après percolation de saumure.

Néanmoins, on peut s'interroger sur le rôle, voire l'existence, d'un effet de la capillarité dans l'état initial du massif avant toute injection de saumure. Autrement dit, quel est le degré initial de saturation en saumure du massif ?

Deux modèles simples basés sur des hypothèses opposées peuvent contribuer à apporter des éléments de réponse :

- a- le milieu poreux saturé de gaz,

- b- Déplacement d'une frontière mobile séparant deux milieux poreux saturés.

### a- Milieu poreux saturé par le gaz

Dans un premier cas, le milieu est supposé saturé en gaz (ne pouvant réagir chimiquement avec le gaz injecté c'est à dire l'azote). Par conséquent, la phase gazeuse est supposée continue: il n'y a donc pas de pression capillaire. Dans cette phase gazeuse, deux types d'écoulement peuvent être considérés : un écoulement obéissant à la loi de Darcy et un écoulement de type Knudsen. Ce dernier a été étudié lors de la conception des essais in situ.

En outre, on admet une migration essentiellement radiale et isotherme, en géométrie cylindrique. L'historique des pressions appliquées est également pris en compte. Les résultats du modèle, en considérant les deux lois d'écoulement, et les données expérimentales sont proposées sur les figures 16 et 17. Malgré des bons ordres de grandeur, les modèles utilisés ne permettent pas de reproduire exactement les valeurs expérimentales.

Néanmoins, on notera sur la figure 16, que la courbe d'injection calculée pour le deuxième palier de 6 MPa donne des valeurs systématiquement plus élevées que celles mesurées. Au cours de la migration, la force motrice de l'écoulement semble avoir progressivement diminuée : la perméabilité a t-elle diminué ? ou alors la pression capillaire a augmenté ? l'écoulement a-t-il changé de mécanisme ? cette dernière explication nous paraît plus convaincante. En effet, on constate sur la figure 17, que l'ajustement entre le modèle basé sur une loi de Knudsen et les données est davantage satisfaisant. Il faudrait se garder d'en dégager une conclusion rapide, étant donné l'incertitude qui pèse sur les valeurs de paramètres et des modèles mis en jeu (la loi de filtration de Knudsen proposée par Schlich n'a pas été validée sur des essais de perméabilité en laboratoire).
<figure><img src="image_46.png" /></figure>

Figure 16. Masse injectée calculée pour les deux paliers gaz de 4 et 6 MPa en considérant la loi de Darcy. La valeur de la perméabilité intrinsèque donnant la meilleure concordance entre données expérimentales et modèle est de $3 \times 10^{-21} m^2$.

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

Figure 17. Masse injectée calculée pour les deux paliers gaz de 4 et 6 MPa en considérant un écoulement Knudsen. La valeur du coefficient Knudsen choisi est de $3,7 \times 10^{-7} m^2/s$ c'est à dire un ordre de grandeur comparable à celui donnée par Jockwer et Rothfuchs (1995) pour le sel de Asse.

## b- Déplacement d'une frontière mobile séparant deux milieux poreux saturés

Dans le deuxième modèle, le massif est considéré être initialement saturé en saumure à une pression de pores donnée. Les différentes migrations de gaz sont modélisées par le déplacement d'une frontière mobile séparant deux milieux poreux saturés : l'un en saumure et l'autre en gaz. Les hypothèses suivantes ont été utilisées :

i- dans la phase gazeuse, le gaz est à une pression constante, égale à la valeur de la pression de palier Pg. On néglige ainsi les pertes de charges provoquées par les effets de viscosité du fluide devant celles existant dans la phase liquide. Cette approximation reste valable dès lors que la distance de pénétration du front de gaz reste faible, de l'ordre du mètre d'après les quantités injectées et compte tenu de la porosité du sel.

ii- dans la phase liquide (en saumure), un écoulement transitoire obéissant à la loi de Darcy est considéré. Cette hypothèse est justifiée par l'interprétation des essais en saumure.
iii- l'interface entre les deux phases se déplace à une vitesse donnée par l'expression suivante :

$$ \frac{\partial R(t)}{\partial t} = - \left. \frac{k}{\eta \phi} \frac{\text{grad}(Pg - Pc)}{R(t)} \right|_{R(t)} \frac{\eta}{n} $$

R : distance du front au centre du forage; k : perméabilité intrinsèque; η : viscosité dynamique; φ : porosité; Pg : pression dans la zone en gaz; η indique la direction de l'écoulement ici considéré comme radial; Pc : pression capillaire.

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

Figure 18. Masse injectée calculée pour les deux paliers gaz de 4 et 6 MPa en considérant un modèle basé sur le déplacement d'un front de saturation. La valeur de la perméabilité intrinsèque donnant la meilleure concordance entre données expérimentales et modèle est de $2 \times 10^{-19} \text{ m}^2$.

Comme le montre la figure 18, ce modèle ne reproduit pas les données expérimentales. On obtient un ordre de grandeur comparable entre expérience et modèle pour une perméabilité de $2 \times 10^{-19} \text{ m}^2$, soit une perméabilité cent fois plus grande que la valeur calculée pour les essais en saumure ($2 \times 10^{-21} \text{ m}^2$). En utilisant une valeur de perméabilité de $10^{-21} \text{ m}^2$, les masses injectées calculées sont donc cent fois plus faibles que les valeurs mesurées.

On pourra noter que d'après l'équation de la vitesse du front, l'existence d'une pression capillaire sensiblement égale à la pression de palier va réduire considérablement le gradient de pression et donc le moteur de l'écoulement. Pour une perméabilité attendue de $10^{-21} \text{ m}^2$ qui donne des masses injectées calculées cent fois trop faibles, l'effet de la pression capillaire pour ce modèle ne va donc pas contribuer à améliorer la comparaison entre modèle et expérience. Les résultats de ce modèle du déplacement du front de saturation suggèrent que la pression capillaire, si elle est présente initialement, se caractérise certainement par une distribution non uniforme dans les microvides.

## 2.3 Récapitulation des essais au gaz

Le déplacement miscible de gaz dissous dans la phase liquide n'est pas le phénomène prédominant dans les mesures de masse injectée de gaz.

L'application de modèles simples d'écoulement de fluides immiscibles, négligeant volontairement la pression capillaire ne permet pas de reproduire les mesures réalisées *in situ*. L'écoulement de gaz observé n'obéit donc pas à une simple loi de Darcy supposant le milieu comme saturé en fluide percolant. Il semble que la pression capillaire et l'effet Knudsen aient tous les deux un rôle essentiel dont l'importance relative est très difficile à estimer. En effet, ces deux effets jouent en sens inverse : la pression capillaire à tendance à
diminuer la force motrice de l'écoulement alors que les glissements en paroi contribuent à augmenter la valeur de la perméabilité mesurée.

CONCLUSION

L'essai de perméabilité réalisé aux MDPA n'est pas la première tentative de mesure in situ de perméabilité dans le sel, mais constitue, à plusieurs égards, une expérience originale visant à mesurer dans des conditions maitrisées et aussi précises que possibles, l'infiltration de la saumure et du gaz dans un banc de sel relativement pur.

Cette expérience a montré avant tout la difficulté de mesurer des perméabilités aussi faibles, bien en dessous du seuil caractéristique des milieux réputés imperméables (≈10⁻¹⁷ m²). Dans cette condition, nous avons constaté combien est difficile la maîtrise des effets couplés à la température et au comportement viscoplastique (fluage de la paroi), susceptibles de masquer la perméabilité mesurée.

La méthodologie développée a permis de contrôler voire de minimiser les effets couplés à la température et au comportement viscoplastique (fluage de la paroi). Les données obtenues de l'essai *in situ* permettent alors d'affirmer qu'un banc de sel pur, situé assez loin d'ouvrages souterrains est perméable aussi bien au gaz qu'à la saumure.

Concernant les mécanismes d'écoulement, la fuite de saumure est relativement bien interprétée par un modèle basé sur la loi de Darcy (perméabilité intrinsèque de 2.10⁻²¹ m² et pression de pores initiale de 1 MPa), ce qui n'est pas le cas des débits de gaz mesurés. L'influence de la pression capillaire mise en évidence expérimentalement et des mécanismes complexes d'origine microscopiques ("glissements" en paroi du réseau poreux) doit être prise en compte dans la modélisation. La méconnaissance des conditions initiales (degré de saturation, pression de pores initiale, pression capillaire) et des conditions aux limites loin du forage ne nous autorise pas de proposer de manière univoque un seul mécanisme d'écoulement dans le sel.

L'origine de cette perméabilité n'est pas clairement identifiée : endommagement d'origine minière (et ceci malgré la distance importante au mur de la galerie) ? endommagement d'origine tectonique ? rôle des d'impuretés (et ceci malgré leur teneur extrêmement faible) ? Quoi qu'il en soit, la perméabilité du sel gemme au voisinage d'une structure souterraine confirmée par les essais réalisés est un phénomène qui revêt une importance majeure pour la sûreté du stockage. Il doit être considéré dans toute modélisation appropriée.
# CHAPITRE IV

## MODÉLISATION MACROSCOPIQUE
DES PROCESSUS COUPLÉS DANS LE SEL
# INTRODUCTION

Nous avons vu que le sel endommagé au voisinage d'ouvrages souterrains présente une perméabilité mesurable.

Ce phénomène peut être à l'origine de processus couplés agissant sur la sûreté d'un stockage en milieu salifère. Ces couplages sont d'autant plus forts que le sel est soumis à un gradient de température. C'est le cas par exemple d'un stockage de déchets à haute activité en formations salines.

On s'interroge alors sur les conséquences à long terme des processus couplés, sur la manière dont il conviendrait de les analyser et de les modéliser. On s'introduit alors dans le domaine du calcul de structure, qui nécessite d'abord un cadre théorique approprié. L'écriture de ce cadre constitue l'objet de ce chapitre.

# 1. PRINCIPE DE LA DÉMARCHE

## 1.1 Hypothèses Générales

Dans un premier temps, à partir du milieu salifère hétérogène et discontinu à une échelle microscopique, il s'agit de construire un milieu poreux continu équivalent. Chaque point géométrique du nouveau milieu équivalent constitue un système thermodynamique élémentaire compris dans ce qui est souvent appelé un Volume Élémentaire Représentatif (ou VER). La taille de ce VER doit être suffisamment grande pour que les différentes grandeurs physiques mises en jeu soient considérées comme des fonctions continues de la variable d'espace du milieu continu équivalent, et en même temps suffisamment petite devant la plus petite longueur caractéristique du problème étudié. Le choix de la taille du VER dépend donc du matériau et du problème posé (cf. § 1.2).

L'approche utilisée est macroscopique dans le sens où les couplages entre comportement mécanique et processus de transfert sont considérés en premier lieu à l'échelle du système élémentaire de milieu poreux. Le système thermodynamique élémentaire comprend une petite partie de squelette solide (incluant la porosité occluse) et une porosité connectée saturée par une solution saline, elle même constituée d'un seul soluté (cas idéal). Ce système est dit ouvert s'il échange de la matière avec l'extérieur, c'est à dire avec les systèmes élémentaires contigus. En définitive, le milieu poreux continu équivalent est assimilable à une superposition de trois milieux continus en interaction mécanique et chimique : le squelette, le soluté (le sel) et le solvant (l'eau). On s'efforce de prendre en considération les différents modes de transfert de chaleur (diffusion, convection) et de masse (diffusion, convection et cinétique chimique) dans un contexte suffisamment général qui souffre toutefois d'une grande restriction pour ce qui concerne la cinétique chimique : une seule réaction chimique est prise en compte entre la matrice solide et le fluide. Cette hypothèse, évidemment très discutable d'un point de vue géochimique, n'a d'autre but que de simplifier les démarches de modélisation et de calcul.

Dans un deuxième temps, on applique les principes généraux de la Thermodynamique des Milieux Poreux Ouverts (lois de conservation de masse et d'énergie et Second Principe) à un système thermodynamique élémentaire. L'essentiel du travail consiste alors à écrire l'inégalité dite de Clausius-Duhem qui permet d'identifier les variables d'état observables, et d'associer de manière adéquate les bonnes variables d'état aux potentiels thermodynamiques introduits. Ces associations auront pu être effectuées à l'aide du postulat de l'état local qui stipule que l'évolution du système élémentaire est constituée d'une succession d'états d'équilibre complètement caractérisés par les variables d'état. Moyennant une hypothèse de linéarisation physique (petites variations des variables
d'état) et géométrique (petites déformations), ces associations conduisent alors aux lois de comportement couplées du système thermodynamique élémentaire ouvert.

On souligne ici que la déformation observable étant celle du squelette (défini comme la matrice *et* l'espace poreux connecté sans fluide), on privilégiera toujours le mouvement du squelette vis à vis des autres mouvements. C'est pourquoi on s'attachera à faire apparaître les dérivés particulaires *par rapport au squelette* des différentes grandeurs introduites lorsque l'on énoncera les principes de conservation.

Comme on l'aura compris, le modèle proposé repose sur l'existence d'une perméabilité, et d'un Volume Élémentaire Représentatif nécessaire pour construire un milieu continu équivalent. Ces deux hypothèses s'appliquent à la roche salifère endommagée mais plus difficilement à la roche salifère pure en halite et saine située à grande profondeur, loin d'ouvrages souterrains et à l'abri d'efforts tectoniques intenses.

En effet, a notre avis, comme il a été discuté dans les chapitres I et III, il n'existe pas de preuves expérimentales convaincantes de l'existence ou non d'une porosité connectée dans une roche salifère pure en halite, *saine* et loin d'un ouvrage souterrain :

a- la perméabilité du sel sain, si elle existe, n'est pas mesurable avec les dispositifs de perméatrie *in situ* actuels (la limite en perméabilité intrinsèque étant $10^{-21} \text{ m}^2$);

b- nos mesures aux MDPA montent qu'il existe une perméabilité pour un sel relativement pur et située à 4 rayons de la galerie la plus proche. L'origine de cette perméabilité n'a pas été clairement identifiée (endommagement d'origine minière ou tectonique ? rôle des impuretés en teneur extrêmement faible ?)

Une porosité et une perméabilité même extrêmement faibles et non mesurables au laboratoire pourraient néanmoins avoir un rôle certain sur de grandes échelles de temps. Dans le cadre d'une étude de sûreté à très long terme d'un stockage en souterrain, cette éventualité ne peut alors être complètement écartée et mérite d'être étudiée. Par conséquent, faute de preuves suffisantes, on prendra comme hypothèse que le sel sain possède une faible porosité cinématique, tout en gardant à l'esprit le caractère pessimiste et discutable de cette assertion.

## 1.2. Estimation de la taille du VER

Afin de justifier l'approche du milieu poreux équivalent, un ordre de grandeur de la taille du VER est nécessaire. On s'intéressera plus particulièrement à un sel endommagé dont la perméabilité est prouvée expérimentalement. La taille du VER est estimée à l'aide d'un logiciel réalisé par F. Coste au G.3S pour les milieux fracturés dont on rappelle brièvement les différentes étapes de calcul.

Dans une première étape, un réseau de microfissures est généré en deux dimensions de manière aléatoire sur un domaine donné. Par conséquent, les lois statistiques (loi uniforme, exponentielle, normale, ou lognormale pour le logiciel) auxquelles obéissent la longueur, l'épaisseur et l'orientation des microfissures doivent être connues. Dans notre cas, nous avons peu de données statistiques concernant l'orientation : celle-ci est donc supposée uniforme. Les travaux de Lebrun (1993) sur le sel des MDPA nous autorisent à prendre comme longueur de microfissures une valeur comprise entre le millimètre et le centimètre et comme épaisseur moyenne un micron au maximum. Suivant les travaux du même auteur, la densité de microfissure et estimée à $10^5$ microfissures/m$^2$, soit 10 microfissures dans un carré de 10 millimètres de côté.

A un deuxième stade, un gradient hydraulique est appliqué au domaine microfissuré. Si le réseau généré est percolant, une perméabilité correspondante à l'échelle du domaine est obtenue grâce au débit calculé à la frontière du même domaine. Le débit calculé à l'intérieur de chaque microfissure obéit à la loi de Poiseuille.
Il existe une taille critique de domaine à partir de laquelle la perméabilité ne varie plus avec la taille du domaine : cette taille critique est sensiblement celle du VER.

Des calculs ont été effectués avec la densité proposée et une longueur moyenne de 1 cm et une épaisseur de 1 micron (figure 1). On souligne ici que ce dernier paramètre n'influence pas la taille du VER. Afin d'apprécier la sensibilité du paramètre longueur, deux lois statistiques ont été utilisées : loi uniforme et loi normale en faisant varier le rapport moyenne/écart-type (égal à 1 et 0,4 successivement). Les résultats de la figure 1 montrent que :

a- l'ordre de grandeur de la perméabilité d'un sel peu endommagé est retrouvé (de l'ordre de $10^{-21} \text{ m}^2$).

b- la taille du VER est de l'ordre de plusieurs cm (entre 5 et 10 cm). Par conséquent, pour des structures souterraines dont la taille caractéristique est de plusieurs dizaines de mètres, l'approche milieu poreux équivalent est satisfaisante.

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

Figure 1. Évolution de la perméabilité calculée en fonction de la taille du domaine de roche microfissurée.

Afin d'étudier la sensibilité de la longueur de microfissures, plusieurs lois statistiques ont été considérées. La grandeur $l/\sigma$ désigne le rapport de la moyenne (fixée à 1 cm) par l'écart-type.

L'épaisseur a été prise égale à 1 micron. La densité est de $10^5 \text{ microfissures/m}^2$.

Soulignons une nouvelle fois que ces calculs simples ne visent qu'à obtenir un *ordre de grandeur* de la taille du VER. En effet, nous avons conscience que :

a- le problème de la détermination de la connectivité d'un milieu considéré comme microfissuré (aux joints de grains) est fondamentalement à trois dimensions (3D);

b- les lois statistiques utilisées dans les calculs n'ont pas été validées par des observations *in situ*. Par exemple, selon Davy (1994) et Bour (1996), les distributions de la longueur de fissures obéissent plutôt à une loi de puissance (la connectivité du réseau fissuré étant alors très sensible à la puissance de l'exposant de la loi).

## 2 ÉNONCÉ DES PRINCIPES FONDAMENTAUX DE LA THERMODYNAMIQUE DES MILIEUX POREUX OUVERTS

### 2.1. Conservation de la masse

L'application des principes fondamentaux de la thermodynamique des milieux poreux ouverts qui suit, est proche dans sa formulation à celle de Coussy (1991). C'est
pourquoi on ne détaillera pas les calculs tout en invitant le lecteur intéressé à consulter l'ouvrage du même auteur.

Afin de garder un caractère général à cette approche, considérons deux constituants en solution, notés 1 et 2 pouvant réagir chimiquement suivant l'équation :

$$I \leftrightarrow 2$$

Pour chacun des constituants 1 et 2, l'équation de conservation de masse de soluté dans la phase liquide s'écrit :

$$\dot{m}_i = - \text{Div} \underline{M}_i + P_{i \leftrightarrow j} \quad (1)$$

$\underline{M}_i$ : vecteur courant de masse de soluté i dans la phase fluide, $m_i$ : apport de masse de soluté i dans la phase fluide, $dm_i/dt = \dot{m}_i$ : dérivée particulière de $m_i$ par rapport au squelette, Div : l'opérateur divergence, $P_{i \leftrightarrow j}$ : taux de production (ou de disparition) de soluté j à partir du soluté i, dû à la réaction chimique (kg.m³.s⁻¹).

L'équation de conservation de la masse pour le solvant est obtenue de la même façon :

$$\dot{m}_3 = - \text{Div} \underline{M}_3 \quad (2)$$

$\underline{M}_3$ : vecteur courant de masse de solvant dans la phase fluide, $m_3$ : apport de masse de solvant dans la phase fluide, $dm_3/dt = \dot{m}_3$ : dérivée particulière de $m_3$ par rapport au squelette.

Dans le cas d'une interaction géochimique dans le système saumure-halite, un des constituants, celui indicé 1 par exemple, provient de la matrice. L'équation de la conservation de la masse se simplifie à :

$$\dot{m}_1 = P_{1 \leftrightarrow 2} \quad (3)$$

avec :

$$P_{1 \leftrightarrow 2} = \dot{\xi} = \chi \quad (4)$$

$\xi$ : le degré d'avancement volumique de la réaction et $\chi$ : la vitesse d'avancement de la réaction.

## 2.2 Conservation de la quantité de mouvement. Théorème de l'énergie cinétique.

On considère un élément de volume $d\Omega$ du milieu poreux constituant le système élémentaire associé aux particules suivantes : le squelette, les constituants en solution et le solvant. En prenant comme vitesse virtuelle la vitesse du squelette solide et dans le cadre de l'hypothèse des petites perturbations, le Principe des Travaux Virtuels permet d'écrire la conservation de l'énergie cinétique sous la forme suivante en description lagrangienne liée à la vitesse du squelette solide :

$$\frac{DK}{Dt} + \underbrace{\int \Pi : \dot{\underline{\underline{\Delta}}} d\Omega}_{\text{Puissance des efforts internes}} - \sum_{i=1}^{3} \underbrace{\int M_i^t P \underline{\underline{\underline{\gamma}}}_i}_{\text{Puissance des efforts d'inertie}} = \underbrace{\int \underline{\underline{V}} r F d\Omega}_{\text{Puissance des efforts externes}} + \underbrace{\int \underline{\underline{T}} \underline{\underline{V}} da}_{\text{Puissance des efforts externes}} \quad (5)$$

$\underline{\underline{\Pi}}$: tenseur des contraintes de Piola-Kirchhoff, $\dot{\underline{\underline{\Delta}}}$: dérivé particulière dans le mouvement du squelette du tenseur de Green-Lagrange, : produit tensoriel; $\underline{\underline{P}}$: gradient de transformation $P_{ij}=\partial X_i / \partial x_j$; $X_i$: coordonnées du squelette dans la configuration initiale de $\Omega$ ; $x_j$: coordonnées du squelette dans la configuration actuelle de $\Omega$ ; $\underline{\underline{\gamma}}_i$: vecteur accélération attaché à la particule i ; $\underline{\underline{V}}$: vitesse du squelette; $\underline{\underline{F}}$: densité des forces de volume (la gravité), r : densité du milieu salifère, T : densité des efforts de surface.
Pour alléger les calculs, la puissance des efforts d'inertie ne sera pas réécrite ultérieurement. Sa contribution sera d'ailleurs négligée dans le cadre de l'hypothèse de linéarisation physique adoptée dans la suite (§2).

DK/Dt signifie la dérivée matérielle de l'énergie cinétique qui se définit par :

$$ \frac{DK}{Dt} = \frac{D}{Dt} \int_{\Omega} k d\Omega = \frac{d}{dt} \int_{\Omega} (1-\phi) \rho^{sk} \underline{V}^2 d\Omega + \sum_{i=1}^{3} \frac{d^i}{dt} \int_{\Omega} \phi_i \rho_i \underline{V}_i^2 d\Omega \quad (6) $$

k : la densité d'énergie cinétique, $\underline{V}$ : vitesse du squelette; $\underline{V}_i$ : vitesse du constituant i; $\phi$ : la porosité totale, $\phi_i$ : la fraction volumique du constituant i avec $\phi=\phi_1+\phi_2+\phi_3$; $d/dt$ la dérivée particulaire par rapport au squelette, $d^i/dt$ la dérivée particulaire par rapport au constituant i.

## 2.3 Conservation de l'énergie

La conservation de l'énergie du système élémentaire est donnée par :

$$ \frac{DE}{Dt} = \frac{DE}{Dt} + \frac{DK}{Dt} = \int_{\Omega} \underline{V} r F d\Omega + \int_{\partial\Omega} \underline{T} \underline{V} ds + \sum_{i=1}^{3} \int_{\partial\Omega} -p v_i \underline{M}_i \underline{n} ds + \int_{\Omega} R d\Omega - \int_{\partial\Omega} \underline{q} \underline{n} ds $$

E : énergie totale du système élémentaire, E : énergie interne, K : énergie cinétique, p : pression du fluide, $\underline{v}_i$ : volume massique partiel du constituant i, $\underline{M}_i$ : vecteur courant de masse associée au constituant i, R : chaleur interne générée par unité de volume; q : flux de chaleur; n : normale extérieure au domaine $\Omega$.

En introduisant la densité d'énergie interne e et la densité d'enthalpie massique du constituant i ($\underline{h}_i$), la conservation de l'énergie peut être réécrite à l'aide du théorème de l'énergie cinétique (5) :

$$ \dot{e} = \underline{\underline{\Pi}}: \dot{\underline{\underline{\Delta}}} + \sum_{i=1}^{3} [\underline{\underline{M}}_i \underline{\underline{t}} \underline{\underline{P}} \underline{\underline{F}} - \text{Div}(\underline{h}_i \underline{\underline{M}}_i)] + R - \text{Div} \underline{q} \quad (7) $$

P : gradient de transformation $\underline{P}_{ij}=\partial x_i / \partial x_j$; X$_i$ : coordonées du squelette dans la configuration initiale; x$_j$ : coordonnées dans la configuration actuelle.

## 2.4 Le Second Principe de la Thermodynamique. L'Inégalité de Clausius-Duhem

Le Second Principe de la Thermodynamique pour le système élémentaire considéré s'écrit :

$$ \frac{DS}{Dt} \geq \int_{\Omega} \frac{R}{T} d\Omega - \int_{\partial\Omega} \frac{q}{T} \underline{n} ds \quad (8) $$

S : l'entropie du système élémentaire d$\Omega$, T : température absolue.

En considérant la densité d'entropie s du système élémentaire et la densité massique d'entropie du constituant i ($\underline{s}_i$), l'inégalité (8) peut être réécrite localement après multiplication par la température T :

$$ T \dot{s} + T \sum_{i=1}^{2} \text{Div}(\underline{s}_i \underline{M}_i) - \frac{q}{T} \underline{\underline{\text{Grad}}} T + R - \text{Div} \underline{q} \geq 0 \quad (9) $$

L'expression de R-Divq déduite de (7) est alors introduite dans (9). En considérant la densité d'énergie libre $\Phi=e-Ts$, on obtient alors l'inégalité dite de "Clausius-Duhem" :
$$ \underline{\underline{\Pi: \dot{\underline{\underline{\Delta}}} - s\dot{T} - \dot{\Phi} + T \sum_{i=1}^{3} \mathrm{Div}(s_i \underline{\underline{M_i}}) + \sum_{i=1}^{3} [\underline{\underline{M_i}}' \underline{\underline{P}} \underline{\underline{F}} - \mathrm{Div}(h_i \underline{\underline{M_i}})] - \frac{q}{T} \underline{\underline{\mathrm{Grad}}}} T \geq 0} $$

En écrivant le potentiel chimique pour chacun des constituants : $\mu_i=h_i-Ts_i$ et en remplaçant $\mathrm{Div}\underline{\underline{M_i}}$ par son expression dans les équations de conservation de la masse (1) et (2), on peut reformuler l'inégalité de Clausius-Duhem de la manière suivante :

$$ \underline{\underline{\Pi: \dot{\underline{\underline{\Delta}}} - s\dot{T} - \dot{\Phi} + \sum_{i=1}^{3} \mu_i \dot{m}_i + \frac{P(\mu_2 - \mu_1)}{\Psi_2 \mathrm{Dissipation \ Chimique}} - \frac{q}{T} \mathrm{Grad} T}} \\ \underline{\underline{\Psi_1 \mathrm{Dissipation \ Intrinsèque} + \sum_{i=1}^{3} (-\underline{\underline{M_i}} \mathrm{Grad} \mu_i - s_i \underline{\underline{M_i}}' \underline{\underline{P}} \mathrm{Grad} T + \underline{\underline{M_i}} \underline{\underline{F}})} \geq 0}} \\ \underline{\underline{\Psi_4: \mathrm{Dissipation associée au transport de masse}}} $$

ou encore

$$ \Psi_1 + \Psi_2 + \Psi_3 + \Psi_4 \geq 0 \qquad (10) $$

avec $\Psi_1$: dissipation intrinsèque (car indépendante des gradients spatiaux); $\Psi_2$: dissipation chimique; $\Psi_3$: dissipation thermique; $\Psi_4$: dissipation associée au transport de masse.

A partir de la relation (10), l'écriture de l'inégalité de Clausius-Duhem est retrouvée dans les cas particuliers suivants :

a- le comportement thermo-mécanique d'un matériau monophasique (voir par exemple, Lemaitre et Chaboche, 1988), en faisant $\underline{\underline{M_i}}=0$ et $m_i=0$ dans la relation (10) :

$$ \underline{\underline{\Pi: \dot{\underline{\underline{\Delta}}} - s\dot{T} - \dot{\Phi} - \frac{q}{T} \mathrm{Grad} T \geq 0}} $$

b- le comportement thermo-mécanique d'un matériau poreux saturé par une seule phase, non réactive avec le squelette (voir par exemple, Coussy, 1991) :

$$ \underline{\underline{\Pi: \dot{\underline{\underline{\Delta}}} - s\dot{T} - \dot{\Phi} + \dot{m} - \frac{q}{T} \mathrm{Grad} T + \underline{\underline{M}} \mathrm{Grad} g_m - s_m \underline{\underline{M}} \mathrm{Grad} T + \underline{\underline{MF}} \geq 0}} $$

La nature réactive du liquide saturant fait donc intervenir une dissipation supplémentaire de nature chimique :

$$ \Psi_2 = P_{I \leftarrow 2}(\mu_2 - \mu_1) = A\chi $$

En ayant posé : $P_{I \leftarrow 2} = \chi$ et $A = \mu_2 - \mu_1$

où $\chi$ et $A$ sont des grandeurs respectivement équivalentes à une vitesse volumique d'avancement et à une affinité de la réaction chimique. A l'équilibre chimique, on a bien sûr l'égalité entre les potentiels chimiques du réactif 1 et du produit 2 :

$$ \mu_2 = \mu_1 \Rightarrow \Psi_2 = 0 $$

Dans le cas particulier d'une intéraction géochimique fluide-roche, le système élémentaire de milieu poreux est dit fermé vis-à-vis du constituant 1 attaché à la matrice, ainsi (Cosenza et Ghoreychi, 1993) :
$$ \mathrm{Div} \underline{M}_1 = 0 \Rightarrow \dot{m}_1 = \xi $$

$\xi$ : variable traduisant le degré d'avancement volumique de la réaction.

## 2.5 Formulation générale des lois de comportement

### a- Les équations d'état

Le postulat de l'état local stipule que l'état d'énergie des systèmes élémentaires constituant le milieu poreux équivalent est complètement caractérisé à un instant donné par la donnée des variables d'état des systèmes élémentaires. Les gradients spatiaux et temporels n'intervenant pas pour définir cet état, ce postulat implique donc que les évolutions des systèmes élémentaires constituent des successions d'états d'équilibre.

En outre, les gradients spatiaux (GradT et Gradμᵢ) étant quelconques, ils sont donc pris égaux à 0. A partir du postulat de l'état local et de cette dernière hypothèse, l'inégalité (10) se réduit à :

$$ \Psi_I = \underline{\underline{\Pi}} : \underline{\underline{\Delta}} - s \dot{T} - \dot{\Phi} + \sum_{I=1}^{3} \mu_i \dot{m}_i + A \chi \geq 0 \quad (11) $$

La dissipation $\Psi_I$ est dite "intrinsèque" car elle ne dépend pas des gradients spatiaux (GradT et Gradμᵢ). Ainsi, cette dissipation est relative au seul système élémentaire étudié et non aux autres systèmes élémentaires qui lui sont contigus.

D'après le postulat de l'état local, les potentiels thermodynamiques notamment la densité d'énergie libre et les potentiels chimiques ne dépendent que des variables d'état. Ici, les variables d'état choisies sont le tenseur des Green-Lagrange $\underline{\underline{\Delta}}$, la température T, les variables mᵢ décrivant les apports de masse des constituants en solution et une variable décrivant l'apport de masse de solvant par unité de volume initial notée m₃. Par conséquent, on peut écrire :

$$ \Phi = \Phi(\underline{\underline{\Delta}}, T, m_1, m_2, m_3, \dots, \underline{\underline{\alpha}}) \quad (12) $$

avec $\underline{\underline{\alpha}}$ une variable interne scalaire ou tensorielle traduisant généralement l'état interne de la matière (microstructure cristalline, configuration de microfissures et cavités par exemple). Dans le cas particulier d'un milieu poreux, la variable interne peut traduire la variation de porosité causée par les déformations plastiques du matériau ou par les phénomènes de dissolution-cristallisation.

La dérivée particulière de la densité d'énergie libre dans le mouvement du squelette est alors donnée par :

$$ \dot{\Phi} = \frac{\partial \Phi}{\partial \underline{\underline{\Delta}}} : \underline{\underline{\Delta}} + \frac{\partial \Phi}{\partial T} \dot{T} + \sum_{i=1}^{3} \frac{\partial \Phi}{\partial m_i} \dot{m}_i + \frac{\partial \Phi}{\partial \underline{\underline{\alpha}}} : \dot{\underline{\underline{\alpha}}} \quad (13) $$

En substituant cette équation dans l'inégalité (11), on obtient :

$$ (\dot{\underline{\underline{\Pi}}} - \frac{\partial \Phi}{\partial \underline{\underline{\Delta}}}) : \underline{\underline{\Delta}} - (S + \frac{\partial \Phi}{\partial T}) \dot{T} + \sum_{i=1}^{3} (\mu_i - \frac{\partial \Phi}{\partial m_i}) \dot{m}_i - \frac{\partial \Phi}{\partial \underline{\underline{\alpha}}} : \dot{\underline{\underline{\alpha}}} + A \chi \geq 0 \quad (14) $$

Une hypothèse supplémentaire de découplage sur les dissipations $\Psi_1 + \Psi_3 + \Psi_4 \geq 0$ et $\Psi_2 \geq 0$ conduit de la même façon aux inégalités suivantes :
$$ (\underline{\underline{\dot{\Pi}}}-\frac{\partial\Phi}{\partial\underline{\underline{\Delta}}}):\dot{\underline{\underline{\Delta}}}-(S+\frac{\partial\Phi}{\partial T})\dot{T}+\sum_{i=1}^{3}(\mu_{i}-\frac{\partial\Phi}{\partial m_{i}})\dot{m}_{i}-\frac{\partial\Phi}{\partial\underline{\underline{\alpha}}}:\dot{\underline{\underline{\alpha}}}\geq 0 \quad (15a) $$

et $A\chi \geq 0$

avec également : \hfill (15b)

$$ -\frac{\partial\Phi}{\partial\underline{\underline{\alpha}}} : \dot{\underline{\underline{\alpha}}} \geq 0 \qquad (16) $$

On introduit alors la définition d'un paramétrage normal relativement à une variable d'état donnée par Coussy (1991) :

> "Un paramétrage est dit normal relativement à une variable d'état du système élémentaire si, dans les évolutions réelles du système élémentaire, on peut considérer des variations de la variable d'état en question indépendamment des variations des autres variables d'état du paramétrage".

Le paramétrage normal relativement à la température T consiste donc à supposer qu'il existe une variation dT indépendante des autres variations $d\underline{\underline{\Delta}}$ $dm_i$ (i=1,2,3) et $d\underline{\underline{\alpha}}$ qui peuvent être alors prises égales à 0. L'inégalité (15a) devant être, dans ces conditions, toujours vérifiée et en admettant que l'entropie S est indépendante de la vitesse $\dot{T}$, on a alors l'équation d'état :

$$ S = -\frac{\partial\Phi}{\partial T} \qquad (17) $$

En supposant d'une part un paramétrage normal pour les autres variables d'état $\underline{\underline{\Delta}}$, $m_i$ et $\underline{\underline{\alpha}}$ et d'autre part une indépendance des grandeurs actuelles $\underline{\underline{\Pi}}$ et $\mu_i$ vis-à-vis des vitesses $\dot{\underline{\underline{\Delta}}}$, $\dot{m}_i$ et $\dot{\underline{\underline{\alpha}}}$, on obtient de la même façon les équations d'état suivantes :

$$ \underline{\underline{\Pi}} = \frac{\partial\Phi}{\partial\underline{\underline{\Delta}}} \qquad (18); \qquad \mu_i = \frac{\partial\Phi}{\partial m_i} \quad (i=1,2,3) \qquad (19a,b,c) $$

Par conséquent, on associe aux cinq variables d'état $\underline{\underline{\Delta}}$, T, $m_i$ (i=1,2,3), les cinq "forces thermodynamiques" $\underline{\underline{\Pi}}$, T et $\mu_i$ (i=1,2,3) via les équations d'état (17), (18) et (19a,b,c).

## b- les lois de conduction

Outre l'inégalité (11), le second principe de la thermodynamique et le postulat d'état local permettent d'écrire :

$$ -\frac{q}{T} \underline{\underline{\underline{Grad}} T} + \sum_{i=1}^{3} \left( -\underline{\underline{\underline{M}}}_i \underline{\underline{\underline{Grad}} \mu}}_i - s_i \underline{\underline{\underline{M}}}_i \underline{\underline{\underline{Grad}} T}} + \underline{\underline{\underline{M}}}_i \underline{\underline{\underline{F}}} \right) \geq 0 \qquad (20) $$

Les lois de comportement supplémentaires, relatives aux processus de transfert (de chaleur et de masse) doivent satisfaire cette inégalité, afin d'obéir aux principes généraux de la thermodynamique. Ces lois portant sur les variables flux $q$ et $M_i$ peuvent être formulées à l'aide d'une approche purement phénoménologique. Une autre voie consiste à supposer une forme mathématique admissible avec l'inégalité (20) et à la vérifier a posteriori expérimentalement. Suivant cette seconde démarche, l'hypothèse du "mécanisme dissipatif normal" (Coussy, 1991) est alors utilisée :
Soit la forme de dissipation non négative suivante :

$$ \Psi = \sum_{i} Q_{i} \text{Grad}X_{i} $$

Au flux $Q_i$, on associe la force thermodynamique $\text{Grad}X_i$.

Soit $D(Q_1, Q_2, ..., Q_n)$ une fonction non négative et convexe de ses arguments telle que $D(0,0,...,0)=0$. L'hypothèse de normalité du mécanisme de dissipation associé aux flux $Q_i$ revient à stipuler l'existence d'un paramétrage et d'une telle fonction, telles que les relations complémentaires recherchées s'écrivent :

$$ \text{Grad}X_i = \frac{\partial D}{\partial Q_i} $$

La fonction $D$ reçoit le nom de potentiel de dissipation.

Ainsi, en supposant un mécanisme dissipatif normal pour les dissipations $\Psi_3$ et $\Psi_4$, les forces thermodynamiques associées aux flux $q/T$ and $M_i$ sont données par :

$$ - \underline{\text{Grad}T} = \frac{\partial D}{\partial q} ; \quad - [\underline{\text{Grad}}\mu_i + s_i \underline{\text{Grad}}T - \underline{F}] = \frac{\partial D}{\partial(M_i)} \quad (31) $$

où $D(q/T, M_1, M_2)$ est le potentiel de dissipation lié à la dissipation $\Psi_3 + \Psi_4$. On peut prendre une forme quadratique, la plus simple, pour exprimer le potentiel $D$ :

$$ D = \frac{1}{2} \sum_{i=1}^{2} \sum_{j=1}^{2} t \underline{M}_i k_{ij}^{-1} \underline{M}_j + \sum_{i=1}^{2} t \underline{M}_i C_i \frac{q}{T} + \frac{1}{2} \frac{t q}{T} \underline{\lambda}^{-1} q \quad (32) $$

avec $\underline{k}, \underline{C}_i$ et $\underline{\lambda}$ des tenseurs symétriques.

Si l'on néglige les effets couplés ($C^k = 0$) et ($k_{ij}=0$ $i \neq j$), on retrouve les formes plus classiques :

$$ q = -\underline{\lambda} \underline{\text{Grad}}T \quad (33) ; \quad \underline{M}_i - \underline{k}_{ii} [(\underline{\text{Grad}}\mu_i + s_i \underline{\text{Grad}}T - \underline{F})] (i=1,2) \quad (34a,b) $$

On reconnaît l'équation (33), la loi de Fourier et l'équation (34), la loi de Fick (avec un terme de transport supplémentaire en $\underline{\text{Grad}}T$). Si les effets couplés ne sont pas négligés, on obtient :

$$ \underline{q} = -\underline{\lambda} \underline{\text{Grad}}T - \underline{\lambda} : \sum_{i=1}^{2} \underline{C}_i^i \underline{M}_i \quad (35); $$

$$ \underline{M}_i = -\underline{k}_{ii} [(\underline{\text{Grad}}\mu_i + s_i \underline{\text{Grad}}T - \underline{F})] - \underline{k}_{ii} : \underline{C}_i^i \frac{\underline{q}}{T} \quad (36) $$

De cette manière, le flux de chaleur est fonction des flux de transport de masse et réciproquement. Un résultat équivalent, établi par Hassanizadeh (1986) est obtenu en utilisant le formalisme de la Thermodynamique des Processus Irréversibles (TPI), étendu aux milieux poreux. Cette approche diffère au moins sur deux points de celle qui a été utilisée dans ce chapitre :

a- Hassanizadeh suit deux étapes préliminaires : écriture des lois de conservation (masse, moment, énergie et entropie) à l'échelle microscopique puis par une technique d'homogénéisation, obtention des lois de conservation et des lois de comportement de chacune des phases à l'échelle macroscopique du VER ;
b- dans le formalisme du même auteur, d'un point de vue mécanique, la déformation de la matrice solide est négligée.

Le même raisonnement avec l'inégalité (15b) permet d'écrire une loi de cinétique chimique du premier ordre reliant la vitesse volumique d'avancement de la réaction à l'affinité de la réaction :

$$ \chi = k_c A \qquad (37) $$

$k_c$ : constante de cinétique chimique.

# 3. CAS PARTICULIER D'UNE LINÉARISATION PHYSIQUE

## 3.1 Hypothèses

### a. Un modèle macroscopique et linéaire

En premier lieu, on rappelle que l'approche utilisée est macroscopique : à cause de la perméabilité du milieu, les transferts en solution peuvent se produire sur des distances beaucoup plus grandes que le grain. Le système est donc ouvert. Néanmoins comme on le verra, les transports de masse en deçà du volume du système élémentaire, contrôlant les variables d'état macroscopique tels que la déformation ou la porosité peuvent être aussi pris en compte grâce à une approche micro-macro adéquate.

En second lieu, le modèle envisagé est linéaire car il repose sur deux linéarisations. La première est géométrique : on se place dans l'hypothèse des petites déformations. La géométrie du milieu considéré reste donc inchangée au cours de l'évolution. Cette linéarisation nous autorise à remplacer d'une part le tenseur des contraintes de Piola-Kirchhoff $\Pi$ par le tenseur des contraintes de Cauchy $\sigma$ et d'autre part, le tenseur de Green-Lagrange $\Delta$ par le tenseur des déformations élastiques $\varepsilon$ (pas de variable interne, $\sigma=0$). La déformation anélastique, symbolisée par le tenseur $\varepsilon^{inel}$ est introduite ultérieurement en considérant l'hypothèse de la partition des déformations :

$$ \underline{\underline{\varepsilon}}^t = \underline{\underline{\varepsilon}} + \underline{\underline{\varepsilon}}^{inel} $$

où $\underline{\underline{\varepsilon}}^t$ est la déformation totale et $\underline{\underline{\varepsilon}}$ est la déformation élastique.

La deuxième linéarisation est physique : seules les petites variations (i.e. au premier ordre) des variables d'état sont prises en compte. Pour alléger les écritures, on considère des propriétés physiques isotropes.

### b- Choix des variables d'état

On prend comme variables d'état : $\underline{\underline{\varepsilon}}$, la déformation élastique, T, la température, m, la masse de saumure (par unité de volume initial de milieu poreux) apportée au système élémentaire (on regroupe de cette manière, les variables apports de masse soluté $m_2$ et de solvant $m_3$), et $\underline{\underline{\xi}}$, la masse de sel (par unité de volume initial de milieu poreux) donnée uniquement par la réaction chimique de dissolution (ou cristallisation) entre la saumure et la matrice salifère à l'échelle du système élémentaire. Par conséquent, à la suite de la démarche du § 2, les lois de comportement à considérer sont :

$$ \text{a- la loi de comportement mécanique : } \underline{\underline{\sigma}} = \frac{\partial \Phi}{\partial \underline{\underline{\varepsilon}}} \qquad (38) $$
$$
b- \text{la loi de comportement thermique : } S = -\frac{\partial \Phi}{\partial T} \qquad (39)
$$

$$
c- \text{ la loi de comportement hydraulique (attachée à la saumure) : } g_m = \frac{\partial \Phi}{\partial m} \quad (40)
$$

$g_m$ étant l'enthalpie massique de la saumure

$$
d- \text{ la loi de comportement associée à la cinétique chimique : } A = \frac{\partial \Phi}{\partial \xi} \qquad (41)
$$

$$
\text{avec le paramétrage de la densité d'énergie libre : } \qquad \Phi = \Phi(\underline{\epsilon}, T, m, \xi) \tag{42}
$$

Ce choix des variables d'état permet ainsi d'intégrer dans le formalisme les deux cas
limites: (1) le cas où n'existe pas de réaction chimique (ξ=0) et (2) le cas où le système
élémentaire est fermé vis à vis du transport de saumure (perméabilité nulle, m=ξ).

- Principe d'obtention des lois de comportement

Dans le cadre d'une modélisation linéaire, il existe a priori deux chemins pour parvenir à l'obtention des lois de comportement. Le premier consiste à prendre une forme quadratique des variables d'état pour la densité d'énergie libre Φ puis la dériver pour obtenir les lois recherchées. Ou au contraire, à partir de considérations physiques simples, les lois sont formulées puis intégrées afin d'obtenir la forme de Φ qui obéisse aux principes de la thermodynamique. C'est cette deuxième approche, plus physique que la première, que nous avons choisie. Dans une première étape, il s'agit donc de déterminer le tenseur des contraintes σ, l'entropie du système S, l'enthalpie masique gm et l'affinité chimique A.

3.2 Loi de comportement mécanique

Dans cette partie, seules les déformations élastiques sont considérées. Les déformations irréversibles sont introduites ultérieurement (§ 4.6)

Considérons le volume d'un système élémentaire Ω dans lequel se produit une
interaction chimique saumure-roche. Dans ce système, le volume de la phase solide est
noté V_s et le volume de la phase liquide, V_{fl}. Une réaction de dissolution ou de
cristallisation à l'échelle du système est susceptible de produire une déformation volumique
du système, consécutive au changement d'état (solide ⇌ liquide) de la matière dissoute ou
cristallisée.

Pour illustrer notre propos, soit dm_s une quantité infinitésimale de masse de sel (comptée positive) issue de la réaction chimique, une dissolution par exemple, dans le système Ω. La réaction amène le système Ω d'un état 1 à un état 2. La variation du volume total au cours de la transformation, dV s'écrit :

$$
dV = (V_2^s + V_2^{fl}) - (V_1^s + V_1^{fl}) = dV^s + dV^{fl} \quad (43)
$$

$dV^s$: variation de volume de la phase solide; $dV^{fl}$: variation de volume de la phase liquide; $V_i^s$: volume occupé par la phase solide à l'état i (i=1,2); $V_i^{fl}$: volume occupé par la phase liquide à l'état i (i=1,2).

La variation de volume de la phase solide, dV<sup>s</sup> est donnée par :
$$dV^s = - \frac{dm_s}{\rho^s}$$

ρs : mase volumique du solide.

La variation de volume de la phase liquide, dVfl :

$$dV^{fl} = m_{l}^{fl} \left( \frac{1}{\rho_{2}^{fl}} - \frac{1}{\rho_{l}^{fl}} \right) + \frac{dm_{s}}{\rho_{2}^{fl}}$$

ρfl : masse volumique du liquide à l'état i (i=1,2); mfl : masse de liquide dans le système à l'état 1.

Le premier terme à droite de l'égalité traduit la variation de masse volumique du liquide au cours de la réaction. Cette variation due essentiellement à la réaction chimique est très faible. En effet, en utilisant les tables de l'ATG (1986), on constate qu'une variation de 10 % de la concentration massique induit un changement de 1 % de la masse volumique de saumure. On négligera donc cet effet dans la suite du chapitre.

Par conséquent, la variation relative du volume du système vérifie la relation suivante :

$$\frac{dV}{V} = \frac{dV^s}{V} + \frac{dV^{fl}}{V} = \left[ \frac{1}{\rho^{fl}} - \frac{1}{\rho^s} \right] \xi \quad (44)$$

ξ : masse dm<sub>s</sub> divisée par le volume du système ; ρ<sub>l</sub><sup>fl</sup> = ρ<sub>2</sub><sup>fl</sup>.

La grandeur ξ est comptée positive (resp. négative) lorsque la masse dm<sub>s</sub> est ajoutée (resp. retirée) à celle de la saumure. Comme ρ<sup>s</sup>>ρ<sup>fl</sup>, on a alors 1/ρ<sup>s</sup> < 1/ρ<sup>fl</sup>. Ainsi, si la dissolution (resp. cristallisation) est prédominante à l'échelle de Ω, le volume du système V aura tendance à croître (resp. décroître).

Par conséquent, en extrapolant ce raisonnement à un système élémentaire de roche naturelle, on suppose qu'il existe une déformation volumique liée à la réaction chimique s'écrivant :

$$\varepsilon^{ch} = \alpha^{c} \xi \quad (45)$$

α<sup>c</sup> : coefficient équivalent à une dilatation chimique.

Dans notre exemple très simple, le terme α<sup>c</sup> est égal à la différence des inverses de masse volumique des deux phases. Néanmoins, on prendra soin de ne pas généraliser l'expression de ce coefficient à un milieu naturel tel qu'une roche salifère dont la teneur en halite est variable : la roche est certainement le siège de multiples réactions chimiques, le bilan des déformations volumiques des différentes phases est complexe.

En outre, on peut se demander si cette déformation d'origine chimique ne correspondrait pas à la déformation transitoire observée sur les essais de fluage lorsque l'échantillon est mis en contact avec de la saumure (Spiers et al., 1986 ; Brodsky et Munson, 1991; Chen, 1994) ?

La façon la plus simple d'introduire le couplage mécano-chimique dans la loi de comportement est de considérer l'hypothèse de la partition des déformations (utilisée en petites déformations) :
$$ \underline{\underline{\epsilon}} = \underline{\underline{\epsilon}}^{ep} + \underline{\underline{\epsilon}}^{th} + \underline{\underline{\epsilon}}^{ch} \quad (46) $$

$$ \text{ou } \underline{\underline{\epsilon}} = \frac{1 + v_o}{E_o} (\underline{\underline{\sigma}} - \sigma_o \underline{\underline{l}}) - \frac{v_o}{E_o} \text{tr}(\underline{\underline{\sigma}} - \sigma_o \underline{\underline{l}}) \underline{\underline{l}} + \frac{b}{K_o} \frac{P - P_o}{3} \underline{\underline{l}} + \alpha_o^t (T - T_o) \underline{\underline{l}} + \alpha_o^c \underline{\underline{\xi}} \underline{\underline{l}} \quad (47) $$

avec respectivement $\underline{\underline{\epsilon}}^{ep}$ la déformation dite poro-élastique donnée par :

$$ \underline{\underline{\epsilon}}^{ep} = \frac{1 + v_o}{E_o} (\underline{\underline{\sigma}} - \sigma_o \underline{\underline{l}}) - \frac{v_o}{E_o} \text{tr}(\underline{\underline{\sigma}} - \sigma_o \underline{\underline{l}}) \underline{\underline{l}} + \frac{b}{K_o} \frac{P - P_o}{3} \underline{\underline{l}} \quad (48) $$

$E_o$ : module d'Young drainé (P=P₀); $v_o$ : coefficient de Poisson drainé; P : pression de pores, P₀ : pression de pores initiale.

$\underline{\underline{\epsilon}}^{th}$, la déformation d'origine thermique :

$$ \underline{\underline{\epsilon}}^{th} = \alpha_o^t (T - T_o) \underline{\underline{l}} \quad (49) $$

$\alpha_o^t$: coefficient de dilatation thermique drainée.

$\underline{\underline{\epsilon}}^{ch}$, la déformation d'origine chimique :

$$ \underline{\underline{\epsilon}}^{ch} = \alpha_o^c \underline{\underline{\xi}} \underline{\underline{l}} \quad (50) $$

$\alpha_o^c$ : coefficient de dilatation chimique drainée (obtenue si P=P₀).

On vérifie que $\underline{\underline{\epsilon}}^{ch}$ est la déformation de l'échantillon, libre de contrainte dans des conditions isothermes et drainées (P=P₀).

En inversant la relation (47), on obtient la loi de comportement mécanique exprimée de la façon suivante :

$$ \underline{\underline{\sigma}} = \sigma_o \underline{\underline{l}} + (K_o - \frac{2}{3}\mu)\underline{\underline{\epsilon}} + 2\mu\underline{\underline{\epsilon}} - b(P - P_o)\underline{\underline{l}} - \alpha_o^c K_o (T - T_o)\underline{\underline{l}} - K_o \alpha_o^c \underline{\underline{\xi}} \underline{\underline{l}} \quad (51) $$

Le terme $\alpha_o^c K_o \underline{\underline{\xi}}$ représente la contrainte à exercer au système lorsqu'une réaction chimique se produit dans une expérience isotherme, en maintenant constante la pression de pores (P=P₀) et la déformation ($\underline{\underline{\epsilon}}$=0).

Par ailleurs, la loi de comportement mécanique s'écrit à l'aide de la théorie de la poro-élasticité de la manière suivante (Coussy, 1991) :

$$ \underline{\underline{\sigma}} = \sigma_o \underline{\underline{l}} + \left(K - \frac{2}{3}\mu\right)\underline{\underline{\epsilon}} + 2\mu\underline{\underline{\epsilon}} - bM \frac{m}{\rho_o} \underline{\underline{l}} - \alpha^t K (T - T_o) \underline{\underline{l}} \quad (52) $$

m : apport de masse **totale** de samure.

### 3.3 Loi de comportement hydraulique

Une façon simple d'établir la loi de comportement dite hydraulique (associée à la pression de pores) est d'imaginer une transformation isotherme, sans apport extérieur de masse de saumure, et en bloquant la déformation ($\underline{\underline{\epsilon}}$=0). L'identification des deux formes de la loi de comportement mécanique (47) et (52) permet d'obtenir :
$$
\left.
\begin{aligned}
\underline{\underline{\sigma}} &= \sigma_o \underline{\underline{l}} - b \frac{M}{\rho_o^{fl}} \xi \underline{\underline{l}} \\
\underline{\underline{\sigma}} &= \sigma_o \underline{\underline{l}} - b(P-P_o) \underline{\underline{l}} - K_o \alpha_o^c \xi \underline{\underline{l}}
\end{aligned}
\right\} \Rightarrow P-P_o = -\frac{\alpha_o^c K_o}{b} \xi + \frac{M}{\rho_o^{fl}} \xi
$$

Cette relation exprime la variation de pression de pores consécutive à la seule réaction chimique, pour la transformation envisagée. Dans le cas général, cette contribution s'ajoute linéairement aux effets du chargement thermo-hydro-mécanique. On a alors :

$$
P-P_o = M(-b \operatorname{tr}\underline{\underline{\epsilon}} + \frac{m^{\text{ext}}}{\rho_o^{\text{fl}}}) + \frac{\alpha_o^t K - \alpha_o^t K_o}{b} (T-T_o) - \frac{\alpha_o^c K_o}{b} \xi + M \frac{\xi}{\rho_o^{\text{fl}}} \quad (53)
$$

$\rho_o^{\text{fl}}$: masse volumique de la saumure; $m^{\text{ext}}$: apport de masse hors réaction chimique.

Ou encore en introduisant la variable m, apport total de masse de saumure :

$$
P - P_o = M(-b \operatorname{tr}\underline{\underline{\epsilon}} + \frac{m}{\rho_o^{\text{fl}}}) + \frac{\alpha_o^t K - \alpha_o^t K_o}{b} (T - T_o) - \frac{\alpha_o^c K_o}{b} \xi \quad (54)
$$

Afin d'établir ultérieurement l'expression de la densité d'énergie libre Φ,
développons la densité d'enthalpie massique gm au premier ordre :

$$
g_m = g_m^o + \frac{1}{\rho_o^{fl}}(P - P_o) + (-s_m^o)(T - T_o) + \mu_s^o(\chi - \chi_o) \quad (55)
$$

-$s_m^o$: entropie massique de la saumure dans l'état de référence. $\mu_s^o$: potentiel chimique massique du soluté dans l'état de référence. $\chi$: concentration massique de soluté.

Le terme $\mu_s^o(\chi - \chi_o)$ constitue une contribution liée à l'apport (ou cession) de masse de sel par la réaction chimique ou par transport diffusif depuis les systèmes élémentaires contigus. En prenant une forme linéaire :

$$
\mu_s^o(\chi - \chi_o) = \mu_s^o F_m \frac{m}{\phi_o \rho_o^{fl}} \quad (56)
$$

$F_m$ : coefficient de proportionnalité ; $\phi_o$ : porosité initiale.

En remplaçant les relations (54) et (56) dans l'expression (55) de gm, on obtient :

$$
\begin{equation}
\begin{split}
g_m = g_m^o & - b \frac{M}{\rho_o^{fl}} \operatorname{tr}\underline{\underline{\epsilon}} + \frac{M}{\rho_o^{fl}} \left( \frac{1}{\rho_o^{fl}} + \frac{F_m}{\phi_o \rho_o^{fl}} \right) m + \\
& \qquad + \left[ (-s_m^o) - \frac{\alpha_o^t K_o}{b} \right] (T - T_o) - \frac{\alpha_o^c K_o}{b} \xi
\end{split}
\tag{57}
\end{equation}
$$

3.4 Loi de comportement thermique

Il s'agit d'établir une relation linéaire de l'entropie de système S en fonction des variables d'état. Dans un milieu poreux, considérons une nouvelle fois une transformation isotherme (T=T₀) et à déformation bloquée (ϵ=0) du système élémentaire fermé volontairement. L'entropie qu'apporte la réaction chimique s'écrit alors :
$$dS = \frac{\Re}{T_o} \xi \qquad (58)$$

Re : chaleur de réaction non drainé.

Cette variation d'entropie s'ajoute aux contributions des autres variables d'état :

$$S = S_o + C_{\varepsilon} \frac{m_o}{T_o} (T - T_o) + \alpha^t K \operatorname{tr}\underline{\underline{\varepsilon}} + m_s^o - \frac{Lm}{T_o} + \frac{\Re}{T_o} \xi \qquad (59)$$

C<sub>ε</sub> : chaleur spécifique du système à déformation constante; m<sub>o</sub> : apport initial de saumure hors réaction chimique; L : chaleur latente par unité de masse de saumure apportée pour maintenir la température constante dans une expérience isotherme et à déformation constante.

## 3.5 Détermination du potentiel thermodynamique Φ

Dans le cadre d'une théorie linéaire, pour que les lois de comportement établies obéissent aux principes de la thermodynamique, il convient de vérifier qu'elles dérivent d'une fome quadratique de la densité d'énergie Φ. Prenons donc la forme suivante :

$$\begin{aligned} \Phi = \Phi_o + \underline{\underline{\varepsilon}}: \underline{\underline{A}}_{II} : \underline{\underline{\varepsilon}} + A_{22}(T-T_o)^2 + A_{33}m^2 + A_{44}\xi^2 + \\ + A_{11}\operatorname{tr}\underline{\underline{\varepsilon}} + A_{22}(T-T_o) + A_{33}m + A_{44}\xi + A_{12}\operatorname{tr}\underline{\underline{\varepsilon}}(T-T_o) + A_{13}\operatorname{tr}\underline{\underline{\varepsilon}}m + \\ + A_{14}\operatorname{tr}\underline{\underline{\varepsilon}}\xi + A_{23}(T-T_o)m + A_{24}m\xi + A_{25}(T-T_o)\xi \end{aligned} \qquad (60)$$

Avec la loi de comportement mécanique (52), on utilise la relation suivante :

$$\underline{\underline{\sigma}} = \frac{\partial \Phi}{\partial \underline{\underline{\varepsilon}}}$$

à l'aide de laquelle on identifie :

$$\underline{\underline{A}}_{II} = \frac{1}{2} \underline{\underline{C}}; \quad A_1 = \sigma_o; \quad A_{12} = -\alpha^t K; \quad A_{13} = -b \frac{M}{\rho_o^{fl}}; \quad A_{14} = 0 \qquad (61)$$

C : tenseur d'élasticité non drainé.

Avec la loi de comportement (57), la relation suivante : $g_m = \frac{\partial \Phi}{\partial m}$ est employée pour identifier les coefficients :

$$\begin{aligned} A_{33} &= \frac{1}{2} \frac{M}{\rho_o^{fl}} \left( \frac{1}{\rho_o^{fl}} + \frac{F_m}{\phi_o \rho_o^{fl}} \right); \quad A_3 = g_m^o; \quad A_{13} = -b \frac{M}{\rho_o^{fl}}; \quad A_{23} = \left[ (-s_m^o) + \frac{\alpha^t K - \alpha_o^t K_o}{b} \right]; \\ A_{24} &= - \frac{\alpha_o^c K_o}{b} \end{aligned} \qquad (62)$$

La loi de comportement thermique (59) et la relation $S = -\frac{\partial \Phi}{\partial T}$ conduisent à :

$$A_{24} = \frac{1}{2} C_{\varepsilon} \frac{m_o}{T_o}; \quad A_2 = -S_o; \quad A_{12} = -\alpha^t K; \quad A_{23} = -s_m^o + \frac{L}{T_o}; \quad A_{25} = -\frac{\Re}{T_o} \qquad (63)$$

L'identification du coefficient A<sub>23</sub> dans (62) et (63) permet d'établir une relation entre la chaleur latente L et les différents coefficients de dilatation thermique :
$$L = \frac{T_o}{\rho_o^{fl}} \left[ \frac{\alpha_o^t K^t - \alpha_o^t K_o^t}{b} \right] \qquad (64)$$

Si la loi de comportement chimique ($A = \frac{\partial \Phi}{\partial \xi}$) dérive bien de la forme quadratique (60), on a alors :

$$A_4 = A_o ; \quad A_{44} = \frac{1}{2} \left. \frac{\partial A}{\partial \xi} \right|_{T, \underline{\underline{\epsilon}}, m} \qquad (65)$$

avec le développement de l'affinité A suivant :

$$A = A_o + \left. \frac{\partial A}{\partial \xi} \right|_{T, \underline{\underline{\epsilon}}, m} \xi - \frac{\Re}{T_o} (T - T_o) - \frac{\alpha_o^c K_o}{b} m \qquad (66)$$

Cette relation appelle deux commentaires :

(1) la loi (66) suppose implicitement que la vitesse de réaction durant la transformation qui amène le système d'un état indicé 0 à l'état actuel, est constante (en effet, le temps n'intervient pas);

(2) on constate que la cinétique chimique est couplée à l'apport de masse ($-\frac{\alpha_o^c K_o}{b} m$).

## 4. ÉCRITURE DES ÉQUATIONS D'ÉVOLUTION DES GRANDEURS D'ÉTAT

Après avoir obtenu les lois de comportement thermodynamiquement admissibles dans le cadre d'une théorie linéaire, il s'agit d'utiliser les principes de conservation et les lois de transport afin d'établir les équations locales d'évolutions des variables d'état (ou mieux de grandeurs d'état mesurables). Dans un calcul de structures, ce sont ces équations d'évolution qui doivent être résolues pour quantifier les grandeurs d'état en tout point de la structure et à chaque instant.

### 4.1 Équation d'évolution hydraulique

Afin d'établir l'équation d'évolution hydraulique, on isole l'apport de masse de saumure m dans la loi de comportement hydraulique (54) :

$$\frac{m}{\rho_o^{fl}} = \frac{1}{M}(P - P_o) + b \operatorname{tr} \underline{\underline{\epsilon}} - \frac{\alpha_o^t K - \alpha_o^t K_o}{Mb}(T - T_o) + \frac{\alpha_o^c K_o}{Mb} \xi \qquad (67)$$

La conservation de la masse de saumure s'écrit :

$$\dot{m} = - \operatorname{Div} \underline{\underline{M}} + \dot{\xi} \qquad (68)$$

$\dot{\xi}$ : l'apport de masse due à la réaction par unité de temps.

M : le flux de masse de saumure, ici supposé obéir à la loi de Darcy :
$$ \underline{M} = -\rho_o^{fl} \frac{k}{\eta} \underline{grad} P \qquad (69) $$

ou encore à la suite du § 2.5.b; dans un cadre plus général et en prenant en compte simultanément tous les gradients :

$$ \underline{M} = -\rho_o^{fl} \frac{k}{\eta} \underline{grad} P - K_c^p \underline{grad} C - K_t^p \underline{grad} T $$

En remplacant les relations (67) et (69) dans (68), on obtient l'équation d'évolution hydraulique, en tenant compte de (64) :

$$ \frac{1}{M} \dot{P} + b \dot{t} \underline{\underline{\epsilon}} - \frac{(\rho_o^{fl})^2 L}{M T_o} \dot{T} + (-1 + \frac{\alpha_o^c K_o}{M b}) \dot{\xi} = \text{Div}(\frac{k}{\eta} \underline{\underline{grad}} P) \qquad (70a) $$

ou encore en tenant compte des termes de couplage dans la loi de Darcy :

$$ \frac{1}{M} \dot{P} + b \dot{t} \underline{\underline{\epsilon}} + \left(1 - \frac{(\rho_o^{fl})^2 L}{M T_o}\right) \dot{T} + \left(-1 + \frac{\alpha_o^c K_o}{M b}\right) \dot{\xi} = \text{Div}\left(\frac{k}{\eta} \underline{\underline{grad}} P + \frac{K_c^p}{\rho_o^{fl}} \underline{\underline{grad}} C + \frac{K_t^p}{\rho_o^{fl}} \underline{\underline{grad}} T\right) \qquad (70b) $$

## 4.2 Équation d'évolution de la masse de soluté

Soit dans le système élémentaire Ω, une solution saline occupant le volume Vp à la concentration volumique C. La masse de sel en solution dans le volume du système vérifie :

$$ M_I = C V_p \qquad (71) $$

En différentiant, en linéarisant, puis en divisant par le volume initiale du système, Vo, on obtient alors une expression de mI, apport de masse de sel en solution par unité de volume initial au cours d'une transformation quelconque :

$$ m_I = C_o \phi_o \frac{dV_p}{V_p} + \phi_o dC \qquad (72) $$

$\phi_o$ : porosité connectée initiale; $C_o$ : concentration volumique initiale.

La variation de volume de pores ou du fluide (puisque le milieu est saturé) a déjà été donnée pour une transformation isotherme (Marsily, 1986; ou Charlez, 1991). Si on ajoute un terme d'origine chimique supplémentaire, on a alors :

$$ \frac{dV^p}{V^p} = \frac{b}{\phi_o K_o} d\sigma + \left( \frac{b}{\phi_o K_o} - \frac{1}{K_s} \right) dP + \alpha_t^{fl} dT + \left( \frac{dV^p}{V^p} \right)^{ch} \qquad (73) $$

$d\sigma$ : variation de la contrainte moyenne par rapport à un état de référence donné. $dP$ : variation de la pression de pores; $K_o$ : module de la compressibilité drainée; $\phi_o$ : porosité initiale; $dT$ : variation de température; $\alpha_t^{fl}$ : coefficient de dilatation volumique de la saumure.

Afin de faire apparaître la déformation du système (matrice + pores vides de fluide), on écrit la loi de comportement mécanique (48), après en avoir pris la trace :
$$d\sigma = K_o\varepsilon - b \, dP - \alpha_o^t K_o dT - \alpha_o^c K_o \xi \quad (74)$$

En remplaçant (74) dans (73) :

$$\frac{dV^p}{V^p} = \frac{b}{\phi_o} \varepsilon + \left( \frac{b-\phi_o}{\phi_o K_o} \right) dP + \left( \alpha_{fl}^t - \frac{\alpha_o^t}{\phi_o} \right) dT - \frac{b}{\phi_o} \varepsilon^{ch} + \left( \frac{dV^p}{V^p} \right)^{ch} \quad (75)$$

En considérant de nouveau les relations :

$$\varepsilon^{ch} = \alpha_o^c \xi \quad \text{et} \quad \left( \frac{dV^p}{V^p} \right)^{ch} = \alpha_{fl}^c dC$$

avec la relation (75), l'équation (72) peut se réécrire :

$$m_I = C_o b \varepsilon + C_o \left( \frac{b-\phi_o}{K_o} \right) dP + C_o (\phi_o \alpha_{fl}^t - \alpha_o^t) dT - C_o b \alpha_o^c \xi + \phi_o (1 + C_o \alpha_{fl}^c) dC \quad (76)$$

La conservation de la masse de sel en solution s'écrit :

$$\dot{m}_I = - \mathrm{Div} \underline{M}_I + \dot{\xi} \quad (77)$$

$\underline{M}_I$ : le flux de masse de sel en solution donné par (Marsily, 1986) :

$$\underline{M}_I = - \underline{\underline{D}} \underline{\underline{grad}} C + C \underline{\underline{U}} \quad (78)$$

D : le coefficient de dispersion cinématique; $\underline{\underline{U}}$ : la vitesse de Darcy ($\underline{\underline{M}}/p_{fl}$ voir 69).

ou encore à la suite du § 2.5.b; dans un cadre plus général et en prenant en compte simultanément tous les gradients :

$$\underline{\underline{M}}_I = - \underline{\underline{D}} \underline{\underline{grad}} C + C \underline{\underline{U}} - K_t^c \underline{\underline{grad}} T$$

En remplaçant (78) et (76) dans (77), on obtient la loi d'évolution recherchée :

$$\phi_o(1+C_o\alpha_{fl}^c)\dot{C} + C_o b \dot{\varepsilon} + C_o\left(\frac{b-\phi_o}{K_o}\right)\dot{P} + C_o(\phi_o\alpha_{fl}^t - \alpha_o^t)\dot{T} \\ = \mathrm{Div}(\underline{\underline{D}}\underline{\underline{grad}}C - \underbrace{C\underline{\underline{U}}}_{\text{terme convectif}}) + (1+C_o b \alpha_o^c)\dot{\xi} \quad (79a)$$

ou encore avec les termes de couplages dans la loi de transport :

$$\phi_o(1+C_o\alpha_{fl}^c)\dot{C} + C_o b \dot{\varepsilon} + C_o\left(\frac{b-\phi_o}{K_o}\right)\dot{P} + C_o(\phi_o\alpha_{fl}^t - \alpha_o^t)\dot{T} \\ = \mathrm{Div}(\underline{\underline{D}}\underline{\underline{grad}}C - \underbrace{C\underline{\underline{U}}}_{\text{terme convectif}} + K_t^c \underline{\underline{grad}}T) + (1+C_o b \alpha_o^c)\dot{\xi} \quad (79b)$$

## 4.3 Équation d'évolution de la température

La première étape consiste à écrire l'équation de conservation de l'énergie et la dissipation intrinsèque pour les variables d'état choisies :
**conservation d'énergie :**

$$ \dot{e} = \underline{\underline{\sigma}} : \dot{\underline{\underline{\epsilon}}}} - \text{Div}(h_m \underline{\underline{M}}) + R - \text{Div}q $$

e : densité d'énergie; hm : enthalpie massique; R : source interne d'énergie négligée dans la suite; q : flux de chaleur apporté par conduction. Les forces de volume ne sont pas pris en compte.

**dissipation intrinsèque :**

$$ \Psi_I = \underline{\underline{\sigma}} : \dot{\underline{\underline{\epsilon}}}} - S\dot{T} - \dot{\Phi} + g_m \dot{m} + A\chi $$

En isolant la dissipation mécanique $\underline{\underline{\sigma}} : \dot{\underline{\underline{\epsilon}}}}$ dans l'expression de la dissipation intrinsèque et en la substituant dans l'expression de la conservation de l'énergie, une forme générale de la conservation de la chaleur peut être obtenue :

$$ \underbrace{T_o \text{Div}[\underline{\underline{M}} s_m]}_{\text{Apport convectif d'entropie}} + \underbrace{\dot{S} T_o}_{\substack{\text{Variation d'entropie} \\ \text{du matériau poreux}}} + \text{Div}q = \psi \quad (80) $$

Le terme de droite de l'égalité représente l'énergie dissipée par chaleur à cause d'irréversibilités d'origine plastique ou chimique ou encore associée aux forces visqueuses du mouvement du fluide. En pratique, elle est souvent négligée devant l'apport des termes à gauche du signe égal.

En remplaçant dans (80), la loi de comportement thermique (59), en utilisant la conservation de la masse de saumure (68) et la loi de Fourier, l'équation d'évolution de la température s'écrit :

$$ (m_o C_\varepsilon + \frac{\rho_o^{fl} L^2}{\eta T_o}) \dot{T} + (\alpha K T - L \rho_o^{fl} b) \dot{\varepsilon} - \frac{\rho_o^{fl} L}{M} \dot{P} + \frac{\Re}{T_o} \dot{\xi} = \text{Div}(\lambda \underbrace{\text{grad } T}_{\text{terme convectif}} - \rho_o^{fl} c^{\underline{\underline{f}} U}) \quad (81a) $$

$c^{\underline{\underline{f}}}$ : capacité calorifique du fluide.

ou encore en introduisant les termes de couplage dans la loi de Fourier :

$$ (m_o C_\varepsilon + \frac{\rho_o^{fl} L^2}{\eta T_o}) \dot{T} + (\alpha K T - L \rho_o^{fl} b) \dot{\varepsilon} - \frac{\rho_o^{fl} L}{M} \dot{P} + \frac{\Re}{T_o} \dot{\xi} = \text{Div}(\lambda \underbrace{\text{grad } T}_{\text{terme convectif}} + K_p^{\underline{\underline{t}} \text{grad } T} + K_c^{\underline{\underline{t}} \text{grad } P} - \rho_o^{fl} c^{\underline{\underline{f}} U}) \quad (81b) $$

## 4.4 Équation de Navier

L'équilibre mécanique se traduit par l'équation :

$$ \text{Div}(\underline{\underline{\sigma}} - \underline{\underline{\sigma}}_o) + m\underline{\underline{E}} = 0 \quad (82) $$

E : efforts massique (dûs en général à la pesanteur).

L'inconnue mécanique que nous avons choisie est le champs de déplacement $\underline{\underline{u}}$. En utilisant la décomposition du tenseur des déformations en fonction de $\underline{\underline{u}}$ et la loi de comportement mécanique (5) rappelée ci-dessous :

$$ \underline{\underline{\sigma}} = \sigma_o \underline{\underline{l}} + (K_o - \frac{2}{3}\mu)\underline{\underline{\varepsilon}} + 2\mu\underline{\underline{\varepsilon}} - b(P - P_o)\underline{\underline{l}} - \alpha_o^c K_o (T - T_o)\underline{\underline{l}} - K_o \alpha_o^c \underline{\underline{\xi}}\underline{\underline{l}} $$

On obtient l'équation dite de Navier :
$$ 
\begin{aligned}
& \left( K_o + \frac{4}{3}\mu \right) \underline{\underline{\text{grad}(div u)}} - \mu \underline{\underline{\text{rot}(rot u)}} - b \underline{\underline{\text{grad}}}(P - P_o) \\
& \quad - 3\alpha_o^t K_o \underline{\underline{\text{grad}}}(T - T_o) - 3\alpha_o^c K_o \underline{\underline{\text{grad}}}\xi + m\underline{F} = 0
\end{aligned} 
\quad (83) $$

## 4.5 Équation de cinétique chimique

Pour compléter l'ensemble des équations d'évolution, il manque finalement une équation de cinétique chimique relative à la variable apport de masse de sel de la solution par unité de volume initial de milieu poreux ξ. Celle établie en (66) est peu exploitable : elle introduit des paramètres difficilement quantifiables. On a donc recours à une loi plus simple qui a le mérite d'être utilisée en pratique. La vitesse de dissolution (ou de cristallisation) d'une quantité de sel m<sub>s</sub> (James, 1992) s'écrit :

$$ \frac{dm_s}{dt} = K_c S(C_s - C) $$

K<sub>c</sub> : vitesse de dissolution (m.s<sup>-1</sup>); S : surface du cristal en contact avec la solution saline (m<sup>2</sup>); C<sub>s</sub>-C : sous-concentration volumique en supposant une cinétique chimique du premier ordre, ce qui est le cas pour la réaction de dissolution du chlorure de sodium (Davion, 1953), C<sub>s</sub> : solubilité de la solution saline (kg.m<sup>-3</sup>).

Cette équation se généralise en terme de concentration volumique de la façon suivante :

$$ \frac{dC}{dt} = K_c \frac{S}{V} (C_s - C) \quad (84) $$

S/V : rapport de la surface du cristal mis en contact avec la solution rapportée au volume de solution (m<sup>-1</sup>).

En introduisant la variable ξ, l'équation (84) se réécrit dans le cas d'une dissolution :

$$ \dot{\xi} = \phi_o \frac{S}{V} K_c (C_s - C) \quad (85) $$

Ce terme est positif dans le cas d'une sous-saturation (C<sub>s</sub>>C) de la solution provoquant alors une dissolution et un apport de masse en solution. Dans le cas d'une sursaturation (C<sub>s</sub><C), le coefficient de cinétique K<sub>c</sub> se trouve changé ; la solution cède alors de la masse à la phase solide.

## 4.6 Récapitulation

En introduisant l'équation de cinétique chimique (85) dans les équations d'évolution (70), (79), (81), et (83), on obtient le système de 4 équations (3 scalaires et 1 vectorielle) à 4 inconnues P, T, C et u :

Équation d'évolution hydraulique :

$$ 
\begin{aligned}
& \frac{1}{M} \dot{P} + b \dot{\varepsilon} - \frac{(\rho_o^{fl})^2 L}{M T_o} \dot{T} + (-1 + \frac{\alpha_o^c K_o}{M b}) \phi_o \frac{S}{V} K_c (C_s - C) \\
& \quad = \text{Div} \left( \frac{k}{\eta} \underline{\underline{\text{grad}P}} + \frac{K_c^p}{\rho_o^{fl}} \underline{\underline{\text{grad}C}} + \frac{K_t^p}{\rho_o^{fl}} \underline{\underline{\text{grad}T}} \right)
\end{aligned} 
\quad (86) $$
## Équation d'évolution de la masse de soluté

$$ 
\begin{aligned}
& \phi_o(1+C_o\alpha_{fl}^c)\dot{C} + C_o b \dot{\epsilon} + C_o\left(\frac{b-\phi_o}{K_o}\right)\dot{P} + C_o(\phi_o\alpha_{fl}^t - \alpha_o^t)\dot{T} \\
& = \operatorname{Div}(D\underline{\underline{\mathrm{grad}}C} - \underbrace{CU}_{\text{terme convectif}} + K_t^c\underline{\underline{\mathrm{grad}}T}) + (1+C_o b\alpha_o^c)\phi_o \frac{S}{V} K_c(C_s-C)
\end{aligned}
\quad (87)
$$

## Équation d'évolution de la température

$$ 
\begin{aligned}
& (m_o C_\varepsilon + \frac{\rho_o^{fl} L^2}{\eta T_o})\dot{T} + (\alpha K T - L \rho_o^{fl} b) \dot{\epsilon} - \frac{\rho_o^{fl} L}{M} \dot{P} + \frac{\Re}{T_o} \phi_o \frac{S}{V} K_c(C_s - C) \\
& = \operatorname{Div}(\lambda \underline{\underline{\mathrm{grad}}T} + K_c^t \underline{\underline{\mathrm{grad}}C} + K_c^t \underline{\underline{\mathrm{grad}}T} - \rho_o^{fl} c^{fl} \frac{k}{\eta} \underline{\underline{\mathrm{grad}}P})
\end{aligned}
\quad (88)
$$

## Équations de Navier

$$ \left( K_o + \frac{4}{3}\mu \right) \underline{\underline{\mathrm{grad}}}(\underline{\underline{\mathrm{div}}} \dot{\underline{\underline{u}}}) - \mu \underline{\underline{\mathrm{rot}}} (\underline{\underline{\mathrm{rot}}}}\dot{\underline{\underline{u}}}) - b \underline{\underline{\mathrm{grad}}} \dot{P} - 3 \alpha_o^t K_o \underline{\underline{\mathrm{grad}}} \dot{T} - 3 \alpha_o^c K_o \underline{\underline{\mathrm{grad}}} (\phi_o \frac{S}{V} K_c(C_s - C)) + \dot{m} \underline{\underline{F}} = 0 \quad (89) $$

NB : Losque les efforts massiques $\underline{\underline{F}}$ ne sont pas négligeables, on pourra remplacer la dérivée $\dot{m}$ par son expression (67).

## Nomenclature :

$\alpha_o^t$ : coefficient de dilatation thermique drainée;

$\alpha_o^c$ : coefficient de dilatation chimique drainée;

$\alpha_{fl}^c$ : coefficient de dilatation chimique de la phase liquide;

$\alpha^t$ : coefficient de dilatation thermique non drainé;

$\alpha^c$ : coefficient de dilatation thermique non drainé;

$\alpha_{fl}^t$ : coefficient de dilatation thermique de la phase fluide;

$b$ : coefficient de Biot;

$C_s-C$ : sous-concentration volumique;

$c^{fl}$ : capacité calorifique du fluide.

$\underline{\underline{\varepsilon}}$ : tenseur des déformations élastiques;

$\underline{\underline{\epsilon}}$ : trace du tenseur $\underline{\underline{\varepsilon}}$;

$\underline{\underline{\phi}}$ : porosité;

$\underline{\underline{\eta}}$ : viscosité dynamique;

$\underline{\underline{k}}$ : perméabilité intrinsèque;

$K_o$ : module de compressibilité drainé;

$K_c$ : vitesse de dissolution;

$K_b^a$ : coefficient de couplage du terme en gradb attachée à la loi de conduction sur la variable d'état a;

$\lambda$ : conductivité thermique;

$L$ : chaleur latente par unité de masse de saumure apportée pour maintenir la température constante dans une expérience isotherme et à déformation constante.

$\mu$ : module de cisaillement,

$m$ : apport de masse de saumure;

$M$ : module de Biot;

$P$ : pression du fluide;

$\rho^{fl}$ : masse volumique du fluide;

$\Re$ : chaleur de réaction non drainé.

$S/V$ : rapport de la surface de la solution mis en contact avec la solution rapportée au volume de solution;

$\sigma_o$ : contrainte initiale;

$T$ : température;

L'indice 0 désigne l'état initial dans lequel se trouvent les variables d'état.
## 4.7 Prise en compte des non linéarités fondamentales

A la lumière de l'étude bibliographique du chapitre I, les non linéarités qui doivent être prises en compte portent sur les propriétés visqueuses du sel et sur la perméabilité.

### a-Cas du comportement mécanique non linéaire

La déformation inélastique $\underline{\underline{\epsilon}}^{inel}$ est introduite par l'hypothèse de la partition des déformations :

$$ \underline{\underline{\epsilon}}^t = \underline{\underline{\epsilon}} + \underline{\underline{\epsilon}}^{inel} $$

où $\underline{\underline{\epsilon}}^t$ est la déformation totale. Il suffit alors de remplacer dans les lois de comportement le tenseur $\underline{\underline{\epsilon}}$ par $\underline{\underline{\epsilon}}^t - \underline{\underline{\epsilon}}^{inel}$. Le tenseur inélastique est obtenu lui par un modèle suivant la nature de la déformation inélastique (déformation plastique, endommagement, déformation par dissolution-cristallisation etc.).

Cette approche est nécessaire lorsqu'une déformation macroscopique est contrôlée par des processus en deçà de la taille du VER. Cette situation est rencontrée dans le cas de la déformation par dissolution-cristallisation contrôlée par des transferts par diffusion en solution à l'échelle du grain (voir § 2.1 du chapitre II). A titre d'illustration, donnons une relation de la forme suivante proposée par Gratier (1983):

$$ \dot{\epsilon}_{vp} = \alpha D_m \bar{c} V_s \frac{\Delta \sigma_n}{RT d^3} $$

$\epsilon_{vp}$ : vitesse de déformation viscoplastique ; $D_m$ : coefficient de diffusion des substances dissoutes en solution ; $\alpha$ : coefficient qui dépend de la forme des grains et de la répartition des contraintes (pour la calcite, l'auteur choisit $\alpha =32$) ; $\bar{c}=cV_1$ : concentration du solide en solution exprimée en rapport volumique ; $V_s$ : volume molaire du sel ; $V_1$ : volume molaire en solution du soluté ; $\Delta\sigma_n$ : différence de contrainte normale entre les zones de dissolution et les zones de cristallisation (en pratique, le déviateur de contrainte) ; $d$ : taille moyenne des grains;

### b- Prise en compte des non linéarités de la perméabilité :

On a pu voir dans l'étude bibliographique que la perméabilité était susceptible d'évoluer dans les deux cas suivants :

a- lorsque la roche s'endommage consécutivement à un chargement mécanique déviatorique. Peach (1991) a pu trouver expérimentalement au cours d'essais triaxiaux une relation de la forme suivante :

$$ k = \left( \frac{\Delta V}{V} \right)^n $$

avec $n > 4$ lorsque les déformations de l'échantillon sont inférieures à 0,2 %, $n < 1$ sinon. Les résultats expérimentaux de Peach auraient pu être interprétés différemment en attribuant l'évolution de la perméabilité à d'autres paramètres mécaniques, la distorsion irréversible par exemple ou le travail des déformations irréversibles déviatoriques ou volumiques. La relation ci-dessus est donc indicative et mérite d'être précisée par des études complémentaires.

b- lorsque la roche est saturée de saumure. Un changement de la solubilité par variation de pression et/ou de température et le processus de "cicatrisation" sous contrainte peuvent engendrer des phénomènes de dissolution-cristallisation colmatant le réseau connecté. Ces effets seront examinés dans le chapitre V.
## 5. ESTIMATION DES TEMPS CARACTÉRISTIQUES

Les équations d'évolution des grandeurs d'état ne peuvent être résolues d'une manière simple et on peut se demander dans quelles conditions, il devient possible de découpler un ou plusieurs processus de transfert (chaleur, masse de soluté ou masse fluide). Pour cela, une approche simple consiste à estimer et ensuite comparer les temps caractéristiques attachés à chacun de ces processus. Or, ces temps dépendent tout autant de la loi de comportement que de la structure et ne peuvent être évalués précisément dans le cas général. Un calcul d'ordre de grandeur est cependant possible. Cette démarche est mise en oeuvre dans le cas du sel gemme en écrivant les équations unidimensionnelles de transfert de masse et de chaleur sur un élément de volume du matériau poreux, représentant un VER. Par conséquent, cette discussion ne s'intéressera pas aux phénomènes locaux en deçà du VER. Par ailleurs, les termes apportés d'une part par la déformation chimique puisque non quantifiables et d'autre part par les termes croisés de couplage dans les lois de conduction, ne sont pas pris en compte.

a- l'équation unidimensionnelle de transfert de masse de soluté est donnée pour un sel ayant la porosité constante $\phi$ :

$$D \frac{\partial^2 C}{\partial x^2} - U \frac{\partial C}{\partial x} = \phi \frac{\partial C}{\partial t} + \phi K_c \frac{A}{V} C$$

C : concentration en sel de la solution (g.l⁻¹) D: coefficient de diffusion de sel en solution dans le milieu poreux (m².s⁻¹), Dans le cas de milieu très faiblement perméable, ce coefficient prend une valeur très proche de 10⁻¹⁰ m².s⁻¹ ; U : vitesse de Darcy (m.s⁻¹), K_c : constante de cinétique chimique proche de 10⁻⁵ m.s⁻¹ dans le cas d'une dissolution; les données d'Harano et al. (1993) suggèrent que les constantes de cristallisation pour un sel sans impuretés sont aussi du même ordre de grandeur ; A/V : rapport de la surface de grains mis en contact avec la solution par unité de volume de solution (en m⁻¹). Ce paramètre peut être estimé par des modèles géométriques de milieux poreux simple. Pour une distribution isotrope de microfissures, A/V= 2/ w où w est l'ouverture moyenne des microfissures prise égale à 10⁻⁶ m.

L'écriture de l'équation précédente implique que le Volume Élémentaire Représentatif (VER) du sel sur lequel on écrit l'équation de conservation de la masse de sel constitue un système ouvert. Cette hypothèse est discutable pour un sel sain pour lequel l'existence d'une porosité connectée n'est pas prouvée.

On peut mettre en évidence trois temps caractéristiques attachés respectivement à la diffusion de masse, à la convection et à la réaction chimique dont on trouvera l'expression et une estimation grossière dans le tableau ci-dessous.

La valeur de la vitesse de Darcy de 10⁻¹³ m.s⁻¹ est calculée arbitrairement pour une perméabilité de 10⁻²⁰ m² et un gradient hydraulique de 10⁻² MPa.m⁻¹ correspondant au gradient hydrostatique. Ce n'est qu'une valeur indicative et dépend bien sûr du problème hydraulique posé.

<table><thead><tr><th>Diffusion</th><th>Convection</th><th>Réaction</th></tr></thead><tbody><tr><td>τ<sub>D</sub> = <sup>L<sup>2</sup> φ</sup>/<sub>D</sub></td><td>τ<sub>c</sub> = <sup>L φ</sup>/<sub>U</sub></td><td>τ<sub>R</sub> = <sup>L V</sup>/<sub>K<sub>c</sub> A</sub></td></tr><tr><td>3 ans<br>en prenant L = 1 mètre (étude du<br>champs proche d'un stockage) et<br>φ=1%</td><td>3000 ans<br>(en prenant U=10<sup>-13</sup> m.s<sup>-1</sup>,<br>L=1 m et φ=1%)</td><td>10 s</td></tr></tbody></table>

En comparant ces temps caractéristiques, on constate que la réaction chimique est extrêmement rapide vis-à-vis de des autres processus de transfert de masse de soluté.
b- En l'absence de couplages thermo-mécaniques, l'équation de transfert de chaleur s'écrit :

$$D_T \frac{\partial^2 T}{\partial x^2} - U \frac{(\rho c)^{fl}}{\rho c} \frac{\partial T}{\partial x} = \frac{\partial T}{\partial t} + \frac{(-\Delta H)}{\rho c} \zeta$$

$D_T$: coefficient de diffusivité thermique du sel égale à $3.10^{-6} \text{ m}^2\cdot\text{s}^{-1}$; $\rho$: masse volumique du sel ($\approx 2,1610^3 \text{ kg}\cdot\text{m}^{-3}$); $c^{fl}$: la masse volumique du fluide ($\approx 1,2 \times 10^3 \text{ kg}\cdot\text{m}^{-3}$); $c$: la chaleur massique du sel ($\approx 880 \text{ J}\cdot\text{kg}^{-1}\cdot\text{K}^{-1}$); $c^{fl}$: la chaleur massique du fluide ($\approx 3200 \text{ J}\cdot\text{kg}^{-1}\cdot\text{K}^{-1}$); $-\Delta H$: l'enthalpie de la réaction (en J.kg⁻¹). Celle-ci dépend de la teneur en impuretés (notamment de l'anhydrite) et de l'état de saturation de la solution. On prendra la valeur moyenne donnée par l'A.T.G. (1986) de $3.10^4 \text{ J}\cdot\text{kg}^{-1}$.

Le paramètre $\zeta$, taux de dissolution/cristallisation est donné par :

$$\zeta = \phi K_c \frac{A}{V} C$$

Où C est la concentration en sel de la solution et vérifie par conséquent l'équation de transfert de la masse de sel.

<table><thead><tr><th>Diffusion</th><th>Convection</th><th>Chaleur de dissolution</th></tr></thead><tbody><tr><td>$\tau_D = \frac{L^2 \rho c}{\lambda_T}$</td><td>$\tau_c = \frac{L}{U} \frac{\rho c}{(\rho c)^{fl}}$</td><td>$\tau_R = \frac{\rho c T_0}{\zeta(-\Delta H)}$</td></tr><tr><td>$\approx \text{jours}$<br>$(L=1\text{m})$</td><td>$3000 \text{ ans}$<br>$(U=10^{-13} \text{ m}\cdot\text{s}^{-1}; L=1\text{m})$</td><td>$\approx \text{heure}$<br>$(\delta C=1\text{g}\cdot\text{l}^{-1})$</td></tr></tbody></table>

La chaleur de dissolution semble avoir un rôle limité sur l'évolution à long terme de la température. Néanmoins, il convient de rester prudent compte tenu des incertitudes qui pèsent sur le paramètre (-ΔH).

c- En négligeant les termes liés aux couplages thermo-mécaniques, l'équation de diffusion hydraulique s'écrit :

$$\frac{k}{\eta} \frac{\partial^2 P}{\partial x^2} = \left( \frac{\phi}{K^{fl}} + \frac{b-\phi}{K^s} \right) \frac{\partial P}{\partial t}$$

où $k$: la perméabilité intrinsèque du sel (pour un sel faiblement endommagé $\approx 10^{-20} \text{ m}^2$); $\eta$: la viscosité dynamique de la saumure ($\approx 1,2 \times 10^{-3} \text{ Pa}\cdot\text{s}$); $b$: le coefficient de Biot ($\approx 0,1$ pour le sel); $K^s$: le module d'incompressibilité de la matrice solide ($\approx 25 \text{ GPa}$); $c^{fl}$: le module de compressibilité du fluide ($\approx 3,7 \text{ GPa}$).

A partir de l'équation (8) et des valeurs numériques estimées, on dégage le temps caractéristique de la diffusion hydraulique suivant :

$$\tau_D = \left( \frac{\phi}{K^{fl}} + \frac{b-\phi}{K^s} \right) \frac{\eta L^2}{k} \approx 12 \text{ jours}$$

Compte tenu de l'incertitude sur la perméabilité intrinsèque, la valeur du temps caractéristique est sensiblement du même ordre de grandeur que celle obtenue pour la diffusion thermique.

On retiendra de cette analyse les trois points suivants :
- les transports par convection interviennent probablement sur des échelles de temps géologiques,

- les temps caractéristiques des processus de diffusion s'inscrivent dans une échelle de temps humaine,

- la réaction chimique est quasi instantanée en regard des durées caractéristiques des deux processus précédents. Ce résultat nous permet de considérer la saumure piégée dans le réseau poreux comme étant en équilibre chimique à chaque instant, à l'échelle du VER.

## CONCLUSION

Dans un objectif de calcul de structures, nous avons établi un cadre théorique permettant de prendre en compte le caractère faiblement perméable et chimiquement réactif du sel gemme. Les dimensions caractéristiques qui nous intéressent (de l'ordre d'un mètre à plusieurs centaines de mètres) sont beaucoup plus grandes que les dimensions caractéristiques des microfissures connectées et des joints de grains (typiquement inférieurs à quelques centimètres). Aussi était-il justifié d'assimiler ici le sel, tout au moins endommagé, à un milieu poreux. Cette assertion est confortée par une estimation de la taille du VER.

Ce cadre s'appuyant sur la thermodynamique des milieux poreux a permis d'obtenir les équations d'évolution de quatre grandeurs d'état mesurables *in situ* : la température, la pression de pores, le déplacement et la concentration en sel de la solution. Dès lors que les nombreux paramètres sont identifiés expérimentalement, le système d'équation proposé peut être résolu par une méthode numérique adéquate. Les résultats d'une analyse simplifiée des temps caractéristiques suggèrent que dans une première approche, pour une étude à moyen terme, les effets de la cinétique chimique et les transferts par convection peuvent être négligés.

Toutefois, ce cadre théorique souffre, avouons le, d'un handicap majeur pour être utilisé facilement en vue d'une modélisation de processus couplés régissant la sûreté à long terme d'un stockage. De nombreux paramètres, difficiles à mesurer au laboratoire ou *in situ* y interviennent et on ne pourrait mettre en œuvre le cadre proposé sans recours à une campagne expérimentale longue, délicate et onéreuse.

Ainsi, aurait-on intérêt à simplifier, chaque fois que cela s'avérerait justifié, le modèle en hiérarchisant les processus et les phénomènes mis en jeu. On obtiendrait ainsi des variantes du modèle proposé relativement simplifiées et faciles à caractériser expérimentalement.

Une telle approche ne serait justifiée que dans un contexte précis. Elle ne pourrait être généralisée à tout phénomène et à toute échelle de temps ou d'espace. Ce propos sera illustré à travers quelques exemples concrets présentés dans le chapitre qui suit.
# CHAPITRE V

## INCIDENCE DES PHÉNOMÈNES COUPLÉS SUR LE COMPORTEMENT DES OUVRAGES SOUTERRAINS
# INTRODUCTION

Qu'il s'agisse d'une mine, d'une cavité de stockage d'hydrocarbures ou d'un stockage de déchets industriels ou nucléaires, la sûreté à long terme est intimement liée à la perméabilité des terrains et à son évolution sous l'effet de couplages entre les processus de transfert de masse et le comportement mécanique.

Nous examinons ici l'influence de la perméabilité (même très faible) et du caractère soluble du sel sur la sûreté de telles structures souterraines. Pour cela, nous nous appuierons autant que possible sur le modèle théorique décrit dans le chapitre précédent : le sel est considéré comme un milieu poreux saturé et réactif.

Dans une première partie, l'effet d'une interaction géochimique sur la perméabilité d'une roche endommagée est étudié. Dans le cadre d'une analyse de sûreté, il est nécessaire de savoir si cet effet favorise une augmentation ou au contraire une diminution de la perméabilité.

Dans une deuxième partie, la perméabilité mesurée expérimentalement au moins au voisinage d'un ouvrage souterrain constitue le point de départ de l'étude. A travers des calculs sur des cavernes de géométrie simple, l'évolution à court et à moyen terme de la pression de fluide (saumure ou hydrocarbures) dans une caverne abandonnée est discutée.

# PARTIE I. ÉVOLUTION DE LA PERMÉABILITÉ D'UN SEL MICROFISSURÉ SOUS L'EFFET D'UNE INTERACTION GÉOCHIMIQUE

A l'aide du cadre théorique précédemment défini et de l'analyse simplifiée des temps caractéristiques, nous allons nous focaliser dans la suite sur l'interaction géochimique entre saumure et sel pouvant modifier les propriétés de transport du sel. Cette discussion est conduite à travers deux origines possibles de cette interaction :

(1) les changements de solubilité de la saumure piégée dans les microfissures suite à une évolution de pression de pores et de température;

(2) une différence de potentiel chimique du sel en solution induite par l'écart de contraintes normales existant entre les joints de grains et les microfissures. Nous pensons que ce dernier micromécanisme est une cause essentielle du phénomène de "cicatrisation sous contrainte" évoqué dans le paragraphe 3.1 du chapitre I.

## 1. INTERACTION GÉOCHIMIQUE SOUS L'EFFET D'UNE VARIATION DE L'ÉQUILIBRE CHIMIQUE DE LA SAUMURE

On s'intéresse aux phénomènes de dissolution ou de cristallisation, pouvant conduire au colmatage des microfissures sous l'effet d'une variation de l'équilibre chimique de la solution consécutive d'un changement d'une des variables d'état (température ou pression de la solution). Ce phénomène conduit donc à une diminution de la porosité connectée et de la perméabilité dans le cas d'une cristallisation et à une augmentation de la porosité connectée et de la perméabilité dans le cas d'une dissolution.

### 1.1 Hypothèses et écriture des équations d'évolution des grandeurs d'état

Outre les hypothèses générales nécessaires à la construction du modèle théorique du milieu poreux, on adopte les simplifications suivantes à la suite des calculs des temps caractéristiques détaillés au paragraphe 5 du chapitre IV:
H1 La réaction chimique est instantanée à l'échelle du VER. Cette hypothèse est justifiée par les faibles valeurs des temps caractéristiques attachés aux réactions de dissolution/cristallisation comparées aux valeurs des autres temps caractéristiques. Dans le formalisme du chapitre 4, l'hypothèse H1 consiste à prendre la variable "degré d'avancement volumique de la réaction" ξ égale à zéro à tout instant.

H2 Les processus chimiques occasionnant des variations de température sont négligés (chaleur de cristallisation ou de dissolution faibles).

H3 Il n'y a pas d'influence de la diffusion hydraulique sur le transfert de chaleur (pas de convection). Cette hypothèse est justifiée pour le sel gemme massif.

H4 On suppose un découplage des effets purement chimiques (dissolution/cristallisation) et des effets Thermo-Hydro-Mécaniques (dilatation thermique des différents constituants, ouverture des microfissures due au chargement hydraulique et thermique) sur l'évolution de la porosité. L'apport de masse de saumure (soluté + solvant) m est alors donné par :

$$ \frac{m}{\rho_o^{fl}} = \frac{m^{thm}}{\rho_o^{fl}} + \phi^{chim} \quad (1) $$

$\rho_o^{fl}$ : masse volumique de la saumure dans les conditions initiales ; $\phi^{chim}$ : variation de porosité uniquement due aux effets de dissolution-cristallisation; $m^{thm}$ : apport de masse due uniquement au chargement thermo-hydro-mécanique (Palciauskas et Domenico, 1982 ; Mc Tigue, 1986 ; Coussy, 1991) :

$$ \frac{m^{thm}}{\rho_o^{fl}} = \frac{P-P_o}{M} + b\epsilon - 3\alpha_m(T-T_o) \quad (2) $$

$$ \frac{1}{M} = \frac{\phi}{K^{fl}} + \frac{b-\phi}{K^s} \quad \text{et} \quad b = 1 - \frac{K}{K^s} \quad (3) $$

P : pression interstitielle ; $P_o$ : pression interstitielle initiale ; T : température ; $T_o$ : température initiale ; $\epsilon$ : déformation volumique élastique (i.e trace des déformations élastiques), $\alpha_m$ : coefficient de dilatation linéique du fluide relativement au squelette solide; M : module de Biot; b : coefficient de Biot; $K^s$ : module de compressibilité du solide; $K^{fl}$ : module de compressibilité de la phase fluide.

Les variations de masse volumique induites par les changements de concentrations étant faibles, elles sont, par conséquent, négligées.

H5 Il existe une zone endommagée autour de l'excavation. Conformément aux mesures in situ sur le site du WIPP (figure 10 du chapitre I), la perméabilité augmente dans cette zone selon la relation :

$$ k(r) = k(r_o) \left( \frac{r}{r_o} \right)^{\alpha} $$

avec $\alpha \approx -20$ et $k(r_o)=10^{-16} m^2$ (voir également figure 8 de ce chapitre).

Les équations (1), (2) et l'équation de la conservation de la masse de saumure permettent d'écrire pour un problème unidimensionnel, l'équation de diffusion de la pression interstitielle (P). En supposant que les grandeurs d'état, la pression, la température et la déformation élastique changent peu loin de l'ouvrage souterrain, l'équation de diffusion de la pression interstitielle prend la forme suivante :
$$ \left( \frac{1}{M} + \frac{b^2}{K + 4/3\mu} \right) \frac{\partial P}{\partial t} + \left( -3\alpha_m + \frac{3\alpha^t K}{K + 4/3\mu} \right) \frac{\partial T}{\partial t} + \frac{\partial \phi^{\text{chim}}}{\partial t} = \text{div} \left[ \frac{k(\phi^{\text{chim}})}{\eta} \text{grad} P \right] \quad (4) $$

k : perméabilité intrinsèque; K : module de compressibilité non drainé; μ : module de cisaillement; P : pression du fluide;

Les hypothèses H2 et H3 permettent de simplifier l'équation de transfert de chaleur :

$$ (\rho C) \frac{\partial T}{\partial t} = \text{div} \left[ \lambda_T \underline{\text{grad}} T \right] \quad (5) $$

$\lambda_T$ : conductivité thermique du sel, $\rho$ : masse volumique du sel, $C$ : capacité calorifique massique du sel.

Pour résoudre les équations d'évolution (4) et (5), il est nécessaire d'adjoindre deux lois complémentaires portant sur les évolutions de la porosité $\phi^{\text{chim}}$ et de la perméabilité $k(\phi^{\text{chim}})$ dues aux phénomènes chimiques.

## 1.2 Écriture des lois complémentaires. Modèle de perméabilité

L'approche présentée s'inspire de travaux antérieurs relatifs à des modèles de perméabilité de roches poreuses et saturées (Dienes, 1982; Gravilenko et Guéguen, 1989; Guéguen et Palciauskas, 1992). Connaissant la distribution moyenne des microfissures et leurs caractéristiques géométriques moyennes, elle consiste à écrire les paramètres macroscopiques : la porosité $\phi$ et la perméabilité intrinsèque $k$ en fonction de la moyenne des microvariables (tableau 1).

<table><caption>Tableau 1. Relations entre paramètres macroscopiques (porosité et perméabilité) et paramètres microscopiques pour deux modèles de réseau poreux.</caption><thead><tr><th></th><th>Modèle de Plans</th><th>Modèle de Tubes</th></tr></thead><tbody><tr><th>Porosité</th><td>$\phi = \frac{2\bar{w}\bar{a}^2}{\bar{l}^3}$</td><td>$\phi = \frac{\pi\bar{r}^2\bar{d}}{\bar{l}^3}$</td></tr><tr><th>Perméabilité</th><td>$k = \frac{4\bar{w}^3\bar{a}^2}{15\bar{l}^3}$</td><td>$k = \frac{\pi\bar{d}\bar{r}^4}{32\bar{l}^3}$</td></tr><tr><th>Relation entre perméabilité et porosité</th><td>$k = \frac{\bar{w}^2\phi}{3}$</td><td>$k = \frac{\bar{r}^2\phi}{8}$</td></tr></tbody></table>

Un tel modèle permet alors d'étudier l'évolution temporelle des propriétés de transport en fonction des paramètres thermodynamiques de la roche (notamment le tenseur des contraintes, la pression de pores, et la température). On s'appuie sur un modèle de perméabilité volontairement simple basé sur les hypothèses suivantes :

H1 le phénomène de dissolution/cristallisation n'influence pas le degré de connectivité du réseau poreux.

H2 La microfissuration est isotrope et diffuse. On envisage là encore pour simplifier, deux types de distribution de microfissures. Dans un premier cas, les microfissures constituent des plans de demi-ouverture variable (w) et de longueur variables (a). Dans un deuxième cas, elles sont constituées par des tubes de longueur (d) et de rayon (r) (figure 1).
Pour les deux géométries envisagées, l'espacement entre chaque microfissure est l. Chaque microvariable w, a, d, r, et l est définie par des lois statistiques. La densité de microfissures est donc $1/\bar{l}^3$, où $\bar{l}$ est la distance moyenne entre microfissures. En l'absence d'informations précises sur la structure poreuse d'une variété donnée de sel gemme, ce modèle simple décrivant la perméabilité d'une façon approchée présente un intérêt méthodologique.

<figure><img src="image_50.png" /><figcaption>Figure 1. Modèles simplifiés de microfissures. A-modèle de plans; B- modèle de tubes</figcaption></figure>

Si l'on suppose que les moyennes des microvariables $\bar{a}, \bar{d}$ et $\bar{l}$ restent constantes, les changements d'état de saturation (c'est-à-dire les variations de la concentration à l'équilibre $\delta C_s$) de la saumure agissent respectivement sur la demi-ouverture moyenne $\bar{w}$, pour le modèle de plans et sur le rayon $\bar{r}$, pour le modèle de tubes.

### Calcul de $\delta\bar{w}$ et $\delta\bar{r}$:

A l'instant $t_0$, le volume de saumure piégée dans la microfissure est en moyenne $\bar{V}_o$; la concentration massique à l'équilibre est $C_o$ et la masse volumique de la saumure est $\rho_o$. A l'instant $t_0+dt$, le volume de saumure piégée dans la microfissure est en moyenne $\bar{V}$; la concentration à l'équilibre est $C$ et la masse volumique de la saumure est $\rho$. Entre ces deux instants, la température et la pression passent respectivement de $T_0$ à $T$ et de $P_0$ à $P$. La conservation de la masse de sel entre les deux instants considérés s'écrit :

$$ \rho_s \delta V_{sel} = C \rho \bar{V} - C_o \rho_o \bar{V}_o \quad (6) $$

$\delta V_{sel}$ : volume dissout ou cristallisé ; $\rho_s$ : masse volumique du solide.

On ne s'intéresse qu'aux variations de $\delta w$ et de $\delta r$ produites par les changements de solubilité, on prend donc $\rho=\rho_o$ (les variations de la masse volumique avec la concentration sont d'ailleurs négligeables pour les faibles variations de solubilité envisagées). En remarquant que $\bar{V}=2\bar{w}\bar{a}^2$ pour le modèle de plans et $\bar{V}=\bar{d}\pi\bar{r}^2$ pour le modèle de tubes, et en ne gardant que les termes du premier ordre, l'équation (6) devient :

#### pour le modèle de plans :

$$ \frac{\delta \bar{w}}{w_o} = \frac{\rho_o C_o}{\rho_s - \rho_o C_o} \frac{\delta C_s}{C_o} \quad (7) $$

#### pour le modèle de tubes :

$$ \frac{\delta \bar{r}}{r_o} = \frac{\rho_o C_o}{\rho_s - \rho_o C_o} \frac{\delta C_s}{2C_o} \quad (8) $$

Les variations de concentration massique à l'équilibre $\delta C_s$ sont données par une loi linéaire approchée (l'ATIG, 1986) :

$$ \delta C_s = A \delta T + B \delta P \quad (9) $$
avec $A=1,25 \ 10^{-4} \,^\circ\text{C}^{-1}$ ($0<T<50 \,^\circ\text{C}$) et $B=7,510^{-5} \,\text{Pa}^{-1}$ ($1<P<20 \,\text{MPa}$). On notera que la concentration massique à saturation $C_s$ est une fonction croissante de la température et de la pression. En remplaçant (9) dans (7) et (8), puis dans l'équation d'évolution de la porosité $\phi$ du tableau 1, on obtient pour les deux modèles la même relation :

$$ \phi^{\text{chim}} = \phi_o \frac{\rho_o}{\rho_s - \rho_o C_o} (A \delta T + B \delta P) \quad (10) $$

Compte tenu des valeurs numériques des coefficients A et B, il est facile de constater dans l'équation (10) que les variations de porosité seront plus sensibles aux changements de température de la saumure plutôt qu'aux changements de pression. On peut également remarquer que les variations de porosité sont indépendantes du modèle de réseau poreux retenu.

Les variations de perméabilité sont calculées à partir des variations de $\delta w$ (ou $\delta r$) et de $\delta \phi^{\text{tot}}$ à l'aide des relations "perméabilité-porosité" données dans le tableau 2.

Le système d'équations auquel obéissent les grandeurs d'état T et P est par conséquent :

$$ \left\{ \begin{array}{l} (\rho C) \frac{\partial T}{\partial t} = \operatorname{div}[\lambda_T \operatorname{grad} T] \\[1em] \left( \frac{1}{M} + \frac{b^2}{K + 4/3\mu} \right) \frac{\partial P}{\partial t} + \left( -3\alpha_m + \frac{3\alpha K}{K + 4/3\mu} \right) \frac{\partial T}{\partial t} + \frac{\partial \phi^{\text{chim}}}{\partial t} = \operatorname{div} \left[ \frac{k(\phi^{\text{chim}})}{\eta} \operatorname{grad} P \right] \\[1em] \phi^{\text{chim}} = \phi_o \frac{\rho_o}{\rho_s - \rho_o C_o} (A \delta T + B \delta P); \quad \delta k = \delta \left( \frac{\bar{w}^2}{3} \phi \right) (\text{mod \textit{èle plans}}); \quad \delta k = \delta \left( \frac{\bar{r}^2}{8} \phi \right) (\text{mod \textit{èle tubes}}) \end{array} \right. \quad (I) $$

Compte tenu de la nature mathématique des relations "perméabilité-porosité" utilisées, les solutions de ce système fournirent les mêmes ordres de grandeur pour les deux modèles de perméabilité envisagés. Le système d'équations (I) est résolu à l'aide d'un schéma implicite de type différences finies en géométrie cylindrique dans les deux cas présentés ci-dessous.

## 1.3 Cas d'une pression constante

On considère le cas d'un forage vertical dans un massif salifère. A une profondeur de 1000 mètres, la roche et la saumure piégée dans les interstices sont supposées être à une température de 50 °C.

A l'instant initial, on injecte dans le forage une saumure saturée en sel à la température maintenue constante et égale à 10 °C. La pression dans le puits à une profondeur de 1000 m est celle de la colonne de saumure, par conséquent égale à 12 MPa. Vu la différence de température entre la paroi du forage et le massif, la température diffuse progressivement dans le massif de sorte que la concentration à saturation de la saumure du massif diminue également. Comme la concentration à saturation baisse, la saumure se cristallise et colmate le réseau poreux: la porosité cinématique et la perméabilité vont alors diminuer. Les résultats d'un calcul transitoire, présentés sur la figure 2 reproduisent les tendances attendues.

Les deux paramètres microstructuraux $\bar{w}_o$ et $\bar{r}_o$, respectivement la demi-ouverture initiale de microfissures et le rayon initial de tubes sont calculés à partir du profil initial de perméabilité donné par la figure 1, et d'une porosité initiale prise égale à 0,1 % (cf. tableau 2). En r=a, à la paroi du forage, la température appliquée est de 10 °C. En r=b=2a, à la frontière entre la zone endommagée et le massif sain, la température est fixée à 50 °C. Cette dernière condition à la limite n'a pas de grand sens physique puisque le massif va se refroidir sur une distance plus grande que r=2a. Cette condition n'est qu'une commodité de
calcul afin d'isoler seulement l'effet géochimique de la température dans la zone endommagée, là où nous savons que la notion de milieu poreux équivalent s'applique.

On peut également comparer les variations de porosité induites par la cristallisation et les variations induites par le chargement thermique données par :

$$ \delta\phi = (b - \phi) \text{tr}\underline{\underline{\epsilon^e}} - 3\alpha_m \delta T \quad (11) $$

avec $\text{tr}\underline{\underline{\epsilon^e}} = 3\alpha_o K_o / (\lambda_o + 2\mu) \delta T$ en condition drainée, $\lambda_o$: premier coefficient de Lamé, $\mu$: module de cisaillement, $\alpha$: coefficient de dilatation thermique, $K_o$: module de compressibilité.

On a reporté sur la figure 3 les variations de porosité induites par le chargement thermomécanique et décrite par l'équation (11).

En dépit du choix volontaire des paramètres visant à majorer les variations de porosité et de perméabilité, celles-ci restent néanmoins faibles, très localisées au voisinage de la paroi du forage, et bien en deçà des variations provoquées par le seul chargement thermo-mécanique. L'effet de cristallisation sur les variations de porosité est maximum pour les plus fortes perméabilités initiales et s'annule en r=b (figure 3). Ces résultats sont indépendants du modèle de réseau poreux choisi. Les écarts calculés entre les deux modèles sont inférieures à 1%.

<figure><img src="image_51.png" /><figcaption>Figure 2. Réduction de la perméabilité en % due à la cristallisation (C) et due au chargement thermomécanique (B) (toutes les deux rapportées à la perméabilité initiale) en fonction de la distance normalisée à la paroi du forage pour un modèle de plans.</figcaption></figure>

A droite, l'ordonnée, en échelle logarithmique, correspond au profil initial de perméabilité (A) proche de celui mesuré sur le site du WIPP (figure 10 du chapitre I). Les conditions aux limites sont $P_0=P_1=12 \text{ MPa}$; $T_1=10^\circ\text{C}$ en $r=a$ en paroi (ou $R=1$, si R est la distance normalisée à la paroi); $T_0=50^\circ\text{C}$ (température initiale et température à la frontière de la zone endommagée $R=2$).

<figure><img src="image_52.png" /><figcaption>Figure 3. Réduction de la porosité normalisée due à la cristallisation et due au chargement thermomécanique en fonction de la distance normalisée à la paroi du forage.</figcaption></figure>
## 1.4 Cas d'une température constante

Dans cette deuxième simulation, on injecte une saumure possédant la température du massif (pas de diffusion thermique). A l'instant initial, on monte en pression dans le forage jusqu'à une valeur de 18 MPa de sorte que la saumure du forage s'introduise dans le massif. La concentration de la saumure à saturation dans le massif augmente progressivement du fait de la diffusion hydraulique et par conséquent, la saumure dissout la paroi du réseau poreux : la porosité cinématique et la perméabilité augmentent.

Les résultats des calculs effectués sont présentés sur la figure 4 avec les mêmes caractéristiques intrinsèques que celles de l'exemple précédent. Là encore, on retrouve les tendances attendues. L'effet de dissolution reste encore très faible en regard du cas à pression constante. Les variations de porosité et de perméabilité induites par le chargement hydro-mécanique sont supérieures, de deux ordres de grandeur, aux variations induites par dissolution.

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

Figure 4. Augmentation de la perméabilité et de la porosité due à la dissolution (%) en fonction de la distance à la paroi du forage pour un modèle de plans. A droite, l'abscisse, en échelle logarithmique, correspond au profil initial de perméabilité. Les conditions aux limites sont P₁=18 MPa et T₁=50°C en r=a en paroi (R=1, si R est la distance normalisée à la paroi); Po=12 MPa et To=T₂=50 °C (température initiale et à la frontière de la zone endommagée en R=2).

## 1.5 Conclusion partielle

Les résultats des deux calculs approchés ci-dessus mettent en évidence que les variations de la porosité et de perméabilité induites par le changement de concentration à saturation lié à une modification de la pression ou de la température sont très faibles et négligeables.

## 2. CICATRISATION SOUS CONTRAINTE

Dans cette partie, on considère le phénomène de "la cicatrisation sous contrainte" évoqué dans le paragraphe 3.1 du chapitre I. On rappelle que Cave et Sutherland (1980) et Peach (1991) ont pu observer expérimentalement une baisse différée de la perméabilité à la saumure d'un échantillon de sel saturé en saumure et soumis à un chargement isotrope.

Les conditions isothermes et drainées des essais de Peach suggèrent que l'évolution différée de la perméabilité n'est pas associée à l'évolution dans le temps d'une grandeur d'état de la saumure : la température ou la pression de pores ; c'est bien un phénomène localisé, qui résulte directement d'une sollicitation mécanique isotrope.

### 2.1 Hypothèses

Par analogie avec le phénomène de densification contrôlée par diffusion (en anglais, diffusional densification) étudié notamment par Pharr et Ashby (1983), le processus étudié
consiste en un transfert de masse depuis les zones de compression au niveau de joints de grains vers des zones avec une porosité macroscopique (microfissures engendrées par la décompression des terrains) qui seraient alors progressivement colmatées. La modélisation ici proposée possède les mêmes étapes que celles de Olivella *et al.* (1992) qui avaient étudiés ce processus sur un agrégat de sel broyé.

A l'échelle du grain, on considère le système fermé au regard du transfert de soluté. Dans ce système, il existe un assemblage de microfissures assimilées à des tubes de longueur $d$ et de rayon $r$ ou à des microfissures planes de demi-épaisseur $w$, de longueur et de largeur $a$. Ces microfissures, contrôlant la perméabilité de l'échantillon, peuvent atteindre une taille de plusieurs microns. Elles sont remplies de saumure à la pression $P$ en communication avec deux joints de grains de longueur et de largeur $\delta x$ disposés symétriquement par rapport à une microfissure. Les relations qui lient les paramètres macroscopiques, perméabilité et porosité, et les paramètres moyens microstructuraux sont aussi valables dans ce contexte (tableau 1).

Conformément aux observations de Schutjens (1991), la structure du contact des deux lèvres de chaque joint de grain est de type "île et canaux" (figure 5). Dans ce modèle, les "îles" constituent des "ponts" de contact solide/solide, là où s'exercent des contraintes normales. Les "canaux" entourant les "îles", sont remplis de saumure et permettent le transport en solution des substances dissoutes. L'épaisseur du joint de grains serait de l'ordre de plusieurs dizaines de nanomètres.

L'utilisation d'un tel modèle pour des applications quantitatives n'est pas aisée : le rapport entre les surfaces de contact solide/solide et les surfaces "mouillées" qui contrôle la répartition des contraintes le long d'un joint de grain est difficile à estimer. Ainsi, dans un souci de simplicité, on admettra qu'il existe sur les joints de grain une contrainte normale moyenne $\sigma$, correspondant à la contrainte hydrostatique appliquée à l'échantillon.

La cicatrisation est supposée être contrôlée par la diffusion en solution. Cette hypothèse est largement admise dans le cas des processus de déformation par dissolution-cristallisation dans le système sel-saumure (Spiers *et al.*, 1990; Schutjens, 1991). En outre, on considérera que la solution dans les microfissures et aux joints de grains est en équilibre chimique à chaque instant.

<figure><img src="image_54.png" /><figcaption>Figure 5. Les deux échelles d'espace à considérer dans le processus de cicatrisation.</figcaption></figure>
## 2.2 Écriture des équations du modèle

Le flux massique de soluté M qui est transporté dans le joint obéit à la loi de Fick :

$$ \underline{M} = -\rho^{fl} D \underline{\underline{\text{grad} C_s}} \qquad (12) $$

$C_s$ : concentration massique en sel de la solution, $D$ : coefficient de diffusion du soluté, $\rho^{fl}$ : masse volumique de la saumure ici supposé constante.

Le potentiel chimique du soluté $\mu_s$ se définit, de manière générale, par :

$$ \mu_s = \mu_s^0 + RT \ln(\gamma C_s) \quad \text{ou} \quad \underline{\underline{\text{grad} C_s}} = \frac{C_s}{RT} \underline{\underline{\text{grad} \mu_s}} \qquad (13) $$

$R$ : constante des gaz parfaits ; $T$ : température ; $\gamma$ : coefficient d'activité du soluté supposé constant en solution.

En remplaçant cette relation dans la loi de Fick, nous obtenons :

$$ \underline{M} = -\rho^{fl} \frac{D C_s}{R T} \underline{\underline{\text{grad} \mu_s}} \qquad (14) $$

Si on suppose que $\text{grad}\mu_s \approx \delta\mu_s/\delta x$ ($\delta x$ : largeur et longueur du joint de grains), on peut réécrire la norme du flux massique M de la façon suivante :

$$ M = \rho^{fl} \frac{D C_s}{R T} \frac{|\delta \mu_s|}{\delta x} \qquad (15) $$

Si toute la masse apportée par les joints de grain précipite sur les parois de la microfissure au cours d'un petit instant $\delta t$, on obtient alors :

$$ 2\bar{w} \bar{a} \delta \bar{a} \rho_s = -M \delta t \delta x^2 $$

$\rho_s$ : masse volumique de sel.

On aura noté que lorsque t augmente, la longueur moyenne de microfissure $\bar{a}$ diminue. L'épaisseur moyenne $\bar{w}$ est supposée inchangée. On peut écrire encore :

$$ \frac{\partial(\bar{a}^2)}{\partial t} = - \frac{M \delta x^2}{\bar{w} \rho_s} \qquad (16) $$

En utilisant (15) et $\delta\mu_s>0$ :

$$ \frac{\partial(\bar{a}^2)}{\partial t} = - \frac{\delta x^2}{\bar{w} \rho_s} \rho^{fl} \frac{D C_s}{R T} \frac{\delta \mu_s}{\delta x} \qquad (17) $$

Il convient d'expliciter la différence de potentiel chimique $\delta\mu_s$, le moteur du transfert de masse entre le joint de grains et la microfissure. La condition d'équilibre, en un point de la surface solide-fluide, est traduite par l'égalité du potentiel chimique du soluté en solution ($\mu_f$) et le potentiel chimique du solide intact sous contrainte ($\mu_s$) (Paterson, 1973; Guéguen et Palciauskas, 1992) :

$$ \mu_f = \mu_s = f_s + \sigma_n V_m $$

$f_s$ : énergie libre molaire, $\sigma_n$ : contrainte normale à la surface solide et $V_m$ : volume molaire du solide.
Par conséquent, si le moteur du transfert de masse résulte principalement du gradient de contrainte normale entre le joint de grains et la microfissure, alors la différence de potentiel chimique du soluté entre ces deux zones peut être approximée par :

$$ \delta\mu_s = \delta\sigma_n V_m \quad (18) $$

On a ainsi négligé dans le bilan énergétique les apports dus aux effets de surface et aux processus de déformation des grains. $\delta\sigma_n$ est l'écart de contrainte normale entre le joint et la microfissure. On peut prendre $\delta\sigma_n = \sigma_I - P$ avec $\sigma_I$: la contrainte principale majeure supposée égale à la contrainte lithostatique $\sigma$. Cette dernière égalité doit être pondérée par une fonction telle que si la porosité $\phi \to 0$ alors $\delta\mu_s \to 0$. De manière simple, on prend cette fonction égale à $\phi$ :

$$ \delta\mu_s = \phi(\sigma - P)V_m \quad (19) $$

En remplaçant (19) dans (17), on obtient :

$$ \frac{\partial(\bar{a}^2)}{\partial t} = - \frac{\rho^{fl} V_m DC_s}{\rho_s RT} \frac{\delta x}{\bar{w}} \phi(\sigma - P) \quad (20) $$

En procédant au même raisonnement pour le modèle de tubes, on obtient :

$$ \frac{\partial \bar{S}}{\partial t} = - \frac{\rho^{fl} V_m DC_s}{\rho_s RT} \frac{2 \delta x}{\bar{d}} \phi(\sigma - P) \quad (21) $$

Avec $\bar{S} = \pi\bar{r}^2$. On utilise encore les relations du tableau 1 qui lient le paramètre macroscopique $\phi$ et le paramètre microstructural $\bar{a}$ pour le modèle de plans et $\bar{r}$ pour le modèle de tubes. L'espacement moyen entre chaque microfissure $\bar{l}$ et la longueur moyenne $\bar{d}$ sont supposés rester inchangés.

$$ \text{modèle de plans :} \quad \frac{\partial \phi}{\partial t} = 2 \frac{\bar{w}}{\bar{l}^3} \frac{\partial (\bar{a}^2)}{\partial t} ; \quad \text{modèle de tubes :} \quad \frac{\partial \phi}{\partial t} = \frac{2\pi\bar{r}\bar{d}}{\bar{l}^3} \frac{\partial \bar{r}}{\partial t} \quad (22) $$

Dans cette écriture, la porosité des joints de grains et son évolution dans le temps sont négligées. En utilisant (20) et (21), on a alors pour les deux modèles de perméabilité :

$$ \frac{\partial \phi}{\partial t} = -2 \frac{\rho^{fl} V_m DC_s \delta x}{\rho_s RT \bar{l}^3} \phi(\sigma - P) \quad (23) $$

Dans la relation (23), la variation de porosité ne dépend plus que des paramètres microstructuraux $\delta x$, la longueur du joint de grain, et $1/\bar{l}$, la densité de microfissures. Elle est indépendante de la nature géométrique (plans ou tubes) des modèles de perméabilité. Introduisons $\alpha$ tel que :

$$ \alpha = 2 \frac{\rho^{fl} V_m DC_s \delta x}{\rho_s RT \bar{l}^3} (\sigma - P) \quad (24) $$

L'équation (23) peut être réécrite de la manière suivante :

$$ \frac{\partial \phi}{\partial t} = -\alpha \phi \quad (25) $$
Les données expérimentales des essais de Peach autorisent à supposer que le coefficient $\alpha$ reste constant au cours du temps. Si on note $\phi_0$ la porosité à l'instant initial, l'équation (25) se résout simplement :

$$ \phi = \phi_0 \exp(-\alpha t) \qquad (26) $$

La loi d'évolution de la perméabilité est alors obtenue en utilisant le résultat (26) et les relations du tableau 1 :

$$ \text{modèle de plans :} \quad k = k_o \exp(-\alpha t) \qquad \text{modèle de tubes :} \quad k = k_o \exp(-2\alpha t) \qquad (27) $$

$$ \text{avec } k_o = \frac{\bar{w}^2}{3} \phi_0 \qquad \text{avec } k_o = \frac{\bar{l}^3}{8\pi d} \phi_0^2 $$

## 2.3 Comparaison avec les données expérimentales

Comparons les données expérimentales de Peach et les résultats obtenus par les équations (27), en prenant les valeurs numériques suivantes :

$p^{fl}=1,210^3 \text{ kg.m}^{-3}$; $p^s = 2,16 \ 10^3 \text{ kg.m}^{-3}$; $V_m=2.69 \ 10^{-5} \text{ m}^3.\text{mol}^{-1}$; $C_s=0,265 \text{ kg.kg}^{-1}$; $R=8,314 \text{ J.mol}^{-1}.\text{K}^{-1}$; $D=10^{-9} \text{ m}^2.\text{s}^{-1}$; $\phi_0=1,5 \text{ %}$, $k=3,2 \ 10^{-17} \text{ m}^2$ pour $P_{eff}=8,5 \text{ MPa}$; $k=4.65 \ 10^{-18} \text{ m}^2$ pour $P_{eff}=18,5 \text{ MPa}$.

Le paramètre Peff est une pression dite effective donnée par Peach : par exemple $P_{eff}=\sigma-Pf=8,5 \text{ MPa}$, où la contrainte lithostatique $\sigma$ est égale à 10 MPa et la pression de pores Pf est égale à 1,5 MPa; ou encore $P_{eff}= \sigma-Pf=18,5 \text{ MPa}$, où la contrainte lithostatique $\sigma$ est égale à 20 MPa et la pression de pores Pf est égale à 1,5 MPa

Pour le modèle de plans, la meilleure concordance entre les données expérimentales et les résultats des calculs est obtenue pour une valeur du rapport $\delta x / \bar{l}^3$ égale à 24000 (figure 6). Ce qui correspond, par exemple pour l'essai à $P_{eff}=8,5 \text{ MPa}$ aux valeurs de $\delta x=2,5 \text{ mm}$ et $\bar{l}=5 \text{ mm}$. A l'aide des relations (27) et des valeurs de $k_o$ et $\phi_0$, l'épaisseur moyenne de microfissures $\bar{w}$ devient alors égale à $10^{-1} \mu\text{m}$.

Pour le modèle de tubes, afin d'avoir la même concordance entre valeurs calculées et valeurs expérimentales, la valeur du rapport $\delta x / \bar{l}^3$ est 35000. Pour un choix de $\delta x=2,5 \text{ mm}$, on obtient une valeur de $\bar{l}$ égale à $2,68 \ 10^{-4} \text{ m}$. A l'aide des relations (27) et des valeurs de $k_o$ et $\phi_0$, on calcule $\bar{d}$, une valeur de longueur moyenne de tubes égale à 8 800 m (essai à Peff=8,5 MPa) et 62 000 m (essai à Peff= $\sigma-Pf=18,5 \text{ MPa}$). Ces résultats suggèrent que le modèle de tubes n'est pas adapté pour représenter simplement la structure poreuse du sel.

Les différences néanmoins constatées sur la figure 6 sont probablement dues aux déformations viscoplastiques des grains qui n'ont pas été prises en compte dans ce calcul. Or ces déformations apparaissent dans le cas du sel gemme, matériau dépourvu de tout seuil d'écoulement viscoplastique, et se poursuivent donc tant que les contraintes déviatoriques règnent localement du fait de la rugosité des microfissures. On conçoit que le fluage des grains entraîne progressivement une cicatrisation de microfissures, d'ailleurs prouvée expérimentalement par Peach (1991) et par Stormont et Daemen (1992) au cours d'essais de fluage sur des échantillons de sel sec. Ce phénomène d'origine mécanique se superpose à une cicatrisation d'origine mécano-chimique et permet de penser que la perméabilité du sel, à l'échelle de temps géologique, doit être pratiquement nulle loin des ouvrages souterrains.
<figure><img src="image_55.png" /><figcaption>Figure 6. Evolution de la perméabilité. Comparaison entre données expérimentales et modèle.</figcaption></figure>

Peff=σ-Pf=8,5 MPa, la contrainte de confinement σ est égale à 10 MPa et la pression de pores Pf est égale à 1,5 MPa.

Peff=σ-Pf=18,5 MPa, la contrainte de confinement σ est égale à 20 MPa et la pression de pores Pf est égale à 1,5 MPa.

# PARTIE II. INFLUENCE DE LA PERMÉABILITÉ SUR LA SÛRETÉ DES STOCKAGES DANS LE SEL

## 1. PHÉNOMÈNES INDUISANT UNE PERTE D'ÉTANCHÉITÉ DU STOCKAGE

### 1.1 Endommagement et fissuration d'origine mécanique

*a- Rappel sur les critères d'endommagement et de rupture*

Les ouvrages profonds sont soumis à un état de contrainte pouvant microfissurer la roche au voisinage de l'ouvrage (voir chapitre I), ou encore pouvant, en dépassant la résistance à la rupture (en compression ou en traction) du matériau, occasionner une fracturation de la paroi.

Dans le premier cas, la microfissuration donne naissance à une zone de roche endommagée autour de l'ouvrage (en anglais *Disturbed Rock Zone*) présentant des propriétés hydrogéologiques particulières : augmentation drastique de la perméabilité depuis la paroi, désaturation en saumure des pores, drainage depuis une zone située à proximité de l'ouvrage vers la paroi de celui-ci (figure 10 du chapitre I).

Un indicateur intéressant de cette microfissuration est le déviateur correspondant au seuil d'accélération de la variation de volume, mis en évidence par des essais de compression triaxiale de révolution. L'assimilation de cette limite de linéarité à l'initiation de l'endommagement est en particulier justifiée dans le cas du sel gemme, dont le comportement viscoplastique n'est associé à aucune variation de volume significative (Munson, 1981 ; Langer, 1981 ; Ghoreychi, 1990 ; Thorel, 1994). Toute dilatance irréversible (anélastique) du sel gemme est donc attribuée à une microfissuration de celui-ci.

Dans le deuxième cas, le dépassement de la résistance à la rupture conduit à l'apparition de fissures macroscopiques (Cook et Rogenthen, 1991) et du phénomène d'écaillage en paroi de galerie.

Aussi, grâce à une mesure fine de variation de volume global des échantillons soumis aux essais triaxiaux de laboratoire, il est possible de présenter, dans le plan de
contrainte moyenne-contrainte déviatorique, deux critères caractérisant deux états d'équilibre distincts (Figures 7, 9 et 10 du chapitre II) :

. Le seuil d'endommagement correspondant à l'initiation d'une fissuration microscopique et donc à une augmentation de perméabilité plus ou moins diffuse et isotrope.

. Le critère de rupture identifié à partir de la résistance maximale. Ce critère correspond approximativement à l'initiation d'une fissuration macroscopique (fracturation). Au-delà de cette frontière, la continuité macroscopique de la matière n'étant pas assurée, la perméabilité ne peut être considérée comme isotrope et doit être représentée par un tenseur au lieu d'un scalaire (cf. chapitre I).

<figure><img src="image_56.png" /><figcaption>Figure 7 : Critères d'endommagement et de rupture du sel gemme (Sel des Mines de Potasse d'Alsace).</figcaption></figure>

Les critères d'endommagement permettent une analyse qualitative de l'évolution de la perméabilité en fonction du tenseur de contrainte (parties sphérique et déviatorique), mais ne sont pas de nature à décrire quantitativement cette évolution. Cette question relève de la loi d'évolution de l'endommagement, corrélée à celle de la perméabilité.

Une telle description servirait au calcul hydromécanique pour tenir compte de la zone endommagée autour des ouvrages souterrains et de son influence sur le régime hydraulique du champ proche.

On peut penser que des relations analogues à celles utilisées par des hydrogéologues (relations reliant la perméabilité et la porosité), peuvent être déterminées expérimentalement pour des matériaux endommagés. La porosité serait alors liée à la trace du tenseur de déformation irréversible et probablement à la partie sphérique de ce tenseur (dilatance ou contractance) caractérisant l'endommagement.

S'agissant des critères d'endommagement et de rupture présentés sur la figure 7, rappelons, comme l'attestent les résultats de la partie I du chapitre II, que les critères obtenus à l'état sec sont significativement différents de ceux déterminés sur le sel
directement en contact avec le fluide sous pression. Il y a donc deux états d'équilibre bien distincts suivant que la stabilité de l'ouvrage souterrain est analysée avec les critères du sel sec ou du sel humide. Cette question fait l'objet du paragraphe qui suit.

## b- Fissuration et fracturation de la paroi d'une cavité saline

On considère le cas d'une cavité souterraine située à la profondeur H où règne un état de contraintes isotrope $\sigma_0$. Le massif est constitué de sel gemme assimilé à un milieu poro-élasto-viscoplastique obéissant à la loi rhéologique de Norton-Hoff. La pression de pores initiales est $P_0$ à la profondeur moyenne de la cavité. Cette dernière est creusée par lessivage. Elle est soumise à une pression de saumure $P_1$.

L'état d'équilibre de la caverne est discuté dans deux cas extrêmes :

a- à court terme en supposant un creusement instantané de la caverne, ce qui revient à admettre un comportement poro-élastique du massif.

b- à long terme, en considérant un état stationnaire de tous les processus, ce qui signifie que le massif est en équilibre thermique avec la caverne (état isotherme) alors que le champs de pression de pores atteint son régime permanent. Le comportement mécanique, lui aussi, est caractérisé par une évolution stationnaire de la déformation du massif et donc de la convergence de la paroi de la cavité tandis que le champs de contraintes tend vers un état asymptotique caractéristique de la loi viscoplastique de Norton-Hoff (Halphen, 1978; Pouya, 1991b).

Afin de cerner l'influence de la géométrie de la caverne, là encore, deux configurations distinctes sont étudiées : la sphère et le cylindre. Les deux structures se caractérisent par un rayon unité alors que le champs de contraintes et de pression de pores initiaux sont vérifiés sur un rayon extérieur b. Ce dernier se trouve suffisamment loin de la paroi de la cavité et permet de se rapprocher du milieu infini. Ce choix n'a d'autre but que de permettre l'existence d'un régime hydraulique permanent dans le cas du cylindre (ce dernier n'existe pas pour un cylindre dans un milieu infini, contrairement à la sphère).

Pour analyser la tenue mécanique de la paroi de ces cavités à court et à long terme, nous avons comparé l'état des contraintes calculés en poro-élasticité et en poro-viscoplasticité avec les critères d'endommagement et de rupture. Ces derniers correspondent soit à un sel sec (critères obtenus sur des échantillons gainés) soit directement en contact avec de la saumure (échantillons non gainés).

Dans ce but, les trois contraintes principales régnant à la paroi des cavités sont d'abord calculées. Il s'agit des contraintes radiales ($\sigma_{rr}$), orthoradiale ($\sigma_{\theta\theta}$) et axiale ($\sigma_{zz}$) ou ($\sigma_{\phi\phi}$) pour les cavités cylindrique et sphérique respectivement.

Trois variables scalaires représentant trois invariants du tenseur des contraintes sont ensuite déterminées :

- la contrainte moyenne : $\sigma_m = \frac{1}{3} \text{tr} \underline{\underline{\sigma}}$;

- la contrainte déviatorique de Von Mises $\sigma_{eq} = \sqrt{\frac{3}{2} \underline{s} \cdot \underline{s}}$ ($\underline{s}$ : tenseur déviateur des contraintes);

- une variable adimensionnelle parfois nommée géométrie des contraintes $J_m$ donnée par :
$$J_m = \sin 3\theta = \frac{m(9-m^2)}{(3+m^2)^{3/2}}$$

$\theta$: angle de Lode; $m$: paramètre de Lode :

$$m = \frac{3s_{II}}{s_{I} - s_{III}}$$

$s_I$, $s_{II}$ et $s_{III}$ : contraintes déviatoriques principales de s$_{II}$ où $s_I < s_{II} < s_{III}$.

La variable $J_m$ permet de vérifier si la paroi de la cavité se trouve plutôt en compression ou en extension. En effet, dans un triaxial de révolution :

$$J_m = \frac{\sigma_1 - \sigma_3}{|\sigma_1 - \sigma_3|}$$

$\sigma_1$ : contrainte axiale ; $\sigma_3$ : contrainte radiale. Les contraintes en compression sont négatives.

Ainsi, pour un essai en compression (resp. en extension), le paramètre $J_m$ est égal à -1 (resp. +1).

Les tableaux 2a et 2b rassemblent les différentes composantes du tenseur de contraintes et des variables liées aux invariants de contraintes.

La comparaison entre l'état de contraintes à la paroi et les critères d'endommagement et de rupture est présentée dans les tableaux 3a et 3b. Afin de faciliter la confrontation des deux types de critères, les structures sont soumises à une pression de fluide en paroi qui correspond à la colonne de fluide à la profondeur considérée.

La comparaison permet de constater :

a- les risques d'endommagement et de rupture dépendent fortement des états initiaux de contraintes, de la pression de pores mais aussi bien sûr de la pression du fluide en paroi qui joue le rôle bénéfique d'un soutènement ;

b- les risques d'endommagement sont plus précoces en extension qu'en compression (29 MPa par rapport à 32 MPa dans le tableau 3a). Ils sont d'autant plus importants que le déviateur est fort à la paroi et que la contrainte moyenne y est faible. Ces conditions sont en général réunies à court terme lorsqu'une faible pression de fluide s'exerce sur la paroi. A long terme, la relaxation de l'état de contraintes résultant du comportement viscoplastique du sel gemme entraîne une diminution du déviateur à la paroi, ce qui améliore la tenue mécanique de la cavité;

c- les risques d'endommagement et de rupture sont nettement accrus lorsque le critère humide est considéré, c'est à dire lorsqu'est prise en compte la perméabilité initiale du milieu ;

d- le fait que l'état des contraintes en paroi des cavités dans l'analyse poro-élastique n'atteint jamais le critère de rupture du sel sec utilisé doit être considéré avec prudence car ce résultat tient essentiellement à l'hypothèse admise dans l'interprétation des résultats des essais triaxiaux axisymétriques. Ceux-ci ont été généralisés à trois dimensions en considérant que l'écart entre les contraintes axiale et latérale pouvait être confondu avec la

contrainte déviatorique équivalente au sens de Von Mises ($\sigma_{eq} = \sqrt{\frac{3}{2}\frac{s_{II}}{s_{I}}}$).
<table><caption>Tableau 2b. Composantes du tenseur des contraintes et des invariants en poro-vicoplasticité (cf. annexe V)</caption><thead><tr><th></th><th>Cavité cylindrique</th><th>Cavité sphérique</th></tr></thead><tbody><tr><td>σrr</td><td>σrr = -P<sub>1</sub> + (σ<sub>o</sub> + P<sub>1</sub>)(1 - 1/r<sup>2/n</sup>)</td><td>σrr = -P<sub>1</sub> + (σ<sub>o</sub> + P<sub>1</sub>)(1 - 1/r<sup>3/n</sup>)</td></tr><tr><td>σθθ</td><td>σθθ = σ<sub>o</sub> + (σ<sub>o</sub> + P<sub>1</sub>)(2/n - 1)(1/r<sup>2/n</sup>)</td><td>σθθ = σ<sub>o</sub> + (σ<sub>o</sub> + P<sub>1</sub>)(3/2n - 1)(1/r<sup>3/n</sup>)</td></tr><tr><td>σ<sub>zz</sub> ou σ<sub>φφ</sub></td><td>σ<sub>zz</sub> = σ<sub>o</sub> + 2v<sub>o</sub>(σ<sub>o</sub> + P<sub>1</sub>)(1/n - 1)(1/r<sup>2/n</sup>) - b(1 - 2v<sub>o</sub>)(P<sub>1</sub> - P<sub>o</sub>)</td><td>σ<sub>φφ</sub> = σ<sub>θθ</sub></td></tr><tr><td>σ<sub>m</sub> en paroi</td><td>σ<sub>m</sub> = 1/3[(1 + 2/n) + 2v<sub>o</sub>(1/n - 1)](σ<sub>o</sub> + P<sub>1</sub>) - P<sub>1</sub> - b/3(1 - 2v<sub>o</sub>)(P<sub>1</sub> - P<sub>o</sub>)</td><td>σ<sub>m</sub> = -P<sub>1</sub> + σ<sub>o</sub> + P<sub>1</sub>/n</td></tr><tr><td>σ<sub>eq</sub> en paroi</td><td>σ<sub>eq</sub> = 3/2(s<sub>r</sub><sup>2</sup> + s<sub>θ</sub><sup>2</sup> + s<sub>z</sub><sup>2</sup>)<br>s<sub>r</sub> = -1/3[(1 + 2/n) - 2v<sub>o</sub>(1 - 1/n)](σ<sub>o</sub> + P<sub>1</sub>) + 1/3b(1 - 2v<sub>o</sub>)(P<sub>1</sub> - P<sub>o</sub>)<br>s<sub>θ</sub> = -1/3[(1 - 4/n) - 2v<sub>o</sub>(1 - 1/n)](σ<sub>o</sub> + P<sub>1</sub>) + 1/3b(1 - 2v<sub>o</sub>)(P<sub>1</sub> - P<sub>o</sub>)<br>s<sub>z</sub> = -2/3(1 - 2v<sub>o</sub>)(1/n - 1)(σ<sub>o</sub> + P<sub>1</sub>) - 2/3b(1 - 2v<sub>o</sub>)(P<sub>1</sub> - P<sub>o</sub>)</td><td>σ<sub>eq</sub> = -3/2(σ<sub>o</sub> + P<sub>1</sub>)/n</td></tr><tr><td>J<sub>m</sub></td><td>J<sub>m</sub> > 0<br>P<sub>1</sub> > P<sub>o</sub> (n<4)</td><td>1</td></tr></tbody></table>
<table><caption>Tableau 3a : Analyse poro-élastique des cavités cylindriques et sphériques. Les constantes C<sup>e</sup> et C<sup>r</sup> sont les cohésions attachées respectivement aux critères d'endommagement de rupture (de type Tresca), toutes les deux obtenues en condition humide.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Critère en condition sèche (pas de perméabilité)</th><th colspan="2">Critère en condition humide</th></tr><tr><th>Cavité cylindrique</th><th>Cavité sphérique</th><th>Cavité cylindrique</th><th>Cavité sphérique</th></tr></thead><tbody><tr><td>Évolution du déviateur σ<sub>eq</sub> avec la profondeur (contrainte initiale σ<sub>o</sub>)</td><td>σ<sub>eq</sub>=√3(σ<sub>o</sub>+P<sub>1</sub>)</td><td>σ<sub>eq</sub>=-3/2(σ<sub>o</sub>+P<sub>1</sub>)</td><td>σ<sub>eq</sub>=√3(σ<sub>o</sub>+P<sub>1</sub>)</td><td>σ<sub>eq</sub>=-3/2(σ<sub>o</sub>+P<sub>1</sub>)</td></tr><tr><td>Seuil d'endommagement</td><td>σ<sub>eq</sub>=(-σ<sub>o</sub>/0.04)<sup>1/2</sup></td><td>σ<sub>eq</sub>=(-σ<sub>o</sub>/0.05)<sup>1/2</sup></td><td>σ<sub>eq</sub>=2C<sup>e</sup> (C<sup>e</sup>=5 MPa)</td><td>σ<sub>eq</sub>=2C<sup>e</sup> (C<sup>e</sup>=5 MPa)</td></tr><tr><td>Seuil d'endommagement en paroi</td><td>32 MPa (≈1830 m)</td><td>29 MPa (≈1950 m)</td><td>10 MPa (≈580 m)</td><td>10 MPa (≈660 m)</td></tr><tr><td>Critère de rupture</td><td>σ<sub>eq</sub>=-1.43σ<sub>o</sub>+10.3</td><td>σ<sub>eq</sub>=-1.5σ<sub>o</sub>+3</td><td>σ<sub>eq</sub>=2C<sup>r</sup> (C<sup>r</sup>=15 MPa)</td><td>σ<sub>eq</sub>=2C<sup>r</sup> (C<sup>r</sup>=15 MPa)</td></tr><tr><td>Critère de rupture en paroi</td><td>jamais</td><td>jamais</td><td>30 MPa (≈1730 m)</td><td>30 MPa (≈2000 m)</td></tr></tbody></table>

<table><caption>Tableau 3b : Analyse poro-viscoplastique des cavité s cylindriques et sphériques. Les constantes C<sup>e</sup> et C<sup>r</sup> sont les cohésions attachées respectivement aux critères d'endommagement de rupture (de type Tresca), toutes les deux obtenues en condition humide.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Critère en condition sèche (pas de perméabilité)</th><th colspan="2">Critère en condition humide</th></tr><tr><th>Cavité cylindrique</th><th>Cavité sphérique</th><th>Cavité cylindrique</th><th>Cavité sphérique</th></tr></thead><tbody><tr><td>Évolution du déviateur σ<sub>eq</sub> avec la profondeur (contrainte initiale σ<sub>o</sub>)</td><td>-1/3(P<sub>1</sub> + σ<sub>o</sub>) × √(7 + 16ν<sub>o</sub>(ν<sub>o</sub> - 1)) (n=3)</td><td>σ<sub>eq</sub>=-3/2n(σ<sub>o</sub>+P<sub>1</sub>)</td><td>-1/3(P<sub>1</sub> + σ<sub>o</sub>) × √(7 + 16ν<sub>o</sub>(ν<sub>o</sub> - 1)) (n=3)</td><td>σ<sub>eq</sub>=-3/2n(σ<sub>o</sub>+P<sub>1</sub>)</td></tr><tr><td>Seuil d'endommagement</td><td>σ<sub>eq</sub>=(-σ<sub>o</sub>/0.04)<sup>1/2</sup></td><td>σ<sub>eq</sub>=(-σ<sub>o</sub>/0.05)<sup>1/2</sup></td><td>σ<sub>eq</sub>=2C<sup>e</sup> (C<sup>e</sup>=5 MPa)</td><td>σ<sub>eq</sub>=2C<sup>e</sup> (C<sup>e</sup>=5 MPa)</td></tr><tr><td>Seuil d'endommagement en paroi</td><td>86 MPa (≈13 600m)</td><td>88 MPa (≈17 600 m)</td><td>10 MPa (≈1570 m)</td><td>10 MPa (≈2000 m)</td></tr><tr><td>Critère de rupture</td><td>σ<sub>eq</sub>=-1.43 σ<sub>o</sub> +10.3</td><td>σ<sub>eq</sub>=-1.5σ<sub>o</sub>+3</td><td>σ<sub>eq</sub>=2C<sup>r</sup> (C<sup>r</sup>=15 MPa)</td><td>σ<sub>eq</sub>=2C<sup>r</sup> (C<sup>r</sup>=15 MPa)</td></tr><tr><td>Critère de rupture en paroi</td><td>jamais</td><td>jamais</td><td>30 MPa (≈4720 m)</td><td>30 MPa (≈6000 m)</td></tr></tbody></table>
## 1.2 Fissuration par traction résultant d'un refroidissement

Le massif encaissant les ouvrages de stockage peut être soumis à un refroidissement rapide ou lent suivant la nature du stockage. Ainsi, l'aérage de l'infrastructure d'un stockage, l'injection d'un fluide à température ambiante dans une caverne souterraine, la détente "adiabatique" d'une caverne résultant d'une décompression brutale du gaz, ou le refroidissement de déchets nucléaires dans la phase post-paroxysme thermique, n'occasionnent pas un refroidissement aussi rapide et intense que celui d'un stockage cryogénique. Toutefois, le phénomène de refroidissement bien que lent peut aboutir, dans le cas d'un milieu viscoplastique tel que le sel gemme, à des contraintes de traction pouvant conduire à une rupture des parois des cavités (Ould Amy et Ghoreychi, 1992).

Ce phénomène a été déjà observé in situ lors d'une expérimentation dans les Mines de Potasse d'Alsace (Ghoreychi, 1991) : une coupure brutale de chauffe a provoqué des fractures ouvertes de 2 mm, le long de plusieurs sondages soumis ou non à une pression de remblai de sel broyé.

La rupture par refroidissement des ouvrages dans le sel s'explique par le fait qu'après creusement de cavités souterraines, les contraintes déviatoriques se relaxent en raison de la viscoplasticité du massif. Tout refroidissement du massif conduit alors, à court terme, à des contraintes orthoradiales (et aussi axiales pour des galeries) de traction, le massif réagissant instantanément en élasticité. La variation de contrainte ($\Delta\sigma$) liée à une variation de température $\Delta T$ peut être estimée à $\Delta\sigma=E\alpha\Delta T$, soit de l'ordre de 1MPa/°C pour le sel gemme (typiquement $E = 25000$ MPa et $\alpha = 4.10^{-5}\text{C}^{-1}$).

Cette valeur, relativement élevée, est à comparer à la résistance à la traction uniaxiale du milieu qui est, en général, faible pour des géomatériaux. Elle est de l'ordre de 2MPa pour le sel gemme.

Il convient de préciser que dans la pratique, pour une cavité saline, les contraintes de traction résultant d'un refroidissement sont plus faibles que celles estimées en élasticité.

En effet, la relaxation de contrainte dans le sel est d'autant plus rapide que le déviateur est élevé. Les contraintes de traction règnent donc dans une zone peu étendue autour de l'ouvrage et tendent progressivement vers un état de compression. Cette transition est cependant ralentie par la viscosité du sel qui croît à basse température.

On conçoit que l'amplitude et l'extension de la zone fracturée par refroidissement doivent être étudiées en comparant deux temps caractéristiques : celui de la diffusion thermique et la constante de temps de relaxation de l'état de contrainte au voisinage de l'ouvrage.

## 1.3 Fracturation hydraulique résultant d'une pression capillaire

Un milieu géologique peut être perméable à un liquide ou à un gaz, mais demeurer imperméable au gaz si le milieu est préalablement saturé par le liquide. Cette situation est rencontrée lorsque la pression de gaz est inférieure à un certain seuil, appelé "pression de seuil" (de l'anglais threshold pressure), liée au phénomène capillaire.

Selon la loi de Laplace, cette pression est proportionnelle à la tension superficielle entre le liquide et le gaz, et inversement proportionnelle au rayon de courbure du ménisque formé par des bulles de gaz dans le capillaire. Or, le rayon de courbure dépend lui-même de dimensions des pores. De ce fait, toute chose identique par ailleurs, la pression de seuil est d'autant plus élevée que la taille des pores est petite. Cette dernière varie du nanomètre ($10^{-9}$m) au micromètre ($10^{-6}$m) pour le sel gemme. Compte tenu d'une valeur de tension superficielle de 0,07 N/m, le seuil de pression estimé appartient à l'intervalle $0,3-300$ MPa.
En dessous de la pression de seuil, le sel n'étant pas perméable au gaz, la pression augmente dans la cavité de stockage consécutivement à un dégagement de gaz provenant de déchets (du fait de la radiolyse et/ou de la corrosion des conteneurs de déchets radioactifs), ou à une convergence progressive de la paroi d'une caverne de stockage de gaz fermée.

Au cas où la pression de seuil dépasserait l'état de contrainte initial dans le massif, une fracturation hydraulique de la paroi pourrait se produire. Sa propagation instable pourrait compromettre la sûreté à long terme du stockage. Ce scénario est étudié aux Etats-Unis dans le cadre du projet "WIPP" de stockage de déchets radioactifs au Nouveau Mexique (Webb et Davies, 1993).

## 2. ÉVOLUTION D'UNE CAVERNE ABANDONNÉE DANS UN MASSIF SALIFÈRE FAIBLEMENT PERMÉABLE

Dans cette partie, on s'intéresse aux conséquences de la présence d'une faible perméabilité du massif salifère sur l'évolution d'une structure souterraine particulière : une caverne abandonnée et remplie d'un liquide (hydrocarbure ou saumure). Cette faible perméabilité "en grand" peut avoir plusieurs origines : endommagement mécanique du champs proche ; présence de lits d'impuretés constituant des chemins préférentiels d'écoulement, etc..

Dans une première partie, on se focalise sur l'évolution de la caverne à court terme (sur quelques mois). Puis, après avoir discuté de la stabilité d'une zone endommagée éventuellement présente autour du stockage, l'évolution à moyen terme du système cavité-massif (sur quelques dizaines d'années) est examinée.

### 2.1 Évolution à court terme

L'évolution spontanée d'une caverne fermée est gouvernée à court terme par les phénomènes suivants (Bérest et al, 1979; Boucly et Legreneur, 1980 ; Bérest, 1989; Brouard et Bérest, 1995):

- le fluage ; le sel se comporte comme un fluide visqueux : un faible écart à la distribution hydrostatique des contraintes provoque un écoulement visqueux du sel. Ce dernier engendre le fluage de la paroi de la caverne ;

- l'expansion thermique ; le liquide présent dans la caverne est généralement plus froid que le massif encaissant. Si la caverne est fermée, l'expansion thermique due au réchauffement du liquide induit une augmentation de la pression de ce liquide.

Ces deux phénomènes contribuent à une augmentation de la pression de fluide dans la caverne. La perméabilité est alors considérée comme un phénomène secondaire à court terme. Pour s'en convaincre, prenons les résultats de mesures in situ obtenus sur une cavité de Hauterives (Drôme) située à une profondeur moyenne de 1600 m. Sur ce site, la température dumassif à la profondeur donnée est de 61 °C ; la température initiale de la saumure stockée dans la caverne est estimée à 26 °C. Un calcul a été réalisé par Bérest et al. (1979) afin d'apprécier l'importance relative des phénomènes mis en jeu. Dans ce calcul, la cavité est assimilée à une sphère ; l'expansion thermique est prise en compte ainsi que la dissolution complémentaire provoquée par l'augmentation de la pression du fluide. La comparaison satisfaisante entre les calculs et les mesures suggèrent que l'expansion thermique est prépondérante sur tous les autres processus, notamment le fluage (figure 9).
<figure><img src="image_57.png" /></figure>

Figure 9. Comparaison entre l'évolution de la pression mesurée au cours du temps et l'évolution donnée par un modèle semi-analytique.

Rayon moyen de la caverne : 47 m; Profondeur moyenne : 1600 m ; température de la roche : 61 °C; température de la saumure initiale : 26 °C.

Si le massif est perméable, nous montrons que l'évolution de la pression dans la caverne est donnée par une solution semi-analytique de la forme (cf. annexe IV) :

$$ \delta P = f_1(t) \otimes \delta P' + f_2(t) \otimes \delta T $$

$\delta P'$ : évolution de la pression de la caverne en l'absence de fuite dans le massif; $\delta T$ : évolution de la température de la caverne ; $\otimes$ : le symbole de convolution de Stieljes tel que :

$$ f \otimes g = \int_{0}^{t} f(t - \tau)dg(\tau) $$

$f_1$ et $f_2$, fonction du temps sont données par :

$$ f_1(t) = \frac{2}{\sqrt{\Delta}} \operatorname{Im} \left[ s_1 e^{s_1^2} \operatorname{erfc}(-s_1 \sqrt{t}) \right] \quad \text{et} \quad f_2(t) = - \frac{1}{\tau_h} \frac{\tau_{hT}}{\sqrt{\tau_T - \tau_h}} \frac{1}{T_o} \frac{2}{\sqrt{\Delta}} \operatorname{Im} \left[ s_1 e^{s_1^2} \operatorname{erfc}(-s_1 \sqrt{t}) \right] $$

Im: l'opérateur partie imaginaire; le complexe $s_1$ donnée par :

$$ s_1 = \frac{-a + i\sqrt{\Delta}}{2} ; \quad \Delta = a^2 - 4b \qquad a = \frac{\sqrt{\tau_h}}{\tau_h'} ; b' = \frac{1}{\tau_h} $$

Les constantes de temps $\tau_h$, $\tau_T$, $\tau_{hT}$ et $\tau_{h'}$ s'écrivent :

$$ \tau_h = \frac{r_o^2}{D_h} ; \quad \tau_T = \frac{r_o^2}{D_T} ; \quad \tau_{hT} = \frac{r_{hT}^2}{D_{hT}} ; \quad \tau_{h'} = \frac{r_o^2}{D_{h'}} $$

Les termes de diffusivité sont :

$$ D_h = \frac{k}{\eta} \frac{1}{M + \frac{b^2}{\lambda_o + 2\mu}} ; \quad D_T = \frac{\lambda_T}{\rho_s C_s} ; \quad D_{hT} = \frac{1}{T_o} \frac{k}{\eta} \frac{1}{-3\alpha_m + \frac{3\alpha_o K_o b}{\lambda_o + 2\mu}} ; \quad D_{h'} = \frac{1}{K^{fl}} + \frac{3}{M} \frac{1}{\mu} $$

k : perméabilité; $\eta$ : viscosité dynamique; M : module de Biot; b : coefficient de Biot; $\lambda_o$ : coefficient de Lamé drainé; $\mu$ : module de cisaillement; $\lambda_T$ : conductivité thermique du sel; $\rho_s$ : masse volumique du sel; Cs : capacité calorifique du sel; To : température initiale; $\alpha_m$ : coefficient de dilatation linéique du fluide relativement au squelette solide; $\alpha_o$ : coefficient de dilatation linéique du squelette drainé; $K_o$ : module de compressibilité drainé; $K^{fl}$ : module de compressibilité de la saumure.
Sur la caverne EZ 53, Durup (1991) a pu mesuré une perméabilité "en grand" de 6. 10-20 m² sur un découvert de Gaz de France (cf. chapitre I, § 2.2). En prenant des valeurs de perméabilité du même ordre ($k \le 10^{-19}$ m²), on peut constater que l'effet de la perméabilité sur l'évolution de la pression de la caverne se fait peu sentir. Ce constat va en faveur de l'hypothèse d'étanchéité des stockages couramment admise par l'ingénieur sur des durées courtes, de quelques années.

## 2.2 Évolution à moyen et à long terme

Dans le cadre d'une étude de sûreté à moyen et à long terme, l'hypothèse d'étanchéité du stockage doit être revue. En effet, il convient d'étudier l'influence d'une perméabilité, même très faible, quelle que soit son origine. D'autre part, dans le cas où existerait une zone endommagée autour de la caverne, engendrée par la décompression des terrains lors du creusement ou du lessivage, le comportement et les conséquences d'une telle zone doivent être analysés, tant du point de vue de la sûreté vis-à-vis des fuites que du comportement mécanique.

### a- Stabilité de la zone endommagée

Plusieurs éléments convergent pour attester la stabilité mécanique de cette zone (si elle existe dans le contexte d'une cavité réalisée par dissolution). D'une part, les observations sur plusieurs mois de cavités abandonnées (par exemple, à Hauterives ou à Vauvert, Bérest, 1989) ne montrent pas de fuites significatives dans le massif qui traduirait la propagation instable d'une zone endommagée. Ces observations peuvent d'ailleurs être étayées par les calculs de structure simples qui suivent.

Nous supposons que la caverne est de forme sphérique, de rayon $r_0$. Son centre se trouve à une profondeur $H_0$, où règne un état de contraintes initiales hydrostatique égal à $\sigma_0$. A court terme, on considère le massif comme poro-élastique. Pour modéliser l'endommagement d'une partie du massif, le formalisme utilisé est celui de la plasticité; un point du massif se trouve dans la zone endommagée si son état de contrainte satisfait un critère de type Tresca, noté $F(\underline{\sigma})$, comme le suggère les résultats expérimentaux obtenus sur le sel en contact avec le fluide (voir figure 7):

$$F(\underline{\sigma}) = \sigma_r - \sigma_{\theta} + 2C$$

$\sigma_r$: contrainte radiale; $\sigma_{\theta}$: contrainte orthoradiale; C: cohésion estimée à 5 MPa pour le critère d'endommagement d'après les données expérimentales.

Dans la zone endommagée considérée comme poreuse, on suppose les pertes de charges négligeables. Cette hypothèse s'appuie sur la dépendance de la perméabilité en fonction de la distance à la paroi de la galerie suivant une loi de puissance (figure 8). Le calcul montre alors que la frontière $x$ de la zone endommagée et le massif sain obéit à l'équation suivante (voir annexe) :

$$\ln \frac{x}{r_0} = -\frac{1}{4C} \left\{ \sigma_0 + P_1 + \frac{4}{3}C - \frac{2}{3(1-\nu_0)} \left[ b(1-2\nu_0)(P_1-P_0) + \alpha_0 E_0 (T-T_0) \right] \right\}$$

$\sigma_0$: contrainte géostatique initiale; $P_1$: pression à l'intérieur de la caverne; $P_0$: pression de pores initiale, $b$: coefficient de Biot; $\nu_0$: coefficient de Poisson drainé; $\alpha_0$: coefficient de dilatation thermique drainée; $E_0$: module d'Young drainé, $T$: température en $x$; $T_0$: température initiale.
<figure><img src="image_58.png" /></figure>

Figure 8. Évolution de la perméabilité et de la pression normalisée (P-P₁)/(P₀-P₁) en fonction du rayon normalisé r/a dans la zone endommagée.

r₀ : rayon de la caverne; P₁ : pression en paroi ; Po : pression de pores initiale; le profil de perméabilité est choisi d'après les données de Stormont et al. (1991) (cf. figure 10, chapitre I) :

$$k(r) = k(r_o) \left(\frac{r}{r_o}\right)^\alpha \quad \text{avec } \alpha \approx -20$$

(La DRZ ne dépasse pas 2r/r₀ conformément aux observations des mêmes auteurs).

On imagine une petite variation de la pression interne dP₁ et une petite variation de la température dT amenant la frontière de la zone endommagée en x₀+dx. Vu la relation précédente, l'incrément dx vérifie alors :

$$\frac{dx}{x_0} = -\frac{1}{4C} \left[ 1 - \frac{2b(1-2v_0)}{3(1-v_0)} \right] dP_1 + \frac{\alpha_0 E_0}{6C(1-v_0)} dT$$

Comme le coefficient de Poisson v₀ et celui de Biot b sont respectivement plus petits que 0,5 et 1, on peut vérifier qu'une augmentation de la pression dans la caverne (dP₁ >0) et une diminution de la température dans la zone endommagée (dT<0) ne provoque pas une propagation de la frontière en x₀ (dx<0).

Or, suite à la fermeture du puits d'accès, on constate sur les sites une augmentation de la pression de la caverne consécutivement à la résorption du déséquilibre thermique (la saumure étant plus froide que la roche, elle a tendance à se dilater en se réchauffant), et au fluage (Bérest, 1985). La pression de la caverne augmentant et la température du massif au voisinage de la caverne baissant faiblement dans un premier temps (la température s'équilibre au bout de plusieurs mois, Bérest, 1990), il est donc raisonnable de penser que la frontière ne se propage pas à court terme.

Ce résultat est confirmé par un calcul sur une géométrie sphérique, en tenant compte du comportement viscoplastique du sel. L'approche étant qualitative et afin de mener les calculs analytiques jusqu'au bout, nous avons pris les hypothèses supplémentaires suivantes :

(a) il n'y a pas de déformation volumique significative dans la zone endommagée ($v_0=0,5$),

(b) l'exposant de la loi de fluage est pris égal à 1 (loi de Maxwell).

(c) le seuil de l'écoulement viscoplastique est négligeable devant le seuil d'endommagement. Cette hypothèse est couramment admise pour le sel qui possède un seuil viscoplastique théoriquement nul (Langer, 1981; Vouille et al., 1981; Waversick et Preece, 1981).

(d) le chargement en paroi est décroissant jusqu'à une valeur tf correspondant à la fin du creusement (ou du lessivage) et donc à la fermeture de la caverne.
(e) la zone viscoplastique n'est pas perturbée par le champs de pression interstitielle et le champs de température.

Le calcul est conduit de la même façon que celui de Tchiyep (1995) qui avait considéré un modèle de viscoplasticité avec rupture, en symétrie cylindrique. La frontière de la zone endommagée s'écrit (annexe III) :

$$ \dot{P}_1 + \frac{E_o}{\eta_s^{-1}} P_1 = -\frac{E_o}{\eta_s^{-1}} \sigma_o - 4C \frac{\dot{x}}{x} - 4C \frac{E_o}{\eta_s^{-1}} Ln \frac{x}{r_o} $$

x : frontière entre le domaine élastique et la zone endommagée ; P₁ : pression dans la caverne; η_s : coefficient de fluidité du potentiel viscoplastique dans la loi de fluage; E_o : module d'Young drainé ; σ_o : contrainte géostatique initiale ; C : cohésion ; r_o : rayon de la caverne.

A l'instant tf, la chute de la pression dans la caverne est terminée, ce qui correspond à la fin du creusement; la caverne est alors isolée. On note avec l'indice o, les valeurs des grandeurs d'état, pression et température, et de la frontière x correspondant à cet instant. Pendant une durée infinitésimale dt, la pression dans la caverne prend la valeur $P_1^o + dP_1$ ($dP_1 > 0$). La frontière x entre la zone endommagée, elle aussi, évolue en $x_o+dx$. Si $dx>0$, la zone endommagée se propage; si $dx<0$, elle s'arrête. Cet incrément $dx$ vérifie :

$$ dP_1 + \frac{E_o}{\eta_s^{-1}} (P_1 + \sigma_o + 4C Ln \frac{x_o}{r_o}) = -4C \frac{dx}{x_o} $$

On prend : $P_1=12$ MPa; $\sigma_o=-22$ MPa; $C=5$MPa.
On vérifie alors que pour $dP_i > 0$, on a bien $dx<0$.

Ainsi, les deux calculs effectués en poro-élasticité et poro-viscoplasticité sont en accord et suggèrent que tant que la pression de la caverne croît, la zone endommagée ne se propage pas. On supposera dans la suite que ce résultat reste vraisemblablement vérifié pour l'évolution à moyen terme (plusieurs dizaines d'années).

## b- Application à une structure souterraine

La caverne est toujours de forme sphérique, de rayon $r_o$. La température du massif est $T_o$ à cette profondeur. Le massif au voisinage du stockage est considéré comme poreux, saturé et faiblement perméable. La pression de pores est supposée égale à $P_o$ dans le massif vierge. On considère en outre que la pression et la température dans la caverne sont respectivement $P_1$ et $T_1$.

Le sel gemme présente un comportement thermo-élasto-viscoplastique caractérisé par le module de Young "drainé" $E_o$, le coefficient de Poisson "drainé" $v_o$ et le fluage stationnaire décrit par la loi de Norton-Hoff :

$$ \dot{\varepsilon}^{vp} = \eta_s(T) (\sigma^{eq})^n \frac{\partial \sigma^{eq}}{\partial \varepsilon} $$

$\eta_s(T)$ : fluidité du sel, fonction de la température (classiquement selon la loi d'Arrhénius : $\eta_s(T)=\eta_0 \exp(-K/T)$)

; $\sigma^{eq}$ : contrainte équivalente de Von Mises $\sigma^{eq} = \sqrt{\frac{3}{2} \underline{S} : \underline{S}}$ ($\underline{S}$ : déviateur de contrainte) ; n : exposant du critère viscoplastique ; K : constante d'activation.
Nous admettons par ailleurs (hypothèse justifiée dans ce contexte) qu'à moyen
terme, tant que l'hypothèse des petites perturbations est vérifiée, le système évolue vers un
état stationnaire caractérisé par :

1°- Une température de la saumure dans la caverne identique à la température du massif, considérée comme uniforme. Un raisonnement sur les temps caractéristiques montre qu'une telle situation est rencontrée au bout d'une trentaine d'années pour une caverne de 300 000 m³ (Bérest, 1990).

2°- Une contrainte déviatorique vérifiant la relation suivante (voir annexe V) :

$$
\sigma_{rr} - \sigma_{\theta\theta} = f(T, n, K, r_0) \frac{\sigma_0 - P_l}{r^{3/n}}
$$

$f(T,n,K,r_0)$ : fonction de la température, du rayon de la cavité, de l'exposant et de la constante d'activation du critère de viscoplasticité :

$$
f(T, n, K, r_0) = \frac{\exp[K/(nT)]}{2} \left[ \int_{r_0}^{\infty} \frac{\exp[K/(nT(r))]}{u^{3/n+1}} du \right]^{-1}
$$

On remarque qu'à l'état stationnaire, le déviateur et plus généralement l'état de
contrainte ne dépendent pas du champ de pression interstitielle dans le massif.
Le calcul complet permet d'obtenir :

$f = \frac{3}{2n} r_o^{3/n}$

3°- Une pression du fluide dans la caverne et dans le massif, éventuellement distinctes.

Dans ces conditions, la pression limite P₁ dans la caverne obéit aux deux équations suivantes :

a- la conservation de la masse de fluide à la paroi ;

b- la diffusion de la pression interstitielle dans le massif.

a.- Conservation de la masse fluide à la paroi :

$$
\begin{align*}
\dot{M} &= \underbrace{V_0 \rho} \quad + \quad \underbrace{\rho_0 \dot{V}} \\
&\stackrel{\text{Variation de la densité du fluide}}{\Rightarrow \text{équation d'état du fluide (a1)}} \\
&= \underbrace{-S_0 q} \quad \stackrel{\text{Fuite dans le massif}}{\Rightarrow \text{loi de conduction du fluide}} \\
&\stackrel{\text{en milieu poreux (a3)}}{\Rightarrow}
\end{align*}
$$

(28)

M : masse de saumure dans la cavité, ρ₀ : masse volumique de la saumure, S₀ et V₀ : respectivement surface et volume initiaux de la caverne, q : débit de fuite à la paroi, V : volume actuel de la caverne.

al - Équation d'état du fluide :

A l'état stationnaire $\dot{T}=\dot{P}=\dot{\sigma}=0$, de même que : $\dot{p}=0$
## a2 - Équation de comportement du massif :

Pour une sphère, nous avons : $\dot{V}=4 \pi r_0^2 \dot{r} = 4 \pi r_0^2 \dot{u}_0$. Avec $\dot{u}_0$: vitesse de convergence de la paroi donnée par l'expression suivante :

$$ \frac{\dot{u}}{r_0} = - \frac{1 - v_o}{E_o} \dot{\sigma}_{\theta\theta} + \frac{v_o}{E_o} \dot{\sigma}_{rr} - \frac{1}{2} \dot{\epsilon}_{vp} $$

Compte tenu de la loi de comportement (expression de $\dot{\epsilon}^{vp}$) et de la forme du déviateur à l'état stationnaire :

$$ \frac{\dot{u}}{r_0} = -\frac{1}{2} \eta_s (T_1) \left(\frac{3}{n}\right)^n \left(\frac{\sigma_o + P_1}{-2\sigma_n}\right)^n $$

La variation de volume de la cavité au cours du temps vérifie alors :

$$ \dot{V} = - \frac{S_o}{2} \eta_s (T_1) \left(\frac{3}{n}\right)^n \left(\frac{\sigma_o + P_1}{-2\sigma_n}\right)^n \quad (29) $$

## a3 - Loi de conduction du fluide :

En considérant que la loi de Darcy est vérifiée pour le sel, nous avons alors en coordonnées sphériques :

$$ q = - \rho_o \frac{k}{\eta^{fl}} \left[ \frac{\partial P}{\partial r} \right]_{r_o} $$

k : perméabilité intrinsèque, $\eta^{fl}$ : viscosité dynamique du fluide.

Pour expliciter l'expression du débit de masse de fluide $q$ à la paroi, il convient d'avoir l'expression de la pression interstitielle $P$ dans le massif.

## b - Équation de diffusion de la pression interstitielle dans le massif

A l'état stationnaire $\dot{T}=\dot{P}=\dot{\sigma}=0$, l'équation de diffusion du fluide interstitielle n'est pas couplée au comportement thermo-mécanique et se simplifie :

$$ \frac{1}{r^2} \frac{d}{dr} \left[ r^2 \frac{dP}{dr} \right] = 0 \Rightarrow \frac{dP}{dr} = \frac{P_1 - P_o}{r^2} r_o $$

D'où

$$ q = \rho_o \frac{k}{\eta^{fl}} \left[ \frac{P_1 - P_o}{r_o} \right] \quad (30a) $$

Ou encore si le massif n'est perméable que dans une couronne de roche (i.e. dans une zone endommagée), d'épaisseur $r_e-r_0$, autour de la caverne :
$$ \frac{1}{r^2} \frac{d}{dr} \left[ r^2 k(r) \frac{dP}{dr} \right] = 0 \Rightarrow \frac{dP}{dr} = \frac{P_1 - P_o}{k(r) r^2} \frac{1}{r_e} \int_{r_o}^{r_e} \frac{dr}{k(r) r^2} $$

D'où

$$ q = \frac{\rho_o P_1 - P_o}{\eta^{fl} r_o^2} \frac{1}{r_e} \int_{r_o}^{r_e} \frac{dr}{k(r) r^2} \qquad (30b) $$

En remplaçant les équations (29) et (30) dans l'équation de conservation de la masse fluide (28), on obtient une équation non linéaire dont l'inconnue $P_1$ est la pression de la cavité :

(si $k$=constante dans le massif)

$$ \frac{1}{2} \eta_s(T_1) \left(\frac{3}{n}\right)^n \left(\frac{\sigma_o + P_1}{-2\sigma_n}\right)^n = \frac{k}{\eta^{fl}} \frac{P_1 - P_o}{r_o} \qquad (31a) $$

(Si $k$=$k(r)$ dans une zone endommagée)

$$ \frac{1}{2} \eta_s(T_1) \left(\frac{3}{n}\right)^n \left(\frac{\sigma_o + P_1}{-2\sigma_n}\right)^n = \frac{1}{\eta^{fl}} \frac{P_1 - P_o}{r_o^2} \frac{1}{r_e} \int_{r_o}^{r_e} \frac{dr}{k(r) r^2} \qquad (31b) $$

Prenons, comme le suggère les résultats *in situ* du WIPP (Stormont *et al.*, 1991) un profil de perméabilité de la forme : $k(r) = k(r_o) \left(\frac{r}{r_o}\right)^\alpha$ (avec $\alpha \approx -20$). L'équation (31b) se réécrit :

$$ \frac{1}{2} \eta_s(T_1) \left(\frac{3}{n}\right)^n \left(\frac{\sigma_o + P_1}{-2\sigma_n}\right)^n = \frac{k_o}{\eta^{fl}} \frac{P_1 - P_o}{r_o} \left[ \frac{1+\alpha}{1-\left(\frac{r_o}{r_e}\right)^{\alpha+1}} \right] \qquad (31c) $$

Remarque :

Les équations (31a) et (31c) permettent de retrouver deux situations limites :

- le massif est complètement étanche : ($k \to 0$) ; la pression dans la caverne tend alors vers l'état de contrainte géostatique initial dans le massif ($P_1 \to -\sigma_o$)

- le massif est poro-élastique et ne présente pas de comportement différé : ($u=0$) ; la pression de la caverne tend alors vers la pression de pores initiale ($P_1 \to P_o$)

Ces équations peuvent être résolues pour différentes valeurs de perméabilité en faisant varier d'une part les facteurs intrinsèques (coefficient de viscosité, exposant de la loi de fluage) et d'autre part les facteurs extrinsèques (rayon et profondeur de la cavité).
## Influence des facteurs intrinsèques

- On retrouve les deux cas limites évoqués : plus la perméabilité est élevée, plus la pression limite $P_1$ tend vers la pression interstitielle initiale du massif $P_0$. De même, plus la perméabilité est faible, plus la pression limite $P_1$ tend vers la pression lithostatique des terrains $\sigma_0$ (Ghoreychi et Cosenza, 1993);

- Dans le cas d'une perméabilité constante du massif (comprise entre $10^{-22} \cdot 10^{-18} \text{m}^2$), la pression limite $P_1$ est d'autant plus élevée que l'exposant n de la loi de fluage stationnaire et la fluidité $\eta_s$ du sel sont faibles.

- Dans le cas d'une présence d'une zone endommagée (équation 31c), la pression limite $P_1$ est d'autant plus proche de la pression de pores $P_0$, que la perméabilité en paroi $k_0$ est grande (le rayon de la zone endommagée $r_e$ étant maintenu constant) et que l'extension de la zone endommagée est importante (perméabilité $k_0$ étant maintenue constante).

## Influence des facteurs extrinsèques

Pour une perméabilité donnée, la pression limite $P_1$ est d'autant plus élevée que le rayon $r_0$ et la profondeur de la cavité sont grands (figure 9).

<figure><img src="image_59.png" /><figcaption>Figure 9 : Évolution de la pression stabilisée dans la caverne à l'état stationnaire pour différentes valeurs de perméabilité. (Profondeur : 1000m, Pression interstitielle initiale : 12MPa, Pression lithostatique : 20MPa, $r_0$ : rayon de la cavité; $1 \text{ Darcy} = 10^{-12} \text{ m}^2$)</figcaption></figure>

## Discussion

Ce calcul simple met bien en évidence que la pression de la saumure développée à moyen terme (quelques dizaines d'années) au sein d'une cavité de stockage abandonnée dépend de trois paramètres essentiels : la nature du sel gemme (aptitude au fluage, perméabilité), la profondeur et la taille de la caverne.
L'influence de cette dernière tient exclusivement à la nature de la condition à la limite en paroi. L'écriture du bilan de masse à la paroi (équation 28) fait intervenir d'une part la variation de volume de la caverne causée par le fluage et d'autre par la fuite dans le massif à travers la surface de la caverne. Par conséquent, le rapport surface/volume de la cavité constitue un paramètre essentiel du calcul.

L'analyse de l'état de contrainte en paroi permet de discuter du risque de rupture éventuelle de la caverne aussi bien en compression qu'en traction :

- le risque de rupture par compression est d'autant plus élevé que la perméabilité est grande. En effet, une augmentation de la perméabilité conduit à une diminution de la pression limite de la caverne et donc à une diminution de la contrainte moyenne et une augmentation du déviateur. Cette situation contribue à se rapprocher du critère de rupture du sel.

Dans la pratique, ce critère n'est pas atteint puisque le déviateur limite à la paroi est très faible :

$$ \sigma_{rr} - \sigma_{\theta\theta} = -\frac{3}{2n}(\sigma_0 + P_1) $$

Le déviateur limite est d'autant plus faible que l'exposant du critère de viscoplasticité n généralement compris entre 3 et 5 est élevé.

- le risque de rupture par traction est en revanche d'autant plus faible que le matériau est perméable. Une diminution de la perméabilité conduisant à réduire le déviateur, la fracturation hydraulique pourrait se produire à très long terme pour un déviateur nul (la résistance en traction à long terme du sel étant pratiquement nulle).

## CONCLUSIONS

Afin d'apprécier qualitativement et quantitativement l'évolution dans le temps du caractère soluble du sel gemme, et à l'aide de modèles théoriques considérant le comportement poro-élastique, l'influence de deux processus a été discuté :

(a) le phénomène du dissolution/cristallisation consécutif à un déplacement de l'état d'équilibre chimique de la solution saline emprisonnée dans les microfissures ;

(b) le phénomène de cicatrisation sous contrainte qui se manifeste macroscopiquement par une baisse différée de la perméabilité d'un échantillon lorsque celui-ci est soumis à un chargement hydrostatique.

Les résultats des modèles montrent que l'ordre de grandeur de la perméabilité n'est pas modifié par les transferts de masse de soluté, et que la cicatrisation est susceptible de restaurer la perméabilité initiale du sel sain. Ces résultats vont dans le sens d'une plus grande sécurité vis-à-vis de migrations éventuelles d'éléments vers la biosphère, dans le contexte d'u stockage souterrain en formations salines.

La prévision de la sûreté à long terme d'un stockage souterrain nécessite également une analyse appropriée de l'évolution de l'étanchéité du milieu géologique au cours du temps. L'influence de sollicitations mécaniques, thermiques et hydrauliques sur l'évolution de la perméabilité du milieu doit être considérée. Ce trait a été bien mis en évidence par le calcul couplé effectué sur une cavité sphérique représentant une cavité de stockage abandonnée remplie de saumure ; l'évolution est de la pression de la caverne à moyen terme (sur une dizaine d'année) est gouvernée non seulement par les propriétés rhéologiques du sel mais également par la perméabilité (naturelle ou induite par endommagement mécanique du champs proche) très faible du massif (< 10⁻¹⁹ m²).
CONCLUSIONS ET PERSPECTIVES
## CONCLUSIONS ET PERSPECTIVES

L'évaluation de la sûreté à long terme de stockages en formations géologiques nécessite la prise en compte simultanée du comportement mécanique et des propriétés de transfert de masse, notamment la perméabilité. Dans le cas du sel gemme, ce couplage est rarement considéré, le sel étant supposé en principe imperméable, comme l'atteste l'existence de poches de fluides (gaz ou saumure) piégées en son sein, sous fortes pressions.

Toutefois, de récents résultats expérimentaux obtenus en laboratoire et in situ remettent en question cette "imperméabilité". Le sel possède une perméabilité mesurable au moins dans deux cas :

a- lorsque le sel est endommagé consécutivement à un chargement mécanique (par exemple, lors de l'excavation d'un ouvrage souterrain) ou thermo-mécanique, la microfissuration induite provoque alors une augmentation de la perméabilité de plusieurs ordres de grandeur.

b- lorsque le matériau présente une teneur en impuretés significative (argiles et anhydrite notamment), ces d'impuretés peuvent constituer des chemins préférentiels d'écoulement des fluides.

Dans le contexte de la sûreté à long terme d'un stockage souterrain, ces résultats soulèvent les questions suivantes :

- dans quelle condition l'endommagement mécanique et les impuretés organisées en un réseau connecté peuvent donner lieu à une perméabilité mesurable pour un sel naturel ?

- quelles sont les conséquences de l'existence d'une très faible perméabilité sur la rhéologie du sel et sur le comportement mécanique à long terme d'un stockage ?

C'est à l'étude de ces questions que ce mémoire s'est consacré. Sa contribution peut se diviser en trois volets :

* *Sur le plan expérimental, nous avons proposé une méthodologie d'étude et mis au point un dispositif expérimental original, pour mesurer dans un contexte difficile, dans une galerie d'une mine de sel, des très faibles perméabilités (inférieures à 10⁻²⁰ m²). Les résultats expérimentaux obtenus à partir d'un essai in situ aux Mines de Potasse d'Alsace (MDPA) ont permis d'affirmer que le sel gemme est perméable au gaz et à la saumure même assez loin des ouvrages souterrains.*

Les faibles débits volumiques de saumure mesurés sont relativement bien reproduits par la loi de Darcy. Les valeurs de la perméabilité intrinsèque et de la pression de pores initiale obtenues sont égales à 2 10⁻²¹ m² et 1 MPa respectivement. Concernant l'essai au gaz, deux points ont été mis en évidence :

- après percolation de saumure dans le massif, la pression capillaire fait diminuer d'un ordre de grandeur le débit massique de gaz.

- la loi de Darcy n'est pas bien adaptée à l'écoulement du gaz dans une formation salifère. L'effet Klingenberg (ou phénomènes de "glissements en paroi" des molécules de gaz) et la présence éventuelle d'une pression capillaire initiale semblent jouer un rôle significatif. En vue d'applications dans le domaine de la sûreté du stockage, la modélisation quantitative (et bien sûr macroscopique) de l'effet Klingenberg constitue une perspective naturelle de ce travail.
Par ailleurs, des essais triaxiaux, sous des pressions de confinement élevées jusqu'à 60 MPa montrent que le sel gemme présente un endommagement plus précoce et une résistance à la compression plus faible, lorsqu'il est directement en contact avec le fluide de confinement. Les critères d'endommagement et de rupture dans le plan "contrainte moyenne-contrainte déviatorique" obtenus sur les échantillons de sel humide (non gainés) se trouvent bien en dessous de ceux obtenus sur un sel sec (échantillons gainés) et ne dépendent pratiquement pas de la contrainte moyenne.

Ces résultats de laboratoire expliquent l'existence d'une perméabilité même pour les échantillons soumis à une forte pression de confinement :

- avant le seuil d'endommagement, le fluide pénètre bien dans le sel à cause de la perméabilité initiale de l'échantillon et ceci malgré le fort chargement hydrostatique appliqué préalablement, durant environ une journée, et les pressions de confinement élevées. La pression du fluide assiste le processus d'initiation de la microfissuration.

- après le seuil d'endommagement, la percolation liée à la microfissuration produite par le chargement déviatorique assiste également, par une augmentation de la pression de pores, le processus de propagation des microfissures.

La dégradation des critères d'endommagement et de rupture observée en condition humide est accentuée lorsque le liquide en contact est de la saumure. Les réactions physico-chimiques dans le système halite-saumure au sein de l'échantillon sont donc reliées au mécanismes d'endommagement et de rupture. Les investigations expérimentales sur ce couplage mécano-chimique mériteraient d'être poursuivies afin de mieux cerner les mécanismes microscopiques mis en jeu.

En outre, des essais de relaxation et de fluage réalisés également sur des échantillons non gainés révèlent que la percolation d'un fluide au cours de la microfissuration (après le seuil d'endommagement) influence non seulement le critère de rupture mais également les propriétés rhéologiques du sel : la viscosité apparente se trouve diminuée.

* *Sur le plan de la modélisation, un cadre théorique macroscopique, a été proposé, considérant le sel gemme comme un matériau poreux au sein duquel une réaction chimique entre le fluide et la matrice a lieu. Malgré l'indétermination qui pèse sur une partie des paramètres introduits et le caractère simplifié de la cinétique chimique, il est désormais possible grâce à ce cadre théorique :*

- d'analyser le couplage entre la déformation macroscopique du matériau et les phénomènes de dissolution-cristallisation produits par la percolation d'un fluide réactif ;

- et donc d'envisager ultérieurement une extension de ce travail qui consiste à prendre en compte les effets cinétiques associés à plusieurs réactions chimiques.

Notons au passage que le cadre théorique proposé, fondé sur les principes thermodynamiques, est assez général : il n'est pas limité au seul cas du sel gemme. Nous nous sommes efforcés de considérer l'ensemble des processus de transport tout en étant conscient que certains (la convection, par exemple) sont rarement rencontrés dans le sel.

Il faut aussi préciser que cette approche, bien que convaincante d'un point de vue théorique, souffre du nombre important de paramètres intervenants dont la détermination reste difficile, délicate et onéreuse au moins pour certains d'entre eux. C'est la raison pour laquelle nous avons procédé, d'une façon simplifiée (et probablement discutable), à une estimation des temps caractéristiques des différents processus. Nous avons ainsi suggéré de découpler, chaque fois que ceci est justifié, un ou plusieurs processus pour simplifier le modèle, alléger le calcul et faciliter l'analyse des résultats.