Modélisation géochimique et numérique des interactions entre des solutions solides et une solution aqueuse

Un logiciel de transport réactif : ARCHIMEDE 

Le logiciel de transport réactif ARCHIMEDE (Algorithme de Résolution de la Chimie Hétérogène IMplicite en 1D) simule les transformations géochimiques dans un système constitué d’une solution aqueuse circulant à travers une roche formée de composés définis, hors d’équilibre. Cette roche occupe un volume élémentaire, auquel est associé une seule maille dans une modélisation unidimensionnelle. Ce code de calcul a été con¸cu à l’Institut Fran¸cais du Pétrole en 1998 par Cassou afin d’améliorer le logiciel NEWKIN (Bildstein [1998]), lui-mˆeme basé sur le modèle KINDIS (Madé [1991], Clément [1992], Fritz [1981]) du Centre de Géochimie de la Surface de Strasbourg : les équations de transport sont exprimées en fonction de la vitesse de Darcy et les équations aux dérivées partielles sont discrétisées de fa¸con totalement implicite. Depuis le dernier rapport de Cassou [1998], Archimède a évolué. Nous allons décrire dans ce chapitre les systèmes d’équations choisis dans la dernière version de ce logiciel afin de modéliser les altérations d’un milieu poreux, constitué de minéraux à composition fixe, parcouru par une solution aqueuse. Etendu à deux dimensions et couplé avec le module de transport polyphasique SARIP ´ ch (Le Thiez [1989] [1992]), qui tient compte aussi bien de la diffusion/dispersion que de l’advection, des variations de la masse volumique du fluide et de la perméabilité, Archimède engendre la troisième version du logiciel de transport réactif DIAPHORE (Cassou [2000]), qui modélise l’écoulement d’un fluide à travers une roche altérée par des réactions chimiques, selon un schéma de discrétisation implicite à deux dimensions (figure 1.1 page suivante). Archimède permet de prendre en compte un gaz dissous à pression partielle constante. Cependant la modélisation n’est pas très élaborée et ne concerne que le gaz carbonique. Comme ce problème est indépendant de celui des solutions solides, on ne l’évoquera pas. Fig. 1.1: Historique des logiciels de transport réactif Archimède et Diaphore. davantage ici. Le lecteur intéressé se reportera à l’annexe A. La résolution des équations décrivant les transformations géochimiques de la roche en déséquilibre avec un fluide requiert le calcul préliminaire de la spéciation initiale en solution. Archimède détermine donc les molalités initiales des espèces aqueuses en cherchant la solution du système de spéciation initiale, constitué des lois d’action de masse des réactions de spéciation et du bilan de masse des éléments. Les réactions au sein de la solution aqueuse étant rapides, celle-ci est supposée à tout instant en équilibre thermodynamique. Archimède évalue les valeurs des fractions volumiques des minéraux et des molalités des espèces aqueuses à chaque pas de temps, à partir du système dynamique, composé des lois de conservation de la masse des minéraux, de celles des éléments et des lois d’action de masse des réactions de spéciation, discrétisé de manière implicite. Ces deux systèmes non linéaires sont résolus par la méthode de Newton-Raphson, qui permet de les linéariser par approximations successives et de calculer une solution approchée.

Equilibre chimique

 Avant de simuler l’évolution du système considéré, le logiciel Archimède doit déterminer la spéciation des espèces aqueuses en solution. 

Spéciation en solution aqueuse 

La spéciation est la répartition des éléments chimiques dissous dans l’eau entre différentes espèces aqueuses, ions simples, ions complexes ou molécules neutres (figure 1.2 page 9). Elle dépend des conditions thermodynamiques, comme la température, et des quantités de chaque élément. Considérons A = (E1,. .. , ENf ) l’ensemble de toutes les espèces chimiques1 présentes dans la solution aqueuse. Soit Ne le nombre d’éléments chimiques présents dans le système. On peut associer à chaque élément El une espèce aqueuse particulière Ek le contenant : Ek ⇋ X Ne l ′=1 αl ′k El ′ avec αl k 6= 0 . Cette espèce aqueuse est totalement déterminée si l’on connaˆıt les Ne coefficients αlk (Saaf [1996])23 : −→Ek = (α1k .. . αNek) t . (1.1) On note B l’ensemble des espèces aqueuses associées à chaque élément chimique El : B = (E1,.. . , ENe ) et VB l’ensemble des vecteurs (−→Ek )16k6Ne associés aux espèces contenues dans B : VB = (−→E1 ,.. . , −→ENe ). 1La liste des symboles utilisés ici est donnée dans le glossaire à la fin de ce mémoire. 2Cette description ne permet pas de traiter les réactions d’oxydo-réduction. En effet, il faudrait alors ajouter aux composantes αlk du vecteur −→Ek la charge zk de l’espèce aqueuse (Saaf [1996]). Nous nous sommes attachés dans ce chapitre à décrire la version du logiciel Archimède qui a été mise à notre disposition, dans laquelle les réactions redox ne sont pas prises en compte. Dans ces conditions, la définition (1.1) des vecteurs −→Ek permet de représenter toutes les réactions simulées par Archimède. 3Notation : soit −→V un vecteur; −→V t représente son transposé. 7 8 Un logiciel de transport réactif : ARCHIMEDE Si les espèces aqueuses de B ont été correctement choisies, alors VB forme une base de R Ne (illustration ci-dessous). B = (Ek) 16k6Ne est appelé l’ensemble des espèces aqueuses de base. Les Nd espèces aqueuses n’appartenant pas à B peuvent donc se former à partir des espèces de base. Elles sont appelées espèces subordonnées. Exemple de spéciation Considérons un système (figure 1.2 page ci-contre) dont l’ensemble des espèces aqueuses est A = (OH−, H2O, H +, H4SiO4, H3SiO− 4 , CO2− 3 , HCO− 3 , H2CO3). (1.2) Il contient alors les quatre éléments chimiques O, H, Si et C. En associant à ces éléments les espèces de l’ensemble4 B = (OH−, H +, H4SiO4, CO2− 3 ), on peut construire l’ensemble VB des quatre vecteurs −→Ek de type (1.1), qui correspondent aux quatre premières colonnes du tableau 1.1. Il est facile de vérifier que VB constitue une base de R4 . On a ainsi défini les quatre espèces de base du système A. Les espèces subordonnées sont alors les espèces aqueuses de l’ensemble A qui n’appartiennent pas à l’ensemble B. Remarquons que ce choix d’espèces de base n’est pas unique. Dans le système A, il aurait été possible d’associer aux quatre éléments O, H, Si et C les espèces de l’ensemble B ′ = (H2O, H3SiO− 4 , HCO− 3 , H2CO3), auxquelles correspondent quatre vecteurs −→Ek ′ de type (1.1), dont les composantes sont données par les quatre dernières colonnes du tableau 1.1. L’ensemble VB′ ainsi formé est aussi une base de R4 . Les espèces aqueuses de l’ensemble B′ auraient donc tout aussi bien pu constituer les espèces de base du système A. En pratique, on évite de choisir comme espèces de base des espèces minoritaires du système, ce qui pourrait conduire à des problèmes numériques. Espèces de base Espèces subordonnées OH− H+ H4SiO4 CO2− 3 H2O H3SiO4 − HCO− 3 H2CO3 O 1 0 4 3 1 4 3 3 H 1 1 4 0 2 3 1 2 Si 0 0 1 0 0 1 0 0 C 0 0 0 1 0 1 1 1 Tab. 1.1: Composantes des vecteurs Ek de type (1.1) pour l’ensemble A des espèces aqueuses défini par (1.2). 4Les espèces aqueuses de B correspondent aux espèces de base associées aux éléments O, H, Si et C dans la base de données d’Archimède (cf annexe F tableau 5.11 page 159). Remarquons que, dans Archimède, les espèces de base définies dans la base de données ne sont pas nécessairement identiques à celles utilisées par le programme lui-mˆeme. L’espèce de base associée à l’oxygène, par exemple, est OH− dans la base de données et H2O dans le programme. OH − H3SiO4− HCO3− H2CO3 H+ H2OH4SiO4 CO3 2− O H Si C Fig. 1.2: Un exemple de spéciation en solution aqueuse — Les éléments chimiques O, H, Si et C sont par exemple associés aux espèces aqueuses de base OH−, H+, H4SiO4 et CO2− 3 respectivement, et aux espèces aqueuses subordonnées H2O, H3SiO− 4 , HCO− 3 et H2CO3. Dans le programme Archimède, l’espèce de base associée à l’hydrogène est H+ et celle associée à l’oxygène est H2O. Les autres espèces de base sont choisies parmi les espèces aqueuses prédominantes5 . La réaction de spéciation associée à l’espèce aqueuse subordonnée Ej (j > Ne) s’écrit : Ej ⇋ X Ne k=1 νjk Ek , (1.3) o`u νjk sont les coefficients stœchiométriques des espèces de base Ek dans la réaction de dissociation de l’espèce aqueuse Ej . Le calcul de la spéciation repose sur la loi d’action de masse et donc sur l’hypothèse que la solution est à l’équilibre. La loi d’action de masse de la réaction de spéciation (1.3) de l’espèce aqueuse subordonnée Ej (Ne + 1 6 j 6 Nf ) relie les activités ak des espèces de base en solution et la constante d’équilibre Kj de la réaction : Y Ne k=1 ak νjk aj = Kj . (1.4) 5 cf remarque sur le changement de base automatique page 13 9 10 Un logiciel de transport réactif : ARCHIMEDE a) Activité des espèces aqueuses chargées Les activités des espèces aqueuses Ei (i 6= 1)6 sont déterminées grˆace aux molalités mi (exprimées en mol/kg(H2O)) et aux coefficients d’activité γi , qui dépendent de la force ionique I de la solution : ai = γi(I) mi , 2 6 i 6 Nf . (1.5) La force ionique I (Lewis et Randall [1961]) d’une solution est définie en fonction des charges zi des espèces aqueuses7 : I = 1 2 X Nf i=1 zi 2 mi . (1.6) Le coefficient d’activité est calculé selon la loi de Debye-H¨uckel étendue (Lietzke et Stoughton [1961]), qui n’est applicable que pour des milieux peu concentrés dont la force ionique reste inférieure à 1 mol/kg : log γi = zi 2 A √ I 1 + a 0 i B √ I + C I , (1.7) o`u A et B sont les constantes de Debye-H¨uckel caractéristiques du solvant et de la température, C représente l’écart à la loi de Debye-H¨uckel et dépend aussi de la température, a 0 i est un paramètre de taille de l’espèce aqueuse Ei . b) Activité de l’eau Pour les solutions diluées, l’activité de l’eau est très voisine de 1 et n’intervient pas dans les lois d’actions de masse : aH2O = 1 . En revanche, lorsque la solution devient plus concentrée, l’activité de l’eau diminue. Pour des solutions moyennement concentrées pour lesquelles I < 1, Helgeson [1969] a proposé un modèle qui permet de calculer l’activité de l’eau en fonction de la force ionique I et du coefficient osmotique φ(NaCl) d’une solution de chlorure de sodium, en considérant que l’activité de l’eau dans une solution complexe est la mˆeme que dans une solution de sel binaire NaCl de mˆeme force ionique : ln aH2O = − 0, 03603 I φ(I). (1.8) 6L’indice 1 correspond à l’espèce de base H2O, dont le calcul de l’activité est explicité dans le paragraphe b) (équation (1.8) de la présente page). 7Remarque : z1 = zH2O = 0 donc l’équation (1.6) peut aussi s’écrire : I = 1 2 X Nf i=2 zi 2 mi . 10 1.1 Equilibre chimique 11 Le coefficient osmotique φ(I) dépend de la force ionique I via l’expression : φ(I) = 1 − 2, 303A b1 3 I · 1 + b1 √ I − 2 ln ³ 1 + b1 √ I ´ − 1 1 + b1 √ I ¸ + b2I 2 + 2b3I 2 3 + 3b4I 3 4 . (1.9) A est la constante de la loi de Debye-H¨uckel (1.7) et b1, b2, b3 et b4 sont les coefficients de Fritz [1975]. Ces constantes dépendent uniquement de la température. 

 Bilan de masse des éléments

 L’équation de bilan de masse d’un élément El (1 6 l 6 Ne) s’écrit : cl = X Nf i=1 αli mi , 1 6 l 6 Ne , (1.10) avec – cl (mol/kg(H2O)) la concentration totale de l’élément El en solution, – αli le nombre de moles de l’élément El contenues dans une mole de l’espèce aqueuse Ei , – mi (mol/kg(H2O)) la molalité de l’espèce aqueuse Ei . 

 Système de spéciation initiale 

Le calcul de la spéciation initiale8 de la solution aqueuse en système hétérogène requiert la résolution des Nd lois d’action de masse (1.4) des réactions de spéciation des espèces aqueuses subordonnées. Les variables sont les Nf molalités des espèces aqueuses. Etant donné que le système constitué des lois d’action de masse des réactions de spéciation ne contient que Nd relations, et qu’un système numérique ne peut ˆetre résolu qu’à la condition que le nombre d’équations soit égal au nombre d’inconnues, il faut Nf − Nd = Ne relations supplémentaires à Archimède pour déterminer une solution au problème de spéciation. Les mesures des quantités totales d’éléments chimiques en solution étant imprécises ou incomplètes, on préfère en contraindre certaines, via l’équation de bilan de masse (1.10), et équilibrer la charge de la solution à l’aide d’une espèce inactive du système (ianion), via 8Le programme commence par résoudre le système de spéciation initiale lié à l’hypothèse de l’équilibre homogène au sein de la solution aqueuse initiale. En effet, comme les molalités des espèces aqueuses définies par l’utilisateur ne concernent que les espèces de base, le programme doit déterminer les concentrations des espèces subordonnées dans l’eau initiale. Par la suite, l’évolution au cours du temps de la spéciation dans la solution aqueuse est simulée grˆace à la résolution du système dynamique décrit au § 1.2.4 page 29. 11 12 Un logiciel de transport réactif : ARCHIMEDE l’équation d’électroneutralité de la solution aqueuse : X Nf i=1 zi mi = 0 . (1.11) L’espèce inactive Eianion , dont l’ajout permet d’imposer l’électroneutralité, est spécifiée dans le fichier «base.dat» de la base de données (cf annexe F page 157). Dans nos simulations, il s’agit de l’espèce aqueuse anionique Cl−, qui correspond à l’élément Cl. En effet, le chlore n’intervient pas dans les équilibres hétérogènes si le système ne contient pas de sel comme phase solide. L’oxygène est associé à l’espèce H2O, dont la molalité, exprimée en moles par kilogramme d’eau, est connue : mH2O ≃ a55 = 55, 508 . (1.12) Cependant, afin de faciliter la résolution des systèmes d’équations, la molalité de l’eau mH2O est considérée comme une inconnue, au mˆeme titre que les molalités des autres espèces aqueuses, et l’équation (1.12) comme une relation supplémentaire à vérifier. Le pH de la solution aqueuse présente initialement est fourni par l’utilisateur. Sa valeur permet de déterminer l’activité de l’espèce aqueuse H+, et donc sa molalité : pH = − log aH+ = − log (γH+ mH+ ). (1.13) Les (Ne − 3) relations supplémentaires sont données par les équations de bilan de masse (1.10) des éléments en solution. Dans Archimède, on choisit de fixer les concentrations totales des éléments autres que l’oxygène (l 6= 1), l’hydrogène (l 6= 2) et l’élément associé à l’espèce aqueuse inactive ianion fixée par l’utilisateur (l 6= ianion). Le système qui permet de déterminer la spéciation de la solution aqueuse se compose donc des équations suivantes : – la relation (1.12) associée à l’oxygène fixant la molalité de l’eau, – la relation (1.13) associée à l’hydrogène, – les (Ne − 3) équations de bilan de masse (1.10) fixant les concentrations des autres éléments, – la relation d’électroneutralité (1.11) associée à l’élément ianion, – les Nd lois d’action de masse (1.4) des réactions de spéciation des espèces aqueuses subordonnées.