Brevet FR2894672A1 — publié le 20070615

La présente invention concerne une méthode pour déterminer les processus hydrodynamiques conduisant au transport de constituants chimiques, ainsi que les processus chimiques intervenant dans un système géochimique. En particulier, la présente invention permet de déterminer les capacités de stockage de gaz acides (H2S, CO2, CH4) dans un milieu souterrain, par modélisation 15 des processus de transport et de réactions chimiques en milieu poreux. Toutes les échelles utiles peuvent être concernées : laboratoire, puits, réservoir et étude régionale de bassin. La méthode peut également trouver des applications directes dans les études de dépollution des sols par des procédés multiphasiques ("venting", utilisation de 20 solvants chlorés). Etat de la technique La réduction des émissions de gaz acides tels que le CO2 industriel est devenu un enjeu majeur pour l'environnement, et donc pour de nombreux industriels. Une 25 des solutions pour réduire ces émissions repose sur la séquestration géologique de ces gaz. Le principe consiste à capter le CO2 à son point d'émission (centrale électrique, cimenteries, aciérie ...), à le concentrer et à le transporter vers un site géologique adéquat pour son stockage. Pour cela, des technologies pour capturer et stocker du CO2 sont d'ores et déjà disponibles. Cette filière est largement étudiée depuis quelques années au niveau international, surtout par les industries pétrolières qui y voient un second avantage : la valorisation du stockage géologique du dioxyde de carbone en injectant le CO2 en appui à une opération de récupération assistée de pétrole (EOR ; Enhanced Oil Recovery) ou la récupération assistée du gaz naturel dans les veines de charbon non exploitables (ECBM ; Enhanced Coalbed Methane) et dans les gisements de gaz naturel (EGR ; Enhanced Gas Recovery). Dans ce cadre, l'option envisagée pour le stockage du CO2 est l'utilisation de réservoirs souterrains naturels (réservoirs épuisés de pétrole et de gaz, gisements houillers non exploitables, aquifères salins et cavités salines). Ce type de réalisation est illustrée sur la figure 1, où un puits producteur (PROD) extrait des hydrocarbures du sous-sol, tandis qu'un puits injecteur (INJ) réinjecte du CO2 dans le réservoir. Le coût associé à la mise en place d'un site de stockage étant important, les enjeux économiques et les impacts environnementaux étant également très importants, il est nécessaire, pour mettre en place de telles méthodes de stockage de gaz acides dans des réservoirs géologiques naturels, de parfaitement déterminer les conditions de stockage sur le long terme. Il s'agit principalement d'étudier la migration du gaz au sein du réservoir et les réactions de précipitation et de dissolution pour déterminer la qualité d'un site de stockage au cours du temps Une solution technique à ce problème d'évaluation d'un site de stockage à long terme est apportée par les techniques permettant de modéliser des systèmes géochimiques. Les documents suivants que l'on va citer dans le cours de la description ci-après, illustrent l'état de la technique : [1] G.T. Yeh and V.S. Tripathi "A Critical Evaluation of Recent Developments in Hydrogeochemical Transport Models of Reactive Multichemical Components", Water Resources Research, Vol. 25. No. 1,93-108, January 1989. [2] Y. Le Gallo, O. Bildstein and E. Brosse "Coupled reaction-flow modeling of diagnetic changes in reservoir permeability, porosity and mineral compostion.", 30 Journal of Hydrology 209 (1-4), 366-388, Agust 1998. [3] J.W. Johnson, J.J. Nitao, C.I. Steefel, K.G. Knauss "Reactive transport modelling of CO2 storage in saline aquifers to elucidate fundamental processes, trapping mechanisms, and sequestration partitioning", GS-London Spec. Pub. ms, August 8, 2002. [4] Long Nghiem, Peter Sammon, Jim Grabenstetter and Hiroshi Ohkuma, "Modeling CO2 Storage in Aquifers with a Fully-Coupled Geochemical EOS Compositional Simulator", SPE 89474, April 004. [5] Tianfu Xu and K. Pruess "Modeling multiphase no-isothermal fluid flow and reactive geochemical transport in variably saturated fractured rocks: 1. 10 Methodology", American Journal of Science, Vol. 301, P. 16-33, January 2001. [6] Hui Cao "Development of techniques for general purpose simulators", PhD Stanford University, June 2002. [7] Jacob Bear "Dynamics of Fluids in Porous Media", American Elsevier, 1972. Reprinted with corrections, Dover, New York, 1988 . 15 Les trente dernières années ont vu une forte augmentation du nombre de méthodes permettant de modéliser des systèmes géochimiques, qu'ils soient naturels ou provoqués. On appelle système géochimique un milieu poreux, tel qu'un volume de roche, dans lequel circule au moins une phase aqueuse et éventuellement une 20 phase gazeuse. Ces systèmes sont caractérisés par deux types de processus : - un processus de transport de fluides ; - un processus de réactions chimiques entre les différents éléments que sont le milieu et les phases aqueuse et gazeuse : roche/eau et eau/gaz par exemple. Une modélisation précise de ces systèmes géochimiques permet par exemple 25 de déterminer les mécanismes de piégeage des gaz acides dans des formations géologiques. Une telle modélisation permet également de déterminer les conditions de dépollution de sols par des procédés multiphasiques. Les modèles utilisés doivent être représentatifs du comportement de ces systèmes géochimiques. De tels modèles sont basés sur le principe de conservation de la masse la masse n'est ni créée ni détruite dans le système, mais transférée entre le solide, la phase aqueuse et la phase gazeuse. Il existe deux types de modèle. - Les modèles qui ne prennent pas en compte les processus de transport. Ces modèles sont appelés "modèles de réactions chimiques" ou simplement "modèles batch". - Les modèles qui considèrent à la fois les processus de transport et de réactions chimiques sont appelés "modèle de transport réactif'. On parle alors de modélisation couplée. La définition de tels modèles est complexe et leurs utilisation peut être limitée par les capacités des ordinateurs et/ou le temps de calcul qui croît considérablement en fonction du nombre de constituants à considérer. Les constituants chimiques, comme par exemple CO2, H2O ou NaCl, peuvent être présents soit dans la phase gaz, soit dans la phase eau. Dans un premier temps, les modèles de systèmes géochimiques ont été utilisés pour comprendre les phénomènes de base liés à la chimie de l'eau ou à l'évolution de processus diagénétiques (cf. [2]). Depuis près de dix ans, les modèles se sont tournés vers l'analyse d'impact environnemental de la pollution des sols, de l'émission de biogaz dans les décharges ou l'étude du stockage de déchets radioactifs (cf. [3] [4]). Enfin depuis peu, leur utilisation s'est étendue à l'étude du stockage géologique de gaz acides (H2S, CO2, CH4). Un modèle de transport réactif permet de simuler les processus hydrodynamiques conduisant à la migration et à la dispersion de constituants (le transport) associées aux processus géochimiques (les réactions chimiques). Initialement utilisés pour résoudre des systèmes hydrogéologiques (cf. [2]), les modèles de transport réactif sont de plus en plus utilisés pour simuler l'injection, la migration et le stockage de CO2 dans des couches géologiques (cf. [3] [4]). L'introduction d'un constituant initialement présent sous forme gazeuse (le CO2) oblige à considérer des modèles non plus monophasique (uniquement l'eau), mais diphasiques eau - gaz et donc à augmenter le nombre d'inconnues du modèle. D'un point de vue général, un modèle de transport réactif est une méthode permettant de trouver la solution d'un système d'équations du type : a(0.0 =L(ci)+R(C;) at où Pest la porosité du milieu, ci la concentration du constituant i dans les phases mobiles (fluides) et Ci la concentration totale (fluides + roche). On a ainsi la relation : CL = ci + i, où c l est la concentration immobile (concentration dans la phase immobile : la roche). L(ci) définit l'opérateur de transport du constituant i et R(CZ) l'opérateur de réaction chimique. Le système (1) est difficile à résoudre en raison du nombre important d'inconnues et des fortes non linéarités induites par les termes de réaction. On recense aujourd'hui essentiellement deux approches pour résoudre de tels systèmes : 1. une méthode de résolution séquentielle itérative ; 2. une méthode de résolution simultanée. Chacune de ces approches peut être résolue de façon explicite ou implicite. On parle de résolution explicite quand on cherche la solution d'un système temporel à t + dt à l'aide de l'état t. 15 La première méthode, la méthode de résolution séquentielle itérative, est la plus couramment utilisée pour résoudre ce type de problème (cf. [2]). Elle consiste à résoudre le système d'équations en deux étapes : une étape de transport sans réaction et une étape de réactions chimiques maille par maille (ceci au court d'un même pas de temps). Yeh et Tripati ont démontré (cf. [1]) qu'il faut ajouter une boucle itérative 20 pour être en mesure de minimiser suffisamment les erreurs numériques sur la conservation de la masse. Cette méthode de résolution séquentielle a été développée et analysée pour des modélisations d'écoulement réactif monophasique eau, 2D ou 3D, en milieu poreux. Elle a été étendue récemment à la modélisation de problèmes multiphasiques par Xu et al. (cf [5]). 25 La seconde méthode, la méthode de résolution simultanée, résout le système complet en une seule étape. Cette méthode permet de résoudre des systèmes géochimiques complexes et elle est théoriquement très robuste. En pratique, comme l'ont montré Yeh et Tripathi dans un papier de référence (cf. [1]), elle peut devenir rapidement coûteuse en temps CPU et en place mémoire sur des systèmes 30 complexes. Les modèles élaborés à partir de ces approches étaient jusqu'à présent (1) 6 2894672 réservés à des écoulements monophasiques eau. Des études récentes comme celle de Nghiem et al. (cf. [4]) montrent qu'on peut maintenant envisager une approche de ce type pour des problèmes multi-phasiques compositionnels en 3D mais avec un nombre restreint de constituants, ce qui n'est pas suffisamment précis pour permettre 5 une application dans le cadre de l'évaluation des capacités de stockage de gaz acides. La méthode selon l'invention tire parti des avantages des méthodes séquentielles (pas de temps locaux, parallélisation, taille mémoire réduite) tout en évitant d'ajouter une boucle itérative supplémentaire. Elle répond ainsi aux exigences 10 de l'évaluation des capacités de stockage de gaz acides (H2S, CO2, CH4) d'un milieu souterrain, ce dernier représentant un système géochimique multiphasique complexe. La méthode selon l'invention L'invention concerne une méthode pour déterminer les capacités d'un milieu géologique souterrain à stocker au moins un gaz acide, le milieu comprenant une phase solide et une phase aqueuse. La méthode comporte les étapes suivantes : a) on découpe ledit milieu en un ensemble de sous volumes appelés mailles ; b) on définit un modèle de transport réactif multiphasique, en couplant un sous modèle de type réservoir multiphasique compositionnel , dédié à la détermination des taux de transfert liés au gaz, avec un sous modèle de type transport réactif monophasique , dédié à la détermination des taux de transfert liés à ladite phase aqueuse et à ladite phase solide ; c) on réalise une acquisition de données pour initialiser ledit modèle de transport réactif multiphasique ; d) on détermine, pour une période de temps t donnée, des processus de transport de constituants chimiques et des processus de réactions chimiques au sein dudit milieu géologique, en résolvant par couplage ledit modèle de transport réactif en toute maille dudit milieu, pour obtenir au moins : - les quantités de gaz présent dans ledit milieu géologique, ayant migré hors dudit milieu et en surface, ayant précipité, et dissout dans la phase aqueuse ; - les altérations de la phase solide par précipitation et dissolution. e) on en déduit les capacités dudit milieu géologique à stocker ledit gaz acide. Selon l'invention, on peut résoudre ledit modèle de transport réactif sur la période de temps t, à l'aide d'un schéma d'évolution en temps, dans lequel on 10 découpe cette période de temps en un nombre fini de pas de temps et en réalisant les étapes suivantes : a) on définit la durée d'un pas de temps pour synchroniser les deux sous modèles ; b) on résout le sous modèle réservoir multiphasique compositionnel ; 15 c) on résout le sous modèle ((transport réactif monophasique à l'aide des résultats obtenus à l'étape b) ; d) on modifie la durée dudit pas de temps et l'on réitère à partir de l'étape b) en tenant compte des résultats obtenus à l'étape c), jusqu'à couvrir ladite période de temps t ; 20 Les taux de transfert liés au gaz peuvent comporter au moins l'un des termes suivants : termes sources dudit gaz ; - termes de changement de phase entre le gaz et la phase aqueuse ; termes de flux de gaz au sein dudit milieu. 25 Les taux de transfert liés à ladite phase aqueuse et à ladite phase solide peuvent comporter au moins l'un des termes suivants : termes de réactions de précipitation et de dissolution entre la phase aqueuse et la phase solide ; - termes de réaction de spéciation au sein de la phase aqueuse ; termes de flux de la phase aqueuse au sein dudit milieu ; - termes source de la phase aqueuse. Selon l'invention, on peut déterminer à partir dudit modèle de transport réactif multiphasique au moins l'un des phénomènes suivants : - les déplacements de ladite phase aqueuse et dudit gaz ; les échanges de constituants entre la phase aqueuse et le gaz ; les échanges de constituants entre la phase aqueuse et la phase solide ; - les échanges de constituants entre constituants au sein de la phase aqueuse ; les échanges de constituants avec l'extérieur dudit milieu. Selon un mode de réalisation, on peut résoudre le modèle de type transport réactif monophasique de manière implicite à partir d'un sous modèle de type réacteur géochimique OD ouvert en réalisant les étapes suivantes : - on définit des possibilités d'échanges de constituants entre lesdites mailles, définissant ainsi un graphe de liaison ; - on détermine un ordre de parcours desdites mailles en analysant ledit graphe de liaison, de façon à parcourir les mailles dans l'ordre des potentiels décroissants; - on parcourt chaque maille selon ledit ordre de parcours et on résout ledit modèle de type réacteur géochimique pour chaque maille, à l'aide d'une 20 technique de ligne de courant. Si l'on ne peut déterminer l'ordre des potentiels décroissants, on peut résoudre ledit modèle de type transport réactif monophasique de manière semi-implicite. D'autres caractéristiques et avantages de la méthode selon l'invention, 25 apparaîtront à la lecture de la description ci-après d'exemples non limitatifs de réalisations, en se référant aux figures annexées et décrites ci-après. Présentation succincte des figures - La figure 1 montre une représentation schématique d'une opération de stockage géologique de gaz acide. - La figure 2 décrit de manière schématique le couplage entre un modèle de 5 réservoir compositionnel multiphasique (RESERVOIR) et un modèle de transport réactif monophasique (TREACT). - La figure 3 décrit l'implémentation particulière d'un modèle de transport réactif monophasique (TREACT), à partir d'un modèle de type géochimique (GEOCH) appliqué sur des lignes de courant. 10 - Les figures 4 et 5 illustrent un cas favorable et un cas défavorable pour appliquer la méthode des lignes de courant. Description détaillée de la méthode La méthode selon l'invention permet de déterminer les processus 15 hydrodynamiques conduisant au transport des espèces chimiques ainsi que les processus chimiques intervenant dans un système géochimique. La méthode, non itérative, respecte le principe de conservation de la masse de tous les éléments chimiques présents dans le système. La méthode tire parti des avantages des méthodes séquentielles (pas de temps 20 locaux, parallélisation, taille mémoire réduite) en évitant d'ajouter une boucle itérative supplémentaire. Cette méthode s'appuie sur le couplage de deux sous modèles : un modèle de type transport multiphasique compositionnel non réactif (nommé RESERVOIR) et un modèle de transport réactif monophasique dans un champ de vitesses donné (nommé TREACT). 25 La méthode est décrite ci-après dans le cadre de l'évaluation des capacités de stockage de gaz acides d'un milieu poreux tel qu'un réservoir géologique. Selon ce mode de réalisation, la méthode comporte les cinq étapes suivantes : 1) on découpe le milieu poreux en un ensemble de sous volumes appelés 30 mailles et on définit des liaisons entre les mailles ; 9 2) on établit un système d'équations définissant un modèle de transport réactif multiphasique ; 3) on réalise une acquisition de données pour alimenter le modèle ; 4) on évalue, à partir du modèle de transport réactif multiphasique, les 5 processus hydrodynamiques conduisant au transport des espèces chimiques ainsi que les processus chimiques intervenant dans un système géochimique ; 5) on évalue, à partir du modèle de transport réactif multiphasique, les capacités de stockage de gaz acides du milieu poreux. 10 1) Maillage du milieu poreux et définition des liaisons Selon le principe des méthodes numériques aux volumes finis, le domaine physique étudié (le réservoir souterrain) est simplement découpé en volumes élémentaires de manière à ce qu'ils forment une partition du domaine. En pratique, les mailles sont des polygones. C'est ce qu'on appelle un maillage du domaine. 15 Chaque maille possède ses variables (pression, concentrations, porosité, saturation,...) qui représentent des valeurs moyennes des variables physiques que l'on souhaite observer. Les mailles doivent ensuite être liées entre elles en ajoutant des liaisons entre les mailles. Ces liaisons correspondent aux chemins possibles de la matière dans le 20 réservoir maillé. Les liaisons définissent donc les possibilités d'échanges de matière entre deux mailles : il y a une liaison entre deux mailles si ces deux mailles peuvent échanger de la matière. Ces transferts peuvent se faire d'une maille vers une autre uniquement (liaison orientée), ou entre deux mailles dans un sens ou dans l'autre ; on dit alors que la liaison n'est pas orientée intrinsèquement. 25 On définit donc à ce stade un graphe de liaison (ensemble des liaisons entre les mailles) non orienté. Si l'on peut ordonner toutes les mailles du graphe de liaison (orienté par les flux d'eau) de telle manière qu'une maille ne soit alimentée que par des mailles précédentes, on parle de graphe de liaison orienté sans circuit comme illustré figure 4. Sur cette figure, les mailles sont illustrées par des rectangles, les 30 cercles représentent les points d'injection de gaz, et les numéros représentent les niveaux de chacune des mailles. Ces niveaux sont définis de la façon suivante : une maille ne peut recevoir un flux que d'une maille de niveau inférieur. Ainsi le niveau 1 correspond à des mailles n'ayant que des flux sortant (injection), les mailles de niveaux 2 ne peuvent recevoir des flux que des mailles de niveau 1, ... les mailles de niveau 6 peuvent recevoir des flux de toutes les autres mailles mais aucun flux n'en sort. 2) Définition d'un modèle de transport réactif multiphasique 2.1- Les phénomènes à modéliser Les phénomènes intervenant dans ce modèle et que l'on souhaite modéliser sont : A) Les déplacements de fluides en milieu poreux sous la forme de 2 phases (eau, gaz) : Les déplacements des fluides (eau et gaz) sont régis par la loi de Darcy généralisée couplée à un modèle de pression capillaire. B) Les échanges de matière entre phases et constituants : Les échanges de constituants entre la phase gazeuse et l'eau Le gaz et l'eau sont des phases à composition variable, elles peuvent échanger directement des constituants pour vérifier des contraintes d'équilibre thermodynamique physique (changements de phase). - Les échanges d'éléments entre l'eau et les minéraux La roche est constituée de plusieurs phases minérales pouvant échanger de la matière avec l'eau par des réactions chimiques de surface : précipitation/dissolution. Les échanges d'éléments entre constituants au sein de l'eau Au sein de l'eau, les constituants peuvent se dissocier ou se combiner afin de vérifier 25 des contraintes d'équilibre thermodynamique chimique (réactions de spéciation, loi d'action des masses). C) Les échanges de matière avec l'extérieur du système Le système étudié n'est pas un système fermé. De la matière transite du système vers l'extérieur ou inversement par des puits ou des frontières artificielles du système. 2.2- Le modèle de transport réactif multiphasique complet Le modèle de transport réactif multiphasique complet est basé sur des équations de conservation de la masse pour un grand nombre de constituants dans plusieurs phases. On adopte dans cette description détaillée le point de vue des schémas volumes finis. On écrit les équations de conservation comme un bilan sur un volume élémentaire appelé maille. On fait alors apparaître des termes d'échange sous la forme de flux au travers de liaisons liant les différentes mailles, de termes sources liant une maille à l'extérieur ou de termes de réaction liant les constituants au sein d'un même maille. Afin de simplifier l'écriture des équations du modèle, on introduit quelques 15 notations : ~M : volume de la maille M (m3), : porosité du milieu, Pg : densité molaire de la phase gaz, Sg : saturation de la phase gaz, 20 Y : fraction molaire du constituant i dans le gaz, At" : pas de temps entre deux instants t" et (At" = t"+' ù t" ), Fg cy : flux d'échange de gaz entre deux mailles, Fg,b : flux d'échange de gaz entre une maille et l'extérieur, Tg' : transfert de matière d'un constituant i du gaz vers l'eau, 25 pW : densité molaire de la phase eau, SW : saturation de la phase eau, W : fraction molaire du constituant i dans l'eau FW : flux d'échange d'eau entre deux mailles, FH,,b : flux d'échange d'eau entre une maille et l'extérieur, Tw->g transfert de matière d'un constituant i de eau vers le gaz, Tw->w : transfert de matière d'un constituant i au sein du milieu aqueux, T"s : transfert de matière d'un constituant i de l'eau vers le solide, Ci concentration du constituant i, pin : densité molaire du minéral, fraction volumique du minéral m, Ti"' transfert de matière d'un constituant i d'un minéral vers l'eau, ni : nombre de mole du constituant i (mol.m-3) , JZ,a : terme de flux du constituant i (mol.s-1), Qi b : terme source du constituant i (mol.s"1), R1 : terme de réaction du constituant i (mol.s-1). Remarque. Pour un constituant présent en phase gazeuse et en phase aqueuse on définit deux constituants chimiques différents. Le changement de phase peut ainsi 15 s'écrire comme une réaction entre deux constituants du système. L'équation générale de conservation d'un constituant chimique quelconque intégrée sur un volume de maille M entre deux instants t" et tn+1 peut s'écrire : I44+l _nn)+AtnIJi v +AtnIQi b +OtnR1 = 0 6 b Les quatre termes de l'équation (2) sont nommés respectivement terme d'accumulation, terme de flux, terme source et terme de réaction. Pour fermer l'équation (2) pour tous les constituants i, il faut spécifier ou éliminer les termes de flux, de bord et de réaction. On peut séparer les constituants en 25 trois catégories: constituants gazeux (a(i) = g), constituants dans l'eau (a(i) = w) et constituants minéraux (a(i) =m). Sous les hypothèses A) B) et C) décrites précédemment, on peut détailler les différents termes de l'équation (2) : 20 (2) Jio- =(Pg.Y), .Fg, , Si a(l) = g ~ Qi,b ù l/'g •~ !b •Fg,b Ri = Tg" Ji,o- = (Pw W )6 Fw,Q Si a(i) = w Qi,b = `Pw •Wi lb •Fw,b Ri = T"g +T "w +T "s Si a(i) = m Les termes Fg c , Fg b , Fw a et Fw b sont des flux volumiques (m3. s-I). Vu d'une maille du système, les frontières (indice b) et les liaisons entre mailles (indice o) ont donc une représentation tout à fait similaire. Les termes Tg->w T W )g , T V W, T' et Tm*w sont des termes de réaction, ils dépendant d'autres constituants que i dans la même maille. On ne distingue pas ici les réactions physiques (changements de phase) et les réactions chimiques. Pour Ti'', tous les constituants sont dans la phase eau (réactions homogènes). Pour les autres termes réactifs, les constituants qui interviennent sont dans deux phases (réactions 15 hétérogènes). 2.3-Découpage et couplage Le traitement simultané de tous ces phénomènes entraîne inévitablement la résolution d'un système d'équations de grande taille avec un grand nombre 20 d'inconnues et de sous modèles. Pour résoudre ce problème, la méthode repose sur l'association de plusieurs sous modèles couplés de manière à obtenir une bonne approximation des phénomènes physicochimiques complets décrits ci dessus. On réalise, selon l'invention, un découpage simple de ces phénomènes en deux sous modèles. La résolution se fait alors en deux étapes, et chacune d'elles utilise l'opérateur complet dans lequel on remplace certains termes par des termes prédits en fonction des étapes précédentes. On estime que les termes choisis pour ces simplifications et les prédicteurs proposés permettent d'approcher correctement le problème posé dans le cadre du stockage géologique de gaz acide sans itération supplémentaire. La méthode fournit également le moyen de coupler ces deux sous modèles pour résoudre le problème complet (équation (2)). Cette méthode permet de prendre en compte l'ensemble des phénomènes sans itérations supplémentaires tout en conservant une bonne approximation du phénomène complet. Ce découpage a également pour conséquence de permettre l'utilisation d'une méthode de type lignes de courant qui a de nombreux avantages, expliqués dans la suite. Le découpage des phénomènes physicochimiques est réalisé selon un principe de taux de transfert . Un taux de transfert représente la quantité de matière transférée par unité de temps, d'un point à un autre d'un système géochimique. Ainsi, un taux de transfert correspond soit à un terme de flux , s'il s'agit d'un transfert entre deux mailles connectées par une liaison, soit un terme source , s'il s'agit d'un transfert entre une maille et l'extérieur du réservoir, soit à un terme de réaction s'il s'agit d'un transfert entre constituants au sein d'une même maille. Selon la méthode, les taux de transfert sont découpés en deux parties, comme l'illustre la figure 2: une partie traitée par un premier sous-modèle de type réservoir multiphasique compositionnel (RESER VOIR), une autre par un second sous-modèle de type transport réactif monophasique (TREACT) dans un champs de vitesse donné, les deux sous-modèles étant couplés par un coupleur (COUPL). Sous-modèle de type réservoir multiphasique compositionnel : modèle RESER VOIR 20 Le modèle de type réservoir multiphasique compositionnel , est appelé RESER VOIR. Les taux de transfert traités par ce premier modèle concernent tout les transferts liés au gaz : 1. Termes sources de la phase gazeuse (en entrée (injection) et en sortie 5 (production)) 2. Termes de réaction de type changement de phase gaz-eau 3. Termes de flux de gaz au sein du réservoir Ainsi, pour un volume de maille gM et un constituant i, l'équation de 10 conservation résolue par le schéma numérique s'écrit alors : IMI n+1,(1)p;+1,(1)SS+1,(1)yn+1,(I) _IMI on,(I)Pg,(1)Sg,(1)vn,(1) + Atn,(1) E~ (pg y )n: 1,(1)Ig ++1,(1) + At n,(1) Eb (Pg y )nn+1,(1)Fg +1,(1) +Atn,(1)Tg->w,n+o)) IMI ~n+I,(1)Pnv+1,(1)~,w+1,(])Wn+1,(1) _ + At n,(1) Ev (p w. ).+1,(1)F-1,(1) v1,(1) + At n,(1) E(PwPr; )b+1,(1)rw bl,(1) + Qtn,(1)Tw->g,n+l,(1) + Dtn,(1) (7.w-w,n,(z) + rr w~s,n,(z) ) Les transferts dans l'eau sont en fait calculés deux fois. Une fois dans le cadre 15 des transferts liés au gaz, en tant que prédiction, et une fois par le sous-modèle de type transport réactif monophasique (modèle TREACT). Pour tenir compte des modifications de porosité liées aux réactions chimiques on intègre une correction de volume sous la forme suivante : on+1,(1) _ n,(2) (5) Ces trois équations (3), (4) et (5) constituent le modèle de type réservoir multiphasique compositionnel (modèle RESERVOIR). Pour le gaz : Pour l'eau : r7in,(1)10w,(1) rw,(1)wn,(1) (3) =0 (4) =0 Les entrées de ce modèle RESERVOIR sont les suivantes : • Provenant d'un calcul RESERVOIR (indicées (1)) - O n'(1) : porosité à la date to - Sgn'(1) , Swn'(1) : saturations en gaz et en eau A. la date to Y n'(1) , W n'(1) : fractions molaires de constituants dans les phases eau et gaz à la date to . - pw (l) pg'(') : densités molaires des phase eau et gaz à la date t, • Provenant d'un calcul TREACT (indicées (2)) - n'(2) porosité modifiée déduite à la date to T"w'n'(2), Tw_s'n'(2) : débits réactions chimiques eau-eau et eau-roche à la date t, . Les sorties de ce modèle RESERVOIR sont les suivantes (indicées (1)) : • Sorties standard -Atn'(1) : pas de temps entre to et to+r _ on+1,(1) : porosité à la date to+1 - n+1,(1) S, n+l,(1) : g ,wsaturations à la date to+] - yn+l,(1) W n+l,(1) ~ , , : fractions molaires dans les phases eau et gaz à la date to+1• pw+1,(1) , pg+1,(1) : densités molaires des phase eau et gaz à la date tn+I • Sorties plus spécifiques _ Fn+t,(1) Fn+l,(l) g.v w,Qflux volumiques par liaison phase eau à la date tn tn+1 - Fn+l,(1) Fn+1,(1) : g,b 'w,bflux volumiques par frontière phase eau à la date Twù>g,n+l,(I) T gù>w,n+l,(1) , . . .taux de transfert réactions eau-gaz à la 15 date to+1 Sous-modèle de type transport réactif monophasique : modèle TREACT Ce modèle est appelé TREACT. Les taux de transfert traités par le second modèle, modèle de type transport réactif monophasique , concernent tout les 10 transferts de matière liés à l'eau et à la roche : 1. Termes de réaction de type précipitation/dissolution entre eau et minéraux 2. Termes de réaction de type spéciation au sein de la phase eau 3. Termes de flux d'eau au sein du réservoir 4. Termes source de la phase eau IMI n+1,(1),SW+1,(1)Ct +1,(2) - MI ~n,(1),sw,(1)Cn,(2) +Otn'(1)~ (C.)n+1,(2)r*n+1,(1) ~ r Q w +Atn,(1) Eb (Ci )b+1 (2)Fw,bl,(1) + At n,(1)T w~g,n+I,(1) + At n,(1) (7.wù)w,n+1,(2) + Twù>s,n+l,(2) ) t r Pour l'eau : =0 (6) IMI n+1,(2) Pm ù I MI ~ m(2)Pm + Atn,(1)Tm-mw,n+1,(2) Pour la phase minérale : (7) =0 20 Et l'on ajoute : n+1,(2) = 1 _ E n+1,(2) m m (8) Ces trois équations (6), (7) et (8) constituent le modèle de type transport réactif monophasique (modèle TREACT). Les entrées de ce modèle TREACT sont les suivantes : • Provenant d'un calcul RESER VOIR (indicées (1)) - Atn'(1) pas de temps entre t, et to+1 - n,(1) n+1,(1) : porosités aux dates to et to+1 _ S n,(I) n+l,(1) : w ~ wsaturations d'eau aux dates to et to+1 Fw ~'(1) : flux volumiques par liaison phase eau à la date t,~+1 Fw 61'(1) : flux volumiques par frontière phase eau à la date tä+1 - Tw->g,n+',(1) : taux de transfert réactions eau-gaz à la date tn+1 • Provenant d'un calcul TREACT (indicées (2)) - Cn'(2) : concentrations des constituants dans l'eau (molarités) à la date t, - pm : densité des minéraux fraction volumique minéraux à la date tz Les sorties de ce modèle TREACT sont les suivantes (indicées (2)): • Sorties standard - G,ln+1,(2> : concentrations des constituants dans l'eau (molarités) à la date to+l n+1,(2) : fractions volumiques des minéraux à la date to+1 • Sorties plus spécifiques -T.w-3w,n+1,(2) _ Tw-+s,n+l,(2) - (Dn+l,(2) : taux de transfert réactions eau-eau à la date tä+1 : taux de transfert réactions eau-roche à la date tn+1 porosité modifiée déduite à la date tn+1 - Tm-->w,n+l,(2) roche à la date to+r : détail par minéral des taux de transfert réactions eau- Le modèle de transport réactif multiphasique complet est traité grâce à ces deux sous-modèles, de manière à répartir les phénomènes traités. La connaissance de l'ensemble des termes de transfert, à l'aide de ces deux sous modèles, permet de connaître précisément l'état du système à un instant quelconque t >0, étant donné son état initial à la date t=0. Autrement dit, ces termes de transfert déterminent entièrement l'évolution du système géochimique étudié. 3) Acquisition des données nécessaires au modèle Les données d'entrée du modèle RESERVOIR pour l'initialisation sont obtenues soit par des tests in situ, comme des campagnes sismiques et des tests de puits, ou par des mesures directes sur échantillon (carottes). Ces données sont de 15 trois types : caractérisation du milieu poreux (perméabilité, porosité, faciès...) caractérisation des fluides (données thermodynamiques des fluides) caractérisation croisée (perméabilités relatives) Les données d'entrée du modèle TREACT pour l'initialisation sont des données 20 issues d'expériences de laboratoire et de mesures calorimétriques associées. Les données cinétiques sur les réactions de précipitation/dissolution sont moins bien connues que les données de type thermodynamique elles doivent souvent être ajustées au cas par cas. 25 4) Évaluation du transport de constituants et des processus chimiques De façon générale, la méthode selon l'invention permet de déteuniner d'une part, les processus hydrodynamiques conduisant au transport des constituants chimiques, et d'autre part, les processus chimiques intervenant dans un système géochimique. La méthode selon l'invention est basée sur un schéma d'évolution en temps. On utilise un modèle reproduisant l'évolution au cours du temps de phénomènes physiques précis. Comme on ne peut reproduire une période de temps complète en une fois, on la découpe en un nombre fini de sous temps qu'on appelle pas de temps (évolution en temps discret). La longueur de ces pas de temps n'est pas constante : la période est découpée en pas de temps discrets dont la taille peut changer au cours de la simulation. Chaque étape fait avancer la simulation d'un pas de temps strictement positif et dans une simulation normale, le temps final (Tf) est atteint en un nombre d'étapes inférieur à (Tf ù To )/A""n , avec To le temps initial et 10 Ar" le pas de temps minimal autorisé. Ce pas minimal est choisi pour limiter le nombre d'étape (une étape par pas de temps), et constitue ainsi un critère d'arrêt. Une fois l'initialisation de chaque donnée d'entrée des modèles réalisées (par mesures, par des modélisation antérieures, .), la modélisation des processus de transport et de réactions chimiques au sein du système géochimique, se fait 15 schématiquement en cinq grandes étapes pour chaque pas de temps : A- Synchronisation des deux sous modèles B- Résolution du sous modèle RESERVOIR C- Bilan RESERVOIR et préparation du calcul TREACT D- Résolution du sous modèle TREACT 20 E- Bilan TREACT et préparation du calcul RESERVOIR (pas de temps suivant) Étape A : synchronisation des sous modèles La synchronisation de début de pas de temps consiste à estimer un nouveau pas 25 de temps pour le calcul suivant. Cette étape permet également d'effectuer les opérations nécessaires pour remettre les modèles dans un état leur permettant de refaire les étapes B, C, D, et E. Si le pas de temps minimal est atteint, on arrête le programme en renvoyant une erreur, la simulation a échoué. Étape B : Résolution du sous modèle RESERVOIR Dans cette étape on remplace les termes de réaction chimique (eau-eau et eau-roche) par une prédiction calculée en tenant compte des calculs TREACT précédents 5 (au pas de temps précédent). Le système à résoudre est défini par les équations (3), (4) et (5). A priori, ce problème est de type standard et il peut être résolu par tout sous-modèle de type transport compositionnel multiphasique non réactif (cf. [6]). 10 Étape C : Bilan du modèle RESERVOIR et préparation du calcul TREACT Les données issues de la résolution du sous modèle RESERVOIR sont utilisées pour prédire un nombre important de termes qui apparaissent dans le sous-modèle TREACT. Il s'agit du flux volumique d'eau sur les liaisons et les frontières et des termes de réaction eau-gaz. On utilise également les volumes d'eau et de gaz prédits 15 par le sous-modèle RESER VOIR...DTD: Étape D : Résolution du sous modèle TREACT Cette étape s'apparente à une étape d'écoulement réactif avec une seule phase mobile (l'eau) dans un champs de vitesse donné. Pour une maille M et un constituant 20 i, l'équation de conservation résolue par le schéma numérique s'écrit selon les équations (6) et (7). A priori, ce problème est de type standard et il peut être résolu par tout sous-modèle de type transport réactif monophasique . Cependant, selon un mode de réalisation particulier de l'invention on peut 25 définir un sous-modèle de type transport réactif monophasique à `partir d'un modèle géochimique existant. Une telle résolution du modèle de type transport réactif monophasique (TREACT) dans un champs de vitesse donné est proposée et illustrée sur la figure 3. Selon ce mode de résolution, on effectue des calculs de géochimie maille par maille à l'aide d'un modèle de type géochimique (GEOCI-I) 30 tout en faisant appel à une analyse du graphe de liaisons (FA) pour effectuer ces calculs selon un principe de lignes de courant (LC). Ces calculs fournissent les sorties du modèle de type transport réactif monophasique (TREACT). Plus en détail, le modèle TREACT comporte deux grandes étapes : - L'analyse de flot (FA) qui correspond à la recherche de l'ordre de parcours des mailles sur le graphe orienté. Cette étape prend donc en entrée le graphe orienté par les flux et donne en sortie si oui ou non l'écoulement est de type potentiel, et si l'écoulement est de type potentiel l'ordre à utiliser pour parcourir les mailles. - La distribution du calcul de chimie (LC) : pour chaque maille, on prépare le calcul en récupérant les flux entrants et sortant de cette maille, et on lance le calcul géochimique (GEOCH). On enregistre les résultats. Le modèle GEOCH est appelé autant de fois qu'il y a de mailles, et il travaille sur des données différentes à chaque appel. Cette méthode a deux avantages : d'une part, il existe de nombreux modèles géochimiques et la méthode peut donc s'appliquer de façon indépendante au choix de celui-ci. D'autre part, ce modèle utilise les propriétés spécifiques du problème : c'est une méthode basée sur une analyse d'un graphe de liaison et une distribution des calculs de géochimie sur des lignes de courant. La méthode est décrite ci dessous. Le schéma tel qu'il est écrit ci dessus est a priori un schéma implicite, les termes de flux, les termes sources et les termes de réactions font intervenir des variables à la date tn+1 qui sont des inconnues du problème. Les méthodes de type implicite sont connues pour être inconditionnellement stables. Ceci peut être prouvé sur des exemples académiques de manière rigoureuse. La pratique a prouvé qu'en effet le fait d'écrire un schéma implicite stabilise le calcul dans de nombreux cas. La mise en oeuvre d'un schéma implicite demande d'utiliser une méthode itérative pour la résolution d'un systèmes d'équations non linéaires de grande taille. Habituellement cette méthode demande de dériver l'ensemble des termes des équations de conservation par rapport aux variables principales choisies. L'idée est ici d'utiliser les propriétés du champs de vitesses de l'eau de type champs dérivant d'un potentiel pour éviter cette résolution coûteuse. Comme 30 exemple de champs qui dérive d'un potentiel, on peut citer les champs suivants. Le champs électrique E dérive du potentiel électrique V. En effet, on peut écrire : E = ùgrad(V) Le champs gravitationnel sur terre G dérive également d'une énergie 5 potentielle Ep. En effet, on peut écrire : G = ùgrad (E p ) La vitesse de Darcy Vw qui sert à calculer le flux volumique pour les liaisons et les frontières peut s'écrire sous la forme suivante : Vw = K * kY * (grad (P) + p, * G) ~w avec : Vw vitesse de Darcy de la phase eau K : perméabilité du milieu poreux kr : perméabilité relative 15 ,uw : viscosité P : pression pw masse volumique G : vecteur gravité Elle n'est pas sous la forme : Vr, = grad (Q) , où Q serait un potentiel. 20 Cependant, par une transformation usuelle (cf [7]), lorsque la densité de l'eau est peu variable, on peut écrire : Vw =KH * grad (H) H étant la charge hydraulique et KH la conductivité hydraulique. Par conséquent, même si Vw ne dérive pas directement d'un potentiel, sa 25 direction est proportionnelle à grad (H) qui lui dérive d'un potentiel (la charge 10 hydraulique). C'est suffisant pour admettre que le déplacement se fait dans le sens des potentiels de charge hydraulique décroissants. On parlera par la suite d'un écoulement de type potentiel. On présume les propriétés de faible variabilité de la densité vérifiées et donc que l'écoulement de l'eau est de type potentiel. En pratique, on propose une implémentation implicite du modèle TREACT dans les cas où la propriété suivante est vérifiée : on peut ordonner toutes les mailles du graphe de liaison (orienté par les flux d'eau) de telle manière qu'une maille ne soit alimentée que par des mailles précédentes. Dans la théorie des graphes orientés, on dit que le graphe orienté est sans circuit comme illustré figure 4, où les numéros représentent les niveaux de chacune des mailles (rectangle), et les cercles représentent les points d'injection. La figure 5 en revanche, illustre le cas d'un graphe orienté avec circuit. On vérifie a posteriori que cette propriété est vraie mais on ne calcule pas H ni aucun autre potentiel. Le sens de parcours des mailles est rapidement déterminé et il permet deux choses : -triangulariser la partie du système (bloc extra-diagonal) qui lie les mailles entre elles. Conséquence : la résolution d'un grand système non linéaire se résume à faire un calcul sur une maille puis à remplacer le résultat dans les équations des autres mailles et ceci pour chaque maille dans l'ordre de parcours déterminé. - le regroupement de mailles en niveaux permet également de les traiter en parallèle. Les calculs de toutes les réactions ayant lieu sur les mailles d'un même niveau peuvent être faits de manière indépendante donc sur des processeurs ou des machines différentes. En pratique, sur les cas réels, les niveaux sont très importants (une centaine de mailles), on peut donc distribuer au mieux les calculs. C'est le deuxième avantage de cette méthode. En présence de fortes variations de densité, qui peuvent être occasionnées par une dissolution importante de gaz ou de sels minéraux dans l'eau, et pour conserver un algorithme de distribution de la chimie, on propose d'utiliser uniquement des réactions implicites et des flux explicites. On obtient ainsi un schéma semi-implicite. Les calculs de maille deviennent indépendants et on peut paralléliser l'ensemble du calcul. Étape E : Bilan du modèle TREA CT et préparation du calcul RÉSERVOIR Les données issues de la résolution du sous modèle TREACT sont utilisées pour prédire des termes qui apparaissent dans les équations du modèle RESER VOIR. Il s'agit des variations de porosité et des transferts de type réaction chimique eaueau et eau-roche. Pour calculer la variation de porosité dans le sous modèle TREACT, on effectue une simple différence : n+1,(2) = 1 _ E nm+l,(2) m Pour calculer le bilan des transferts liés aux termes de réaction eau-eau et eau-15 roche (T,.w w n,(2) +Tw,s,n,(2)), on a la même difficulté que lors du calcul des termes de réaction eau-gaz à l'étape C car là aussi on peut supposer que certaines réactions sont instantanées. On procède de la même façon : une fois que les variables au temps tn+1 ont été calculées, on utilise l'équation (6) où tous les autres termes sont connus pour en déduire la valeur du terme recherché. 20 5) Évaluation des capacités de stockage de gaz acides du milieu Appliquée à un réservoir géologique, la méthode permet de déterminer les capacités de stockage en gaz acides de ce réservoir, comme le CO2 par exemple. En effet, les phénomènes modélisés sont : 25 - les déplacement de fluides en milieu poreux sous la forme de 2 phases (eau, gaz) - les échanges de constituants entre la phases gazeuse et Peau (changement de (9) phase) - les échanges de constituants entre l'eau et les minéraux (précipitations, dissolutions) - les échanges d'éléments entre constituants au sein de l'eau (réactions de spéciation) -les échanges de matière avec l'extérieur du système (puits, conditions limites) Ainsi, à l'issue des calculs provenant des différents modèles on connaît les propriétés suivantes : - Atn : pas de temps entre to et to+1 10 - (Dn+1 : porosité à la date to+1 Sgn+l Swn+' : saturations à .la date to+1 n+1 - Y,. date to+1. : fractions molaires dans les phases eau et gaz à la n+1 n+1 n+1 - Pw ~ Pg: densités molaires des phase eau et gaz à la date to+1 15 - Fg a' Fw v1 : flux volumiques par liaison phase eau à la date t, +1 Fg b' , Fg b' : flux volumiques par frontière phase eau à la date to+1 - TW->g,n+1 taux de transfert réactions eau-gaz à la date tn+i concentrations des constituants dans l'eau (molarités) à la 20 date tn+r - gym' : fractions volumiques des minéraux à la date tn+1 T W*W TW>S : taux de transfert réactions eau-eau et eau-roche On+1,(2) porosité modifiée déduite à la date tä+1 On peut alors déterminer, pour tout instant t, la quantité de gaz présente dans le réservoir, la quantité de gaz qui précipite au cours de la période de temps t, la quantité de gaz qui se dissout dans la phase aqueuse au cours de la période de temps t, la quantité qui migre au sein du réservoir, la quantité qui migre dans le sous-sol hors du réservoir et la quantité qui migre en surface. On peut également déterminer le niveau d'altérations du réservoir par rapport aux réactions de précipitation et de dissolution. Le modèle peut alors être utilisé de deux façons : - Phase d'injection Lors de la phase d'injection de gaz acides dans le réservoir géologique, la méthode permet de définir une stratégie d'injection et donc d'optimiser le remplissage du réservoir. Parmi les paramètres ci-dessus, la porosité et les taux de réaction eau roche 15 peuvent, par exemple, permettre de suivre les colmatages au niveau du puits. - Phase de stockage Lors de la phase de stockage, la méthode est prédictive. Elle permet d'étudier la stabilité du réservoir au cours du temps et donc d'anticiper les problèmes de fuite (impact sur l'environnement). En effet les réactions de précipitations et de 20 dissolution peuvent altérer la qualité du réservoir sur le long terme. La dégradation de la qualité du réservoir peut alors devenir un frein au stockage voir à un arrêt du stockage. Il est donc très important de dételininer la quantité de gaz qui précipite, la quantité qui migre au sein du réservoir, dans le sous-sol ou à la surface. La méthode permet en particulier de localiser et déterminer la quantité de CO2 25 présente dans le réservoir à un instant t, sous forme de gaz, dissous dans l'eau ou précipité (carbonates). Parmi les paramètres ci-dessus, la composition et les volumes des phases (gaz, eau et minéraux) sont particulièrement intéressants pour l'ingénieur. La méthode peut également trouver des applications directes dans les études de dépollution des sols par des procédés multiphasiques ("venting", utilisation de solvants chlorés). La méthode selon l'invention permet donc de déterminer les capacités de stockage de gaz acides (H2S, CO2, CH4) dans un milieu souterrain par modélisation des processus de transport et de réactions chimiques en milieu poreux. Du point de vue de la modélisation, la méthode tire parti des avantages des méthodes séquentielles (pas de temps locaux, parallélisation, taille mémoire réduite) en évitant d'ajouter une boucle itérative supplémentaire. Les avantages sont alors les suivants : diminution du temps de calcul total d'une simulation par rapport à une méthode séquentielle itérative en gardant la possibilité d'effectuer des calculs en parallèle maille par maille ; approche modulaire permettant de créer un modèle complet à partir de sous modèles RESERVOIR et TREACT existants plus limités ; implémentation facilement parallélisable du modèle TREACT basée sur une méthode d'analyse de flot et de distribution de calculs OD de type géochimique sur des lignes de courant ; Du point de vue du stockage de gaz acides la méthode fournit un outil indispensable au spécialiste pour modéliser les processus de transport et de réactions chimiques d'un milieu poreux souterrain, en donnant une description précise du comportement dynamique des fluides dans le milieu géologique tout en simulant la plus grande partie des processus chimiques qui interviennent en phase aqueuse. La méthode effectue un traitement complet du problème compositionel multiphasique représentant un stockage de gaz acide (CO2, CH4, H2S). REVENDICATIONS 1) Méthode pour déterminer les capacités d'un milieu géologique souterrain à stocker au moins un gaz acide, ledit milieu comprenant une phase solide et une phase 5 aqueuse, caractérisée en ce que la méthode comporte les étapes suivantes : a) on découpe ledit milieu en un ensemble de sous volumes appelés mailles ; b) on définit un modèle de transport réactif multiphasique, en couplant un sous modèle de type réservoir multiphasique compositionnel , dédié à la 10 détermination des taux de transfert liés au gaz, avec un sous modèle de type transport réactif monophasique , dédié à la détermination des taux de transfert liés à ladite phase aqueuse et à ladite phase solide ; c) on réalise une acquisition de données pour initialiser ledit modèle de transport réactif multiphasique ; 15 d) on détermine, pour une période de temps t donnée, des processus de transport de constituants chimiques et des processus de réactions chimiques au sein dudit milieu géologique, en résolvant par couplage ledit modèle de transport réactif en toute maille dudit milieu, pour obtenir au moins : les quantités de gaz présent dans ledit milieu géologique, ayant migré 20 hors dudit milieu et en surface, ayant précipité, et dissout dans la phase aqueuse; - les altérations de la phase solide par précipitation et dissolution. e) on en déduit les capacités dudit milieu géologique à stocker ledit gaz acide. 25 2) Méthode selon la revendication 1, dans laquelle on résout ledit modèle de transport réactif sur ladite période de temps t, à l'aide d'un schéma d'évolution en temps, dans lequel on découpe ladite période de temps en un nombre fini de pas de temps et en réalisant les étapes suivantes :a) on définit la durée d'un pas de temps pour synchroniser les deux sous modèles ; b) on résout ledit sous modèle réservoir multiphasique compositionnel ; c) on résout ledit sous modèle ((transport réactif monophasique à l'aide des résultats obtenus à l'étape b) ; d) on modifie la durée dudit pas de temps et l'on réitère à partir de l'étape b) en tenant compte des résultats obtenus à l'étape c), jusqu'à couvrir ladite période de temps t ; 3) Méthode selon l'une des revendications précédentes, dans laquelle lesdits taux de transfert liés au gaz comportent au moins l'un des termes suivants : termes sources dudit gaz ; termes de changement de phase entre le gaz et la phase aqueuse ; termes de flux de gaz au sein dudit milieu. 4) Méthode selon l'une des revendications précédentes, dans laquelle lesdits taux de transfert liés à ladite phase aqueuse et à ladite phase solide comportent au moins l'un des termes suivants : termes de réactions de précipitation et de dissolution entre la phase aqueuse et la phase solide ; termes de réaction de spéciation au sein de la phase aqueuse ; termes de flux de la phase aqueuse au sein dudit milieu ; - termes source de la phase aqueuse. 5) Méthode selon l'une des revendications précédentes, dans laquelle on détermine à partir dudit modèle de transport réactif multiphasique au moins l'un des 25 phénomènes suivants : les déplacements de ladite phase aqueuse et dudit gaz ; les échanges de constituants entre la phase aqueuse et le gaz ; les échanges de constituants entre la phase aqueuse et la phase solide ;20- les échanges de constituants entre constituants au sein de la phase aqueuse ; - les échanges de constituants avec l'extérieur dudit milieu. 6) Méthode selon l'une des revendications précédentes, dans laquelle on résout le modèle de type transport réactif monophasique de manière implicite à partir d'un sous modèle de type réacteur géochimique OD ouvert en réalisant les étapes suivantes : - on définit des possibilités d'échanges de constituants entre lesdites mailles, définissant ainsi un graphe de liaison ; on détermine un ordre de parcours desdites mailles en analysant ledit 10 graphe de liaison, de façon à parcourir les mailles dans l'ordre des potentiels décroissants ; - on parcourt chaque maille selon ledit ordre de parcours et on résout ledit modèle de type réacteur géochimique pour chaque maille, à l'aide d'une technique de ligne de courant. 15 7) Méthode selon la revendication 6, dans laquelle on résout ledit modèle de type transport réactif monophasique de manière semi-implicite lorsque l'on ne peut déterminer l'ordre des potentiels décroissants.