Le logiciel géochimique PhreeqC

Présentation du logiciel PhreeqC

PhreeqC est un logiciel de calcul géochimique de systèmes triphasique (liquide, solide, gaz) développé par l'USGS. C'est un logiciel libre. Les calculs chimiques peuvent être couplés avec des équations de transfert de matière par advection et/ou diffusion. La modélisation du transport s'effectue sous une dimension. Par ailleurs, le logiciel permet également de prendre en compte la cinétique des réactions. Les systèmes d'études sont fermés ou constitués d'un ensemble de systèmes fermés.

La simulation d'un système nécessite deux types de renseignement. Tout d'abord, une base de données référençant les données thermodynamiques des espèces chimiques. Plusieurs bases de données sont utilisables pour le logiciel. Le choix dépend, de la complexité de la simulation à effectuer mais aussi de l'objectif de la simulation. En exemple, la banque de données "iso.dat" renseigne des informations sur les isotopes les plus courants et sera utilisée en cas d'étude isotopique d'un système ou d'une réaction.

Le second renseignement consiste en les instructions qui permettent de réaliser les réactions recherchées. Le logiciel regroupent pas moins d'une vingtaine d'instructions permettant diverses études géochimiques. 

     Capacités du logiciel PhreeqC

A partir des différentes instructions, le logiciel permet de modéliser un grand nombre de réactions possibles:

  • Les réactions d'oxydo-réduction en phase aqueuse,
  • La précipitation ou dissolution de minéraux en solution,
  • La modélisation des intéractions entre plusieurs phases et notamment les réactions de sorption ou complexation sur les surfaces solides.
  • La simulation temporelle des réactions en définissant les constantes cinétiques de celles-ci.

De plus le système offre une certaine liberté pour définir un système:

  • Tout d'abord les bases de données peuvent être modifiées et complétées lorsqu'un composant n'est pas disponible. Nous utiliserons cette possibilité pour introduire certains composés des résidus de bauxite.
  • Les paramètres définissant le système  peuvent être choisi entre le pe, le pH et la concentration des éléments en solution. Le logiciel calcul alors les variables restantes pour atteindre l'équilibre et l'électroneutralité.
  • Le potentiel électrique pe, peut être défini via un couple oxydant réducteur.

    Limites du logiciel PhreeqC

Comme tout logiciel de modélisation PhreeqC possède certaines limites dans la modélisation des réactions chimiques. Le travail effectué par la suite est concerné principalement par les deux limites suivantes :

  • La modélisation des solutions aqueuses est définie implicitement par les informations thermodynamiques entrées dans la bases de données. Le programme utilise un modèle plus ou moins précis, suivant les données bibliographiques disponibles. Les résultats des simulations dépendent alors fortement de la validité et l'exhaustivité des bases de données. Dans notre étude, nous serons limités par l'exhaustivité des bases de données mises à notre disposition. La base de données "llnl.dat" retenue pour nos travaux est la plus complète. Cependant, elle est aussi la moins documentée et testée.
  • Le modèle de transport utilise la méthode des éléments finis avec un schéma explicite qui peut être à l'origine de dispersion numérique des résultats.

Présentation des modèles utilisés dans nos travaux

Notre étude se consacre à l'identification des dissolutions et précipitation ayant lieu au cours du mélange. L'étude des phénomènes de sorption et complexation en surface des minéraux nécessite une calibration expérimentale du logiciel. Par défaut de données, nous ne pouvons pas traiter ces réactions. De même la cinétique des réactions ne sera pas étudiée.

     Modélisation des solutions aqueuses et des intéractions triphasiques

Modélisation des solutions aqueuses

L'instruction "solution" permet de définir une solution aqueuse. Il est possible de choisir la température, le pH, le pe, et les éléments chimiques à introduire. Le logiciel chimique calcule la spéciation chimique des éléments selon le modèle d'équilibre thermodynamique, représenté à l'aide des lois d'action de masse (\ref{lois_masse}) :

\begin{equation}
\label{lois_masse}
K_{i} = a_{i}\prod_{m}^{M} a_{m}^{-C_{m,i}}
\end{equation}

  • $a_{i}$ et $a_{m}$ : activité des l'espèces;
  • $c_{m,i}$ : coefficient stoechiométrique de l'espèce m dans la réaction ramené pour une molécule de l'espèce i
  • $M$ : nombre d'espèce mises en jeu dans la réaction
  • $K_{i}$ :Constante d'équilibre thermodynamique de la réaction

Les activités des espèces en solution sont reliés aux concentrations par le modèle des solutions aqueuses (\ref{equa2}):

\begin{equation}
\label{equa2}
a_{i}=\gamma_{i}C_{i}
\end{equation}

  • $ \gamma_{i}$ : coefficient d'activité de l'espèce i,
  • $C_{i}$ : concentration molaire de l'espèce i dans la solution aqueuse

Le coefficient d'activité de l'espèce i est déterminé suivant un modèle de solution aqueuse. Trois modèles sont ici disponible pour définir la solution aqueuse: le modèle de Davies (\ref{Davies}), utilisé par défaut car il ne nécessite pas d'entrée de données spécifiques à l'espèce; un modèle étendu de Debye Hückel ("Wataq model, \ref{Wateq}), ainsi que le modèle de Pitzer (\ref{Debye}) :

\begin{eqnarray}
log\gamma_{i} &= &-A\cdot z_{i}^2( \frac{\sqrt{I}}{1+\sqrt{I}}- 0.3 I) \qquad I < 0,5 \text{ mol.L$^{-1}$} \label{Davies} \\
log\gamma_{i} &= &-A\cdot z_{i}^2 \frac{\sqrt{I}}{1+rB\sqrt{I}}-\dot{B} I \qquad I < 1\text{ mol.L$^{-1}$}   \label{Wateq} \\
log\gamma_{i} &= &-A\cdot z_{i}^2 \frac{\sqrt{I}}{1+\sqrt{I}} + \sum_{c}\sum_{a}B_{c,a}m_{c}m{a} \qquad I < 1\text{ mol.L$^{-1}$}   \label{Debye}
\end{eqnarray}

  • $z_{i}$ : charge de l'espèce i considéré
  • $I$ : force ionique
  • A, B sont des constantes dépendantes de la température
  • $B_{c,a}$ constante spécifique à la solution
  • $m_{a}$ , $m_{c}$ masse molaire des anions, cations en solution
  • r est un paramètre de dimension représentant l'ion en solution

Dans le cadre de notre étude, il sera important de vérifier la concordance du modèle utilisée avec la force ionique élevée de l'eau de mer étudiée. La figure 1, ci-dessous résume les domaines de validités des différents modèles de solutions aqueuses. Le modèle de Pitzer n'est pas représenté ci dessous. La modélisation établie par ce modèle est globalement supérieure à celles dérivées de la théorie de Debye Huckel, cependant il est plus difficile à appliquer à cause du nombre de paramètres à définir.

Domaine de validité des modèles de solutions aqueuses
Figure 1 : Correspondance entre a force ionique de la solution et les coefficients d'activités des espèces (sources : d'après Garrels et Christ, 1965 et modilié par Merkel et Planer-Friedrich, 2008)

Outre la loi d'action de masse, le logiciel PhreeqC cherchera à satisfaire deux autres bilans chimiques : l'électroneutralité de la solution aqueuse, le bilan de matière pour chaque élément entré dans la solution aqueuse. Dans l'absolu, le système devrait être bouclé pour l'ensemble des éléments du système. Néanmoins, le logiciel offrent la possibilité d'utiliser certains éléments pour satisfaire le bilan en masse du système. Par exemple, le pH ou peut être utiliser pour satisfaire l'électroneutralité.

Intéraction liquide - gaz

Pour prendre en compte les équilibres entre le gaz dissous dans les solutions aqueuses et le gaz présent à l'interface du liquide (par exemple équilibre O2 dissous et pression partiel de O2 dans l'air), la loi de Henry (\ref{Henry}) est utilisée. Elle relie la concentration de l'élément dissous en phase aqueuse à sa pression partielle dans la phase gaz :

\begin{equation}
\label{Henry}
P_{i} = C_{H} a_{i}
\end{equation}

  • $C_{H}$ = Constante de Henry
  • $P_{i}$ = Pression partielle de l'espèce i

Intéraction liquide - solide

PhreeqC permet de représenter deux types de réactions de dissolution / précipitation : les solides sont modélisés soit sous une forme de phases pure, soit sous la forme de solution solide. Les réactions sont alors modélisées par des lois d'action de masse telles que l'équation (\ref{lois_masse}). La modélisation diffère dans le calcul de l'activité du solide :

  • Pour les phases pures, l'activité est considérée unitaire;
  • Pour les solutions solides, l'activité du minéral est définie par sa fraction molaire dans la solution solide.

Lors de la réaction d'un solide, le logiciel dissout intégralement le minéral dans la solution et cherche ensuite à établir l'équilibre d'action de masse. A l'état final, il est possible de calculer l'indice de saturation (SI, \ref{SI}) à l'aide du produit d'activité ionique (PAI, \ref{PAI}) :

\begin{eqnarray}
SI &=& log(\frac{PAI}{K})  \label{SI}  \\
PAI &=& a_{i}\prod_{m}^{M} a_{m}^{-C_{m,i}}  \label{PAI}
\end{eqnarray}

Avec K, la constante d'équilibre de la réaction de dissolution. L'indice de saturation permet de déterminer si la solution est sous saturée ($ SI < 0 $); à l'équilibre ($SI = 0$) ou sursaturée ($SI > 0$). L'équilibre n'est pas obligatoirement atteint.

     Résumé

Le schéma suivant résume la modélisation de la solution aqueuse et des intéractions tri-phasiques :

Schéma de fonctionnement du logicie PhreeqC

Modélisation du transport via PhreeqC

 Le modèle de transport de PhreeqC permet d'intégrer la diffusion et l'advection d'une solution aux calculs de réactions chimiques. Ces transferts sont modélisés à partir de l'équation générale une dimension d'advection diffusion avec terme source/puits :

\begin{equation}
\label{advec_diffu}
\frac{∂C_{i}}{∂t} = -\nu·\triangle C_{i} + D_{L}·\nabla^{2}C_{i} +\sum R_{i}  
\end{equation}

  • $\nu $ : vitesse d'advection de la solution
  • $D_{L} $ : coefficient de diffusion longitudinale
  • $R_{i} $ : terme source/puits représentant les réactions chimiques

        Le modèle impose un coefficient de diffusion ainsi qu'une vitesse d'advection commune à toutes les espèces aqueuses. Concernant les bornes du milieu, trois types de conditions peuvent être imposés : Conditions aux limites de Cauchy, Dirichlet ou Neumann.

Méthode de calcul employée par le logiciel PhreeqC

      L'équation (\ref{advec_diffu}) est résolue par la méthode des différences finies avec un schéma centré dans l'espace et avancée temporellement. Les intéractions chimiques sont calculées séparément des termes de transport à chaque étape. Concrètement le logiciel calcule le transport advectif, puis l'équilibre chimique de la solution et enfin la diffusion. L'équilibre chimique est de nouveau calculé ensuite.

méthode de calcul avec transport

.