FR2823877A1

FR2823877A1 - Methode pour contraindre par des donnees dynamiques de production un modele fin representatif de la repartition dans le gisement d'une grandeur physique caracteristique de la structure du sous-sol

Info

Publication number: FR2823877A1
Application number: FR0105339A
Authority: FR
Inventors: Mokhles Mezghani; Frederic Roggero
Original assignee: IFP Energies Nouvelles IFPEN
Current assignee: IFP Energies Nouvelles IFPEN
Priority date: 2001-04-19
Filing date: 2001-04-19
Publication date: 2002-10-25
Anticipated expiration: 2021-04-19
Also published as: US6662109B2; GB2375635A; FR2823877B1; NO20021824L; NO334984B1; NL1020410C2; US20030028325A1; GB0208563D0; GB2375635B; NO20021824D0

Abstract

- Méthode pour une mise à jour directe, par des données dynamiques de production, d'un modèle géologique fin représentatif de la répartition dans un gisement d'une grandeur physique caractéristique de la structure du sous-sol.- La méthode est basée sur le couplage des techniques d'inversion et de changement d'échelle (up-scaling) afin de permettre l'optimisation des paramètres pétrophysiques d'un modèle de simulation plus grossier issu du modèle géologique fin. Elle comporte le paramétrage direct du modèle géologique fin, le changement d'échelle n'étant qu'un moyen pour obtenir rapidement une approximation des résultats de simulation et de leurs dérivées par rapport à la paramétrisation du modèle géologique fin.- Applications à la détermination d'un schéma d'exploitation qui optimise la production d'un gisement d'hydrocarbures par exemple.

Description

La présente invention concerne une méthode pour contraindre par les données dynamiques de production un modèle géologique fin représentatif de la répartition dans un réservoir hétérogène, d'une grandeur physique caractéristique de la structure du sous-sol, telle que la perméabilité ou la porosité.

L'objectif visé est de fournir aux ingénieurs de réservoir une méthodologie permettant une mise à jour efficace des modèles géologiques au fur et à mesure de l'acquisition des données dynamiques.

Etat de la technique
L'état de la technique auquel on fera référence ci-après est décrit dans les publications suivantes : - Wen, X-H, et al "Upscaling hydraulic conductivities in heterogeneous media: An overview. Journal of Hydrology (183)", ix-xxxii, 1996.

- Renard P., "Modélisation des écoulements en milieux poreux hétérogènes : calcul des perméabilités équivalentes". Thèse, Ecole des Mines de Paris, Paris, 1999.

- G. de Marsily. De l'identification des systèmes hydrologiques. Thèse, Université
Paris 6, Paris, 1976.

- Hu L. -Y. et al. "Constraining a Reservoir Facies Model to Dynamic Data Using a Graduai Déformation Method". VI European Conference on the Mathematics of Oil Recovery, Peebles, 1998.

- Tarantola. A. "Inverse Problem Theory: Method for Data Fitting and Model
Parameter Estimation". Elseiver, Amsterdam, 1987.

Anterion F. et al. "Use of Parameter Gradients for Réservoir History Matching".

SPE 18433, Symposium on Réservoir Simulation of the Society of Petroleum
Engineers, Houston, 1989.

- Wen X.-H. et al. "High Resolution Réservoir Models Integrating Multiple-Well
Production Data". SPE 38728, Annual Technical Conference and Exhibition of the Society of Petroleum Engineers, San Ahtonio, 1997.

- Chu L. et al., Commutation of Sensitivity Coefficients With Application to the
Intégration of Static and Well Test Pressure Data". Eclipse International Forum, Milan, 1994.

Pour conduire l'exploitation d'un réservoir et prévoir son comportement dynamique en fonction des différents scénarii de production, l'industrie pétrolière fait un large usage de simulations numériques de modèles d'écoulements. Les modèles géostatistiques utilisés pour représenter la structure géologique du réservoir (perméabilité, porosité, etc.) requièrent une discrétisation constituée d'un nombre important de mailles pouvant atteindre une dizaine de millions.

Pour être en mesure de réaliser des simulations numériques d'écoulement avec des temps de calcul raisonnables, la pratique usuelle consiste à construire un modèle de simulation grossier, en regroupant des mailles de propriétés différentes et en affectant aux macro-mailles une propriété équivalente calculée à partir des propriétés locales. Cette opération est appelée changement d'échelle.

La caractérisation contrainte de réservoir a pour but de déterminer les paramètres du modèle de simulation afin que ce dernier puisse reproduire l'historique de production du gisement qu'il doit modéliser. Cette phase d'estimation des paramètres est encore appelée calage d'historique de production. Ainsi, le modèle de simulation d'écoulement est-il compatible avec l'ensemble des données statiques et dynamiques disponibles.

Dans la pratique courante, l'estimation des paramètres du modèle de simulation est effectuée par une série d'essais et d'erreurs grâce au simulateur d'écoulement.

Le problème du calage d'historique de production peut aussi se formuler comme un problème de minimisation d'une fonction objectif mesurant l'écart entre l'historique de production observé sur le terrain et les prévisions du simulateur d'écoulement. La minimisation est alors réalisée en utilisant des techniques d'optimisation ou de contrôle optimal.

Une méthode pour prédire, par une technique d'inversion, l'évolution de la production d'un gisement souterrain et notamment d'un gisement recelant des hydrocarbures, est décrite par exemple dans le brevet FR 2 734 069 (US 5,764,515) du demandeur
Dès que les paramètres du modèle de simulation sont ajustés, ce modèle peut être utilisé pour simuler le comportement présent et futur du gisement. On peut alors disposer d'une évaluation des réserves in situ et déterminer un schéma d'exploitation qui optimise la production.

La caractérisation contrainte de réservoir fait donc appel à une multitude de techniques allant de la modélisation géostatistique aux problèmes d'optimisation.

L'introduction des principales techniques mises en #uvre dans le cadre de la méthodologie de couplage inversion et changement d'échelle fera l'objet de la suite de cette section.

Modélisation géostatistique
La géostatistique, dans sa présentation probabiliste, suppose qu'une variable spatiale, comme par exemple la perméabilité, peut être interprétée comme une réalisation particulière d'une fonction aléatoire, définie par sa loi de probabilité en tout point de l'espace. L'utilisation de plus en plus courante de la géostatistique par les compagnies pétrolières conduit à la construction de modèles fins pouvant atteindre un grand nombre de mailles. En effet, la géostatistique permet l'estimation des propriétés pétrophysiques dans l'espace à partir de mesures locales. En toute rigueur, la réalisation du modèle géostatistique doit être menée à l'échelle du support de la mesure, le modèle ainsi obtenu peut alors atteindre plusieurs millions de mailles. La puissance actuelle des ordinateurs ne rend pas envisageable la simulation numérique d'écoulement à l'échelle du modèle géostatistique. Pour réduire le nombre de mailles, il faut les regrouper, ce qui nécessite le calcul des propriétés équivalentes des nouvelles mailles en fonction des propriétés des mailles à petite échelle, une opération que l'on désigne par changement d'échelle ( ou upscaling).

Changement d'échelle
Le calcul de la perméabilité équivalente des milieux poreux hétérogènes est un sujet largement étudié par la communauté des géologues, des ingénieurs de réservoirs et plus généralement des physiciens des milieux poreux.

D'un point de vue mathématique, le processus de changement d'échelle de chaque perméabilité directionnelle peut être représenté par l'opérateur vectoriel F défini par :

k : la perméabilité à l'échelle du modèle géostatistique (dimension Rm);

K : la perméabilité à l'échelle du modèle de simulation d'écoulement (dimension ) ;
L'état de l'art des techniques existantes est dressé par Wen et al. (1997) et Renard (1999) déjà cités. Parmi les techniques connues de changement d'échelle, on peut citer les méthodes algébriques qui proposent des règles analytiques simples pour un calcul plausible des perméabilités équivalentes sans résoudre un problème d'écoulement. On peut opter par exemple pour la méthode connue dite de moyenne de puissance . La perméabilité K du bloc est égale à une moyenne de puissance, aussi appelée moyenne d'ordre w, dont l'exposant w est compris entre -1 et +1 :

Ainsi le problème de calcul de la perméabilité équivalente se ramène à l'estimation de l'exposant w permettant de minimiser l'erreur induite par le changement d'échelle (définie selon un certain critère). Pour les milieux à distribution log-normale isotrope et à faible longueur de corrélation, on sait que :

a étant la dimension en espace ([alpha] = 1, 2 ou 3).

On connaît aussi des techniques numériques de changement d'échelle où, pour calculer la perméabilité équivalente, on est amené à résoudre les champs de pression P et de vitesse v d'un problème d'écoulement local ou global :

Le terme désigne la viscosité du fluide en écoulement.

Paramétrisation
Le problème de la mise à jour des modèles géologiques par les données dynamiques repose sur la résolution d'un problème inverse. Ceci pose naturellement le problème de la paramétrisation du champ de perméabilité afin de rendre envisageable la minimisation de la fonction objectif qui mesure l'écart (au sens des moindres carrées) entre les données dynamiques observées sur le terrain et les résultats de simulation.

La paramétrisation des modèles géostatistiques est un point fondamental pour garantir le succès de la phase d'intégration des données dynamiques dans les modèles géologiques. En effet, cette intégration s'effectue dans le cadre d'un schéma itératif, régi par le processus d'optimisation, et visant à perturber un champ de perméabilité initial représentatif du modèle géostatistique considéré.

Dans l'idéal, le champ de perméabilité final doit non seulement respecter l'ensemble des données dynamiques prises en compte dans la fonction objectif, mais également préserver la cohérence géostatistique du modèle (moyenne, variogramme, etc. ) Le respect des données dynamiques est contrôlé par la fonction objectif dont la valeur est une évaluation de la qualité du calage. Quant à la cohérence des données géostatistiques, c'est la paramétrisation du champ de perméabilité qui va permettre de la gérer.

Une technique connue permettant d'effectuer cette paramétrisation est celle dite des points pilotes qui repose sur le principe de la simulation géostatistique conditionnelle appliquée aux modèles de type gaussien est décrite par exemple par de Marsily (1976) déjà cité.

Une autre technique connue permettant d'effectuer cette paramétrisation est celle dite des déformations graduelles. Telle que décrite par Hu et al. (1998) ainsi que dans les brevets FR 2 780 798, FR 2 795 841 et la demande de brevet FR EN 01/03194, du demandeur, la méthode des déformations graduelles consiste à

écrire qu'une nouvelle réalisation du champ de perméabilité à estimer, supposé de type gaussien, est une combinaison linéaire de réalisations indépendantes de la fonction aléatoire qui le modélise. Le champ de perméabilité * est donné donc par :

(#i )l#i#n : les coefficients de la combinaison linéaire ; et (ki)l#i#n : les réalisations indépendantes du modèle géostatistique considéré.

Afin de préserver les propriétés géostatistiques du modèle, les coefficients 0 doivent satisfaire la contrainte de normalité suivante :

Les coefficients 0 sont estimés afin que le champ de perméabilité résultant k(#) reproduise au mieux les données dynamiques.

Contrairement à la méthode des points pilotes, la méthode des déformations graduelles peut être appliquée de façon locale ou globale. Un respect strict des propriétés géostatistiques du modèle est assuré moyennant le respect de la contrainte de normalité (Equation 6) sans recours à l'introduction d'un modèle a priori dans la fonction objectif.

Fonction objectif
La mise à jour d'un modèle géologique par les données dynamiques repose sur la minimisation d'une fonction objectif qui mesure l'écart entre les données dynamiques observées sur le terrain et les résultats de simulation obtenus pour une valeur fixée des paramètres 8 .

Plusieurs formulations sont possibles pour définir une fonction objectif. La formulation au sens des moindres carrés est la plus couramment utilisée dans le secteur pétrolier. Ainsi la fonction objectif s'écrit comme suit :

avec : dobs :les données dynamiques observées sur le terrain ; D(#) : les résultats de simulation pour la valeur fixée des paramètres # ; Cd : la matrice de covariance sur les observations.

Comme décrit par Tarantola (1987), -une formulation mieux adaptée à la résolution des problèmes inverses mal posés, consiste à ajouter un terme de régularisation (modèle a priori) dans la fonction objectif :

avec : #pri : une estimation à priori des paramètres 19 ; C# : la matrice de covariance sur les paramètres.

Cette dernière formulation de la fonction objectif possède une interprétation probabiliste. En effet, dans le contexte de l'inversion bayésienne, le modèle a priori est donné par une densité de probabilité de loi quelconque.

Pour un modèle a priori de loi gaussienne de moyenne #pri et de covariance C#, cette densité de probabilité s'écrit :

Dans ce même contexte, la probabilité d'obtenir les observations connaissant la valeur des paramètres 19, ou fonction de vraisemblance, peut alors s'exprimer sous la forme suivante :

Lorsque l'opérateur de simulation d'écoulement D est linéaire par rapport aux paramètres #, la densité de probabilité a posteriori est encore de loi gaussienne.

La minimisation de la fonction objectif nécessite le calcul des dérivées des résultats de simulation par rapport aux paramètres à estimer. Soit :

Ce calcul des dérivées, primordial pour conduire le processus de minimisation dans les meilleures conditions, a fait l'objet de plusieurs travaux dont une synthèse est donnée dans l'article de Chu et al. (1994) déjà cité.

Jusqu'ici deux méthodes sont essentiellement rencontrées en industrie pétrolière : les gradients numériques et la méthode des gradients. Compte tenu de ses qualités en termes de stabilité numérique et de rapidité, nous avons adopté la méthode des gradients pour le calcul des dérivées des résultats de simulation par rapport à la paramétrisation du modèle géostatistique fin.

On rappelle que l'usage des lettres minuscules permet de référer au modèle géostatistique fin et l'usage des lettres majuscules permet de référer au modèle de simulation grossier. A titre d'exemple :

k désigne le champ de perméabilité à l'échelle du modèle géostatistique alors que K désigne le champ de perméabilité à l'échelle du modèle de simulation d'écoulement (après changement d'échelle) ; d désigne les résultats de simulation, obtenus à partir du modèle géostatistique fin, alors que D désigne les résultats de simulation obtenus grâce au modèle de simulation grossier (après changement d'échelle).

Calcul des dérivées - Méthode des gradients
La méthode des gradients permet de calculer les dérivées des résultats d'une simulation numérique d'écoulement par rapport à un certain nombre de paramètres qui interviennent dans le modèle de simulation. A titre d'exemple, il est possible de calculer les dérivées des principaux résultats de production (pression, saturation, débit, etc. ) par rapport aux propriétés pétrophysiques (perméabilité, porosité, etc.) assignées à des zones du réservoir.

La méthode des gradients est basée sur la dérivation des équations discrétisées du modèle d'écoulement telle que décrite par Antérion et al. (1989) déjà cité. Ces équations discrétisées se présentent sous la forme d'un système d'équations non-linéaires du type suivant :

# : les paramètres à estimer ;
Uint :l'initialisation des inconnues à simuler. Cette initialisation est calculée à partir des conditions initiales du système d'équations aux dérivées partielles qui modélisent l'écoulement ;
Ur : les inconnues de simulation calculés à l'instant tR ; Un+1 : les inconnues de simulation calculés à l'instant tn+1.

Le système (12) étant non-linéaire, sa résolution est effectuée généralement par la méthode de Newton basée sur des linéarisations successives du système nonlinéaire (12) selon le schéma suivant :

Le calcul des dérivées des résultats de simulation par rapport à la paramétrisation repose sur la dérivation directe du système (12). De cette dérivation résulte un nouveau système linéaire dont les inconnues sont les dérivées #Un+1 ##. Pour chaque paramètre ', ce système s'exprime sous la forme suivante :

La matrice du système linéaire est donnée par le terme :

II s'agit de la matrice de Newton du système (13) à l'itération finale. Le second membre de ce système linéaire est donné par le terme :

La résolution de ce système linéaire (un second membre par paramètre) permet d'obtenir l'ensemble des dérivées des inconnues de simulation par rapport à la paramétrisation souhaitée.

Par dérivation composée, il est possible d'exprimer les dérivées des principaux résultats de production par rapport à la paramétrisation :

Techniques d'optimisation
Les algorithmes d'optimisation non-linéaire permettent de calculer, selon un processus itératif, une valeur #opt des paramètres () qui minimise (localement ou globalement) la fonction objectif que nous souhaitons optimiser.

Les résultats de simulation issus de la distribution ' doivent permettre un meilleur calage des données dynamiques que ceux issus de la distribution initiale k(#(0)). #(0) désigne la valeur des paramètres () utilisés pour initier le processus d'optimisation.

L'itération (k+1) d'un tel algorithme d'optimisation a pour but de déterminer une nouvelle estimation des paramètres 0 selon le principe suivant :

- Calcul d'une direction : la direction s(k) est solution d'un certain problème linéarisé en .La formulation de ce problème linéarisé repose sur les résultats de simulation ainsi que sur leurs dérivées par rapport à la paramétrisation considérée.

Soient :

- Recherche linéaire : le pas t(k) est calculé afin de satisfaire la relation de descente :

Diverses méthodes d'optimisation sont utilisées en industrie pétrolière. Nous citerons à titre d'exemple la méthode de plus profonde descente, la méthode de Fletcher-Powell, la méthode de Levenberg-Marquardt et la méthode de GaussNewton, toutes méthodes qui sont bien connues des gens de l'art.

La mise à jour d'un modèle géologique par les données dynamiques, repose sur la combinaison des différentes méthodes et techniques qui ont été rappelées cidessus. Lorsque le modèle géostatistique est de taille raisonnable, l'inversion peut être effectuée directement sur ce dernier sans recours aux techniques de changement d'échelle. Dans ce contexte la mise àjour s'effectue selon le schéma de la figure 2.

Cependant, lorsque la taille du modèle géostatistique devient trop grande pour être utilisé directement dans le simulateur d'écoulement, l'utilisation d'une technique de changement d'échelle devient obligatoire. L'objectif de cette phase de changement d'échelle est d'effectuer les simulations d'écoulement sur un modèle de simulation de taille réduite (dit modèle grossier), permettant ainsi l'obtention des résultats de simulation dans des temps raisonnables. Dans la pratique courante, le calage s'effectue sur le modèle de simulation grossier et non pas sur le modèle géostatistique. Le principe général de cette inversion est illustré sur la figure 3.

Malheureusement, à la convergence du processus d'optimisation, seul le modèle de simulation est modifié et il est très difficile de remonter au modèle géostatistique fin sous-jacent. En effet, durant le processus d'inversion, la cohérence entre le modèle géologique initial et le modèle de simulation n'est pas maintenue.

Pour pallier ce problème, des techniques de rétrécissement d'échelle dites de downscaling ont été élaborées. Leur objectif consiste à déterminer un modèle géologique compatible avec le modèle de simulation contraint.

Ces techniques de downscaling sont assez lourdes d'un point de vue numérique, en particulier lorsque le modèle géologique est de taille assez grande.

Elles ne permettent pas toujours de remonter depuis l'échelle de simulation jusqu'à l'échelle géologique tout en respectant les contraintes géostatistiques De plus, ces techniques présentent l'inconvénient majeur de ne pas garantir que le modèle géostatistique fin issu de ce rétrécissement d'échelle permette à son tour le respect des données dynamiques (via une simulation d'écoulement sur ce modèle fin ou sur un modèle de simulation après sa mise à l'échelle).

Les nombreuses publications traitant des problèmes de calage des grands modèles géologiques notamment la publication de Wen et al. (1997) déjà citée, mettent en évidence le besoin d'une nouvelle méthodologie pour une mise à jour directe du modèle géologique fin.

La méthode selon l'invention
La méthode de mettre à jour par les données dynamiques de production un modèle géologique fin représentatif de la répartition dans le gisement d'une grandeur physique caractéristique de la structure du sous-sol (la perméabilité ou la porosité des roches du gisement par exemple).

La méthode selon l'invention permet une mise. à jour directe, par des données dynamiques, d'un modèle géologique discrétisé par un maillage fin représentatif de la répartition dans un gisement souterrain, d'une grandeur physique caractéristique de la structure du sous-sol : la perméabilité (k), la porosité ( ), etc. Elle comporte : - une paramétrisation du modèle géologique fin par un facteur de paramétrisation (#) pour obtenir la distribution de la grandeur physique dans ce modèle géologique ;

un changement d'échelle pour la détermination de la distribution de la dite grandeur physique dans un modèle de simulation discrétisé par un maillage grossier ; - la résolution via le modèle de simulation des équations d'écoulement des fluides, pour obtenir des données dynamiques simulées ; et - la détermination des relations analytiques liant les variations des données dynamiques simulées et les variations correspondantes du facteur de paramétrisation(# ).

Suivant un mode de mise en oeuvre, on détermine les relations analytiques liant les variations des données dynamiques simulées et les variations correspondantes du facteur de paramétrisation du modèle géologique fin, en combinant les dérivées des données dynamiques simulées par rapport au facteur de paramétrisation à l'échelle du modèle de simulation et les dérivées du facteur de paramétrisation du modèle de simulation par rapport au facteur de paramétrisation du modèle géologique fin.

De préférence, on procède préalablement au calibrage du modèle de simulation afin de réduire l'erreur induite par le changement d'échelle, en réalisant par exemple les opérations suivantes : - on choisit un modèle géologique fin a priori représentatif du modèle étudié (modèle de calibrage) ; - on détermine directement des premiers résultats de simulation compatibles avec ce modèle a priori ; - on détermine un modèle de simulation par changement d'échelle du modèle géologique fin ;

- on détermine directement des deuxièmes résultats de simulation compatibles avec le modèle de simulation formé, dépendant des paramètres de changement d'échelle (c) et des paramètres de simulation (s) ; et - on ajuste les paramètres de calibrage (c, s) liés au changement d'échelle et à la simulation pour que les résultats de simulation effectués à partir du modèle a priori et du modèle de simulation soient compatibles.

Les données dynamiques sont par exemple des données de production telles que la pression, le rapport volumétrique du gaz à l'huile (GOR) ou la fraction d'eau dans l'huile.

Suivant un mode de mise en oeuvre, la sélection du paramètre de paramétrisation est effectuée par une technique de déformations graduelles ou de points pilotes.

Suivant un mode de mise en #uvre, le changement d'échelle est effectué par une méthode analytique de type moyenne de puissance ou encore par une méthode numérique par résolution d'un problème d'écoulement local ou global.

En d'autres termes, la méthode selon l'invention est composée essentiellement de deux phases indépendantes qui peuvent être utilisées dans un processus itératif : une phase de calibrage et une phase de calage.

L'objectif objectif de la phase de calibrage est de réduire l'erreur induite par le changement d'échelle effectué pour réaliser la simulation d'écoulement. Un bon calibrage garantira la cohérence entre le modèle géologique fin (discrétisé par un maillage fin) et le modèle de simulation en terme d'écoulement. Ceci est primordial pour pouvoir reproduire le calage, déjà obtenu avec le modèle de simulation, en utilisant le modèle géologique fin sous-jacent ou un modèle de simulation plus grossier (modélisé par un maillage à plus grandes mailles) obtenu après une nouvelle mise à l'échelle. La méthode de calibrage proposée est basée sur des techniques de calage d'historique. Les données à caler ne sont plus les données dynamiques

observées sur le terrain mais plutôt les résultats d'une simulation de référence effectuée sur un modèle géologique donné, représentatif du modèle géostatistique étudié. Le calibrage s'effectue en utilisant le modèle de simulation obtenu suite à une mise à l'échelle du modèle géologique de référence.

L'objectif essentiel de la phase de calage est de contraindre par les données dynamiques directement le modèle géologique fin et non pas le modèle de simulation. Pour cela, nous procédons à la paramétrisation directe du modèle géologique fin, Le changement d'échelle est effectué sur le modèle géologique après sa paramétrisation. Cette phase de calage fait appel par exemple au calcul des dérivées des résultats de simulation par rapport à la paramétrisation à l'échelle du modèle géologique fin. Ceci permet d'utiliser un processus d'optimisation classique pour mettre à jour directement le modèle géologique fin.

Présentation sommaire des figures
D'autres caractéristiques et avantages de la méthode selon l'invention, apparaîtront à la lecture de la description ci-après d'un exemple non limitatif de réalisation, en se référant aux dessins annexés où : - la Figure 1 montre un domaine de changement d'échelle; - la Figure 2 montre une mise àjour du modèle géostatistique par inversion directe; - la Figure 3 montre une mise à jour du modèle de simulation après changement d' échelle; - la Figure 4 la Fig. 4 montre les dérivées des résultats de simulation par rapport à la paramétrisation du modèle géostatistique; - la Figure 5 montre un schéma général de la phase de calage; - la Figure 6 montre une validation du gradient de pression par rapport à un paramètre de déformation graduelle (puits de production);

- la Figure 7 montre une validation du gradient de pression par rapport à un paramètre de déformation graduelle (puits d'observation); - la Figure 8 montre un schéma général de la phase de calibrage; - la Figure 9 montre un schéma général de la méthodologie de couplage; - la Figure 10 montre un historique des débits; - la Figure 11montre un modèle géostatistique de référence; - les Figure 12A à 12E montrent un historique des pressions observées aux puits; - la Figure 13 montre un modèle géostatistique initial ; - les Figures 14A à 14E montrent une comparaison entre les pressions de calibrage et les résultats de simulation avant calibrage; - les Figures 15A à 15E montrent une comparaison entre les pressions de calibrage et les résultats de simulation après calibrage; - la Figure 16 montre un modèle géostatistique contraint ; - la Figure 17 montre un tableau comparatif des temps de calcul qui illustre l'intérêt de la méthode.

- - les Figure 18A à 18E montrent une comparaison entre les données de pression et les résultats de la simulation initiale ; - les Figures 19A à 19E montrent une comparaison entre les données de pression et les résultats de simulation après calage; et - les Figures 20A à 20C montrent des histogrammes respectivement des modèles géostatistiques de référence, initial et contraint.

Description détaillée de la méthode
La méthode de couplage inversion et changement d'échelle selon l'invention comporte essentiellement deux phases indépendantes qui peuvent être utilisées dans un processus itératif : une phase de calage, complétée avantageusement d'une phase préalable de calibrage.

L'objectif de la phase de calage est de contraindre par les données dynamiques directement le modèle géologique fin et non pas le modèle de simulation comme on le pratiquait jusqu'alors. A cet effet, nous procédons à la paramétrisation directe du modèle géologique.

Le changement d'échelle est effectué sur le modèle géologique après sa paramétrisation (Figure 4). Cette phase de calage fait appel au calcul des dérivées des résultats de simulation par rapport à la paramétrisation à l'échelle du modèle géologique fin. Ceci permet d'utiliser un processus d'optimisation classique pour mettre àjour directement le modèle géologique fin.

Afin de mettre en évidence les points clés de la méthode, nous détaillerons en premier la phase de calage. La phase de calibrage ainsi que le schéma itératif de la méthodologie seront détaillés plus loin.

1. Phase de calage : inversion & changement d'échelle
Dans la méthode que nous proposons, le changement d'échelle est considéré comme étant une partie intégrante de la paramétrisation et non pas comme étant un outil de pré traitement pour le simulateur d'écoulement, ce qui a pour effet que la mise à jour est effectuée directement sur le modèle géologique fin et non pas sur le modèle de simulation grossier.

1.1. Paramétrisation
La nouvelle paramétrisation que nous proposons est obtenue par la composition d'une paramétrisation classique obtenue les techniques connues dites des points pilotes ou des déformations graduelles, déjà citées, avec une technique de changement d'échelle (Figure 4) :

La mise à jour directe du modèle géostatistique par les données dynamiques de production nécessite le calcul des dérivées des résultats de simulation par rapport à la paramétrisation que nous venons de présenter. Afin de pouvoir calculer ces dérivées, il faut auparavant calculer les dérivées du modèle de simulation par rapport à la paramétrisation (Figure 4) :

Ces dérivées sont obtenues par une technique de dérivation composée (Figure 4). Pour chaque paramètre ' nous avons :

1.2. Méthode des gradients - Dérivation composée
Une fois l'ensemble de ces dérivées calculées, il est possible de déduire les dérivées des résultats de simulation par rapport à la paramétrisation de l'équation 22 :

avec :

Les dérivées des résultats de simulation par rapport à la paramétrisation seront exploitées par l'algorithme d'inversion pour calculer un jeu de paramètre #opt optimal permettant un meilleur calage des données dynamique que le jeu de paramètre initial #(0)
Le schéma général de cette phase de calage est donné par la figure 5.

1.3. Mise en oeuvre informatique: contexte monophasique
L'ensemble des algorithmes permettant la mise en oeuvre de la méthodologie selon les figures 5 et 7 pour les phases de calibrage et de calage ont été développés dans le cadre d'une boucle d'inversion généralisée.

En particulier, dans un contexte monophasique, le calcul analytique des dérivées des pressions en fond de puits par rapport à la paramétrisation du modèle géostatistique a été développé et intégré dans une version de recherche, à un simulateur d'écoulement. Le logiciel est composé essentiellement de deux modules : un module statique pour la préparation du modèle de simulation et un module dynamique pour la simulation d'écoulement.

1.3.1. Module statique
Dans ce module, la paramétrisation du modèle géostatistique fin ainsi que le processus de changement d'échelle sont exécutés. Ce module fournit non seulement le modèle de simulation grossier K(#) (Equation 21) mais également ses dérivées par rapport à la paramétrisation du modèle géostatistique (Equation 22). Une illustration simple de ce module se résume comme suit :
La paramétrisation par la méthode dite des déformations graduelles (un seul paramètre) du modèle géostatistique se résume à :

Le changement d'échelle avec la méthode de moyenne de puissance déjà citée, fournit la distribution de perméabilité du modèle de simulation :

La dérivation composée de l'équation (27) fournit les dérivées des perméabilités du modèle de simulation par rapport à la paramétrisation du modèle géostatistique. Ces dérivées sont données par :

Le terme #/##k (i) résulte de la dérivation de l'équation (26):

Soit:

Les résultats des équations (27) et (30) sont primordiaux pour un calcul analytique des gradients des pressions fond par rapport au paramètre 0. Ils seront stockés pour alimenter le module dynamique.

1. 3.2 Module dynamique
Dans le cas d'un écoulement monophasique d'un fluide peu compressible, une équation discrétisée qui régi l'évolution. de la pression dans le réservoir est donnée par :

La transmissivité TIv(I) entre la maille et une maille voisine v(I) est donnée par :

Les coefficients al , [alpha]2 et a3 sont fonctions de la géométrie des mailles 1 et v (I).

L'application de la méthode des gradients à l'équation (31) (en la dérivant par rapport à #) permet de calculer les dérivées des pressions fond par rapport à ce paramètre :

La résolution de l'équation (33) nécessite, pour chaque paramètre 19, le calcul du second membre correspondant. C'est pour le calcul de ce second membre que les résultats du module statique vont être exploités :

Ce calcul est complètement explicité en exploitant l'équation (32) qui permet

d1r(v(I) dJr(vl) de calculer les termes ao et de .

1.4. Validation numérique
Avant de pouvoir exploiter les dérivées des pressions de fond dans un processus d'inversion, une validation de ce calcul a été réalisée. Le test de validation le plus simple consiste à comparer les résultats issus de la méthode des gradients, telle que nous l'avons développée, à ceux obtenus par des gradients numériques (différences centrées).

Les résultats (Figures 6 et 7) sont présentés sous la forme de courbes montrant l'évolution au cours d'un test d'interférences des dérivées des pressions dans un puits de production et un puits d'observation par rapport au paramètre 0 de déformation graduelle.

Les figures 6 et 7 permettent de valider les résultats des dérivées que nous avons développées dans le simulateur d'écoulement utilisé.

2. Phase de calibrage
L'objectif de cette phase de calibrage est de réduire l'erreur induite par le changement d'échelle lors de la simulation d'écoulement. Un bon calibrage garantira la cohérence entre le modèle géologique et le modèle de simulation en terme d'écoulement. Ceci est primordial pour pouvoir reproduire le calage, déjà obtenu avec le modèle de simulation, en utilisant le modèle géologique fin sous-jacent ou un modèle de simulation obtenu après une nouvelle mise à l'échelle.

La méthode de calibrage proposée est basée sur des techniques de calage d'historique. Les données à caler ne sont plus les données dynamiques observées sur le terrain mais plutôt les résultats d'une simulation de référence dref effectuée sur un modèle géologique donné * (dit modèle de référence), représentatif du modèle géostatistique étudié.

Le calibrage s'effectue en utilisant le modèle de simulation obtenu suite à une mise à l'échelle du modèle géologique de référence. Dans le cadre de la méthode nous avons considéré deux types de paramètres de calibrage :
Des paramètres de calibrage relatifs au changement d'échelle que nous notons c. En effet, pour une technique de changement d'échelle donnée, le modèle de simulation dépend des diverses options que nous retenons pour appliquer la technique (conditions aux limites, exposant de la moyenne de puissance, etc.).

L'ensemble de ces données numériques constitue les paramètres de calibrage, notés c, relatifs au changement d'échelle Le modèle de simulation d'écoulement sera bien entendu une fonction de ces paramètres de calibrage :

Des paramètres de calibrage relatifs à la simulation d'écoulement que nous notons s. En effet, il est tout à fait envisageable de vouloir réduire l'erreur induite

par le changement d'échelle en modifiant quelques paramètres du simulateur d'écoulement. L'objectif de cette modification n'est pas de donner une interprétation physique de l'erreur induite par le changement d'échelle, mais plutôt d'en minimiser les conséquences grâce à certains paramètres contrôlables via le simulateur d'écoulement.

En particulier, le changement d'échelle engendre une modification de l'indice de productivité numérique (IPN) des puits. L'indice de productivité étant une fonction des perméabilités :

Afin de corriger cette modification de l' IPN , il est possible dans certain cas, d'introduire de façon artificielle, un facteur dit de skin dont la valeur sera déterminée par calibrage.

Dans un contexte plus général, les résultats de simulation en calibrage sont donc fonctions des paramètres de calibrage c et S que l'ingénieur de réservoir estime devoir calibrer :

Après avoir réalisé la simulation de référence, il est possible de procéder à l'identification des paramètres de calibrage par minimisation de la fonction objectif suivante :

Les paramètres copt et sopt issus de cette optimisation seront utilisés dans la phase de calage qui succédera à la phase de calibrage. Le schéma général de cette phase de calibrage est donné par la figure 8.

L'objectif de la phase de calibrage, comme on l'a vu, est de construire un modèle de simulation le plus représentatif possible du modèle géologique sousjacent. Pour cela on cherche à déterminer une formule de changement d'échelle (c) ainsi qu'un modèle de simulation d'écoulement (s) optimaux, c'est à dire permettant de réduire l'erreur de simulation induite par le changement d'échelle. La phase de calibrage est effectuée sur un modèle géostatistique donné, dit modèle de référence kref.

A la fin de la phase de calage qui suit la phase de calibrage, il est envisageable d'effectuer une simulation de confirmation sur une partie ou la totalité du modèle géostatistique contraint afin de s'assurer de la qualité du calage. Dans la négative, les deux phases de calibrage et de calage doivent être reprises dans un processus itératifjusqu'à la satisfaction du critère de convergence. Le schéma général décrivant cette méthodologie de couplage inversion et changement d'échelle est donné dans la figure 4.

3. Validation de la méthode
Un test d'interférences synthétique a été construit afin de valider l'utilisation de la méthodologie proposée. Ce test synthétique comprend un modèle géostatistique de référence du champ de perméabilité, supposé représentatif d'un gisement réel, et des données synthétiques de pression produites suite à la simulation d'écoulement de ce modèle de référence.

Lors de l'intégration des données de pression, le modèle géostatistique de référence est supposé inconnu. L'application de la méthodologie, présentée dans la section précédente, permet de construire un modèle géostatistique assurant le calage

des données de pression et préservant les propriétés géostatistiques du modèle de référence.

3. 1. Description du cas de validation 3.1.1. Description générale
Le cas de validation est un réservoir 3D, dont l'extension horizontale est de 4020 m x 4020 m pour une épaisseur constante de 50 m, possédant les propriétés pétrophysiques homogènes suivantes :
Porosité : Q= 0.3 ;
Anisotropie horizontale : ky/kx=1;
Anisotropie verticale : kz/kx=0.1;
Compressibilité de la roche : c = 0.0001 bar-1 ;
Viscosité : p. = 1 cP.

La pression initiale du réservoir est de 200 bars (20 Mpa) pour une saturation initiale en eau de 10%. La compressibilité de l'eau est de 0.0000435 bar-1, celle de l'huile est de 0.000168 bar-1. La compressibilité totale est égale à 0. 000255 bar-1.

Un puits producteur vertical de rayon 7.85 cm et de facteur de skin nul est implanté au centre du réservoir. Son historique de débit est composé d'une période de 7 jours en production à débit constant de 600 m3/j, suivie par une période de 35 jours de remontée de pression à débit nul' (Figure 10).

Quatre puits d'observation verticaux, perforés sur la totalité du réservoir, sont disposés à équidistance autour du puits de production, perforé uniquement sur les 30 premiers mètres du réservoir.

3.1.2. Modèle géostatistique de référence
La perméabilité du réservoir est modélisée par une fonction aléatoire de distribution log-normale de moyenne 300 mD et d'écart type 300 mD. Ce modèle géostatistique est complété par la donnée d'un variogramme sphérique de longueurs de corrélation de 600 m, 300 m et 10 m selon les principaux axes d'anisotropie :
Premier axe : (1,1,0) ;
Deuxième axe : (-1,1,0).

Le troisième axe d'anisotropie est calculé de telle façon que le repère formé par l'ensemble de ces trois axes soit orthonormé.

Le maillage de simulation géostatistique est régulier et formé de :
201 mailles de 20 m selon l'axe des x ;
201 mailles de 20 m selon l'axe des y ;
50 mailles de 1 m selon l'axe des z.

Le nombre total de mailles résultant de cette discrétisation est égale à 2 020 050.

Une simulation géostatistique de référence de ce modèle a été réalisée dans la perspective de l'obtention du champ de perméabilité de référence représentatif du réservoir que nous étudions. Ce champ de perméabilité de référence ainsi que l'emplacement de l'ensemble des puits sont présentés sur la figure 11.

3.1.3. Données synthétiques de référence
Une simulation d'écoulement, sur une période de 42 jours, est conduite sur le modèle géostatistique de référence. L'historique synthétique des pressions (Figure 12) est défini à partir des résultats de cette simulation de référence par la pression de

fond au puits de production ainsi que sa dérivée par rapport au temps ; et la pression de fond aux quatre puits d'observation.

La perturbation émise au puits producteur atteint les puits d'observation avec un retard de l'ordre de 2 jours. Seule la période succédant à ce retard est retenue pour les puits d'observation. En ce qui concerne le puits de production, seule la période en remontée de pression ainsi que sa dérivée par rapport au temps ont été prises en compte dans la formulation de la fonction objectif.

Pour la phase d'inversion, nous nous intéressons uniquement à la caractérisation de la distribution des perméabilités. Le variogramme du modèle de référence sera supposé connu de façon complète (directions principales d'anisotropie et longueurs de corrélation).

Egalement les perméabilités des mailles perforées par les cinq puits sont supposées connues. Ces perméabilités serviront pour le conditionnement du modèle géostatistique fin à chaque itération du processus d'inversion.

3. 2. Modèle géostatistique initial
Le modèle géostatistique de référence est maintenant supposé inconnu. En partant d'une nouvelle réalisation (réalisation initiale), un modèle contraint va être déterminé pour respecter à la fois les propriétés géostatistiques du modèle de référence (en terme de moyenne, écart type, variogramme etc. ) et les données de pression (en terme de calage).

La simulation d'un modèle géostatistique initial (Figure 13), avec un germe aléatoire, permet de générer un champ de perméabilité initial respectant les propriétés géostatistiques du modèle de référence. Cependant, ce modèle initial ne permet pas de respecter les données de pression. Afin de contraindre le modèle géostatistique par les données de pression, nous allons appliquer la méthodologie présentée dans la section précédente à la réalisation initiale.

3. 3 Calage de l'historique des pressions
L'intégration des données de pression dans le modèle géostatistique initial a été effectuée en utilisant : la méthode de moyenne de puissance pour la phase de changement d'échelle ; la méthode des déformations graduelles pour la paramétrisation du modèle géostatistique.

Afin de rendre compte au mieux des valeurs des perméabilités aux puits, intégrées dans le modèle géostatistique via une phase de krigeage, le changement d'échelle n'a pas été pratiqué sur les mailles puits (pour chacune des couches). Le modèle de simulation grossier ainsi obtenu est formé de :
43 mailles selon l'axe des x ;
43 mailles selon l'axe des y ;
10 mailles selon l'axe des z.

Le nombre total de maille résultant est égal à 18 490, soit une réduction de plus de 99% du nombre de maille par rapport au modèle géostatistique. De cette réduction du nombre de mailles résulte une réduction du temps de simulation d'écoulement. Ainsi la simulation d'écoulement sur le maillage géostatistique est de l'ordre de 180 minutes contre 3 minutes seulement sur le maillage grossier (station SUN ULTRA 10 - 440 MHz), soit une réduction d'environ de 99%. Une quantification plus globale des gains (en temps CPU) apportés par la méthodologie sera présentée dans la section suivante.

En ce qui concerne la méthode des déformations graduelles, on a eu recours à une seule chaîne de déformation contenant six réalisations indépendantes du modèle géostatistique.

3.3.1. Phase de calibrage
Au cours de cette phase de calibrage, on cherche à réduire l'écart entre les résultats de simulation sur les maillages fin et grossier et ceci pour le modèle géostatistique initial (Figure 13). Les paramètres de calibrage qui ont été retenus sont: le facteur de skin s pour le puits producteur P1; l'exposant w de la formule de changement d'échelle.

En effet, les résultats de simulation sur le maillage grossier suite à un changement d'échelle en utilisant une moyenne harmonique (w=-1) et un facteur de skin nul au puits producteur (la valeur du facteur de skin physique) sont en léger désaccord avec ceux obtenus sur le maillage fin (Figure 14). Ce désaccord correspond à une valeur de 195 de la fonction objectif.

Au cours du processus de calibrage, il a été possible de réduire cette différence en modifiant les paramètres (s,w) de calibrage. Ainsi en optant pour une moyenne arithmétique pour la phase de changement d'échelle (w=+l) et en introduisant un facteur de skin artificiel de -0.04, les résultats de simulation sur le maillage grossier après calibrage sont en parfait accord avec ceux obtenus sur le maillage fin (Figure 15). La valeur de la fonction objectif est alors égale à 10 et ce à l'issue de 5 itérations.

Pour ce cas test, la valeur calibrée de l'exposant w de la formule de changement d'échelle (w=+1) est un résultat connu. Cependant, nous avons délibérément initialisé cet exposant à la valeur-1 afin de valider la robustesse de la phase de calibrage.

3. 3.2 Phase de calage
Une fois la phase de calibrage terminée, il est possible de procéder à l'intégration des données de pressions dans le modèle géostatistique initial. En effet les résultats de simulation obtenus sur maillage grossier en tenant compte de la phase de calibrage (w=+1, s = -0.004) restent en désaccord avec les données de pressions (Figure 17)
La fonction objectif correspondant à cette simulation initiale est de l'ordre de 112. Afin de diminuer cette fonction objectif, une chaîne de déformation graduelle incluant six réalisations indépendantes du modèle géostatistique a été construite :

Le modèle géostatistique contraint (Figure 16) a été obtenue après 21 itérations avec une fonction objectif égale à 7,5.

Ce modèle géostatistique contraint permet un très bon calage des données de pression, comme en atteste la figure 18. Ainsi, le modèle géostatistique contraint permet non seulement le respect de ces données de pression mais également le respect des propriétés géostatistique du modèle. En particulier, les histogrammes expérimentaux correspondants aux modèles de référence, initial et contraint sont en accord (Figure 19).

3.3.3 Bilan global des gains en temps CPU
Dans le cas test présenté dans la section précédente, environ 10 heures CPU ont été nécessaires pour réaliser le calage. Une mise à jour du modèle géostatistique sans changement d'échelle aurait pu s'avérer plus coûteuse puisqu'une estimation a priori du temps CPU requit par telle approche est de l'ordre de 378 heures. Ce temps CPU est partagé entre la simulation d'écoulement et le calcul des gradients des

résultats de simulation par rapport aux divers paramètres de calage et calibrage. Une description plus détaillée de la répartition des temps CPU est donnée dans le tableau 1.

On voit donc que la méthode selon l'invention permet une mise à jour directe, cohérente et rapide du modèle géostatistique. La paramétrisation directe du modèle géostatistique fin et non pas celle du modèle de simulation grossier, ainsi que le calcul des gradients des résultats de simulation par rapport à cette paramétrisation, a permis la mise à jour directe du modèle géostatistique dans un processus d'optimisation itératif. La paramétrisation du modèle géostatistique par la méthode des déformations graduelles permet de conserver, au cours de ce processus itératif, les propriétés géostatistiques globales du modèle. Ainsi, on est assuré que le modèle géostatistique contraint, issu du processus itératif, reste cohérent Le couplage des techniques de changement d'échelle et d'inversion permet de réduire considérablement le temps CPU nécessaire pour réaliser l'intégration des données de pression dans le modèle géostatistique.

Les tests de validation effectués dans un contexte monophasique, ont mis en évidence la puissance de la méthode proposée pour la mise àjour des grands modèles géologiques par les données dynamiques.

Comme grandeur physique caractéristique de la structure du sous-sol, nous avons pris jusqu'ici comme exemple la perméabilité absolue. Ceci n'est bien sûr pas limitatif. La méthode selon l'invention pour l'intégration des données dynamiques dans les grands modèles géologiques, reste cependant applicable pour d'autres grandeurs physiques, notamment la porosité du milieu envisagé.

Claims

REVENDICATIONS

1) Méthode pour une mise à jour directe, par des données dynamiques, d'un modèle géologique discrétisé par un maillage fin représentatif de la répartition dans un gisement souterrain, d'une grandeur physique caractéristique de la structure du sous-sol, caractérisée en ce qu'elle comporte : - une paramétrisation du modèle géologique fin par un facteur de paramétrisation (#) pour obtenir la distribution de la grandeur physique dans ce modèle géologique ; - un changement d'échelle pour la détermination de la distribution de la dite grandeur physique dans un modèle de simulation discrétisé par un maillage grossier ; - la résolution via le modèle de simulation des équations d'écoulement des fluides, pour obtenir des données dynamiques simulées ; et - la détermination des relations analytiques liant les variations des données dynamiques simulées et les variations correspondantes du facteur de paramétrisation( ).

2) Méthode selon la revendication 1, caractérisée en ce que la détermination des relations analytiques liant les variations des données dynamiques simulées et les variations correspondantes du facteur de paramétrisation du modèle géologique fin, est obtenue en combinant les dérivées des données dynamiques simulées par rapport au facteur de paramétrisation à l'échelle du modèle de simulation et les dérivées du facteur de paramétrisation du modèle de simulation par rapport au facteur de paramétrisation du modèle géologique fin.

3) Méthode selon l'une des revendications 1 ou 2, caractérisée en ce que l'on procède préalablement au calibrage du modèle de simulation afin de réduire l'erreur induite par le changement d'échelle.

4) Méthode selon l'une des revendications 1 ou 2, caractérisée en ce que la phase de calibrage comporte les opérations suivantes : - on choisit un modèle géologique fin a priori représentatif du modèle étudié (modèle de calibrage); - on détermine directement des premiers résultats de simulation compatibles avec ce modèle a priori; - on détermine un modèle de simulation par changement d'échelle du modèle géologique fin ; - on détermine directement des deuxièmes résultats de simulation compatibles avec le modèle de simulation formé, dépendant des paramètres de changement d'échelle (c) et des paramètres de simulation (s) ; et - on ajuste les paramètres de calibrage (c, s) liés au changement d'échelle et à la simulation pour que les résultats de simulation effectués a partir du modèle a priori et du modèle de simulation soient compatibles.

5) Méthode selon l'une des revendications précédentes, caractérisée en ce que la dite grandeur physique est la perméabilité ou la porosité des roches du gisement.

6) Méthode selon l'une des revendications précédentes, caractérisée en ce que les données dynamiques sont des données de production telles que la pression, le rapport volumétrique du gaz à l'huile (GOR) ou la fraction d'eau dans l'huile.

7) Méthode selon l'une des revendications précédentes, caractérisée en ce que la sélection du paramètre de paramétrisation, est effectuée par une technique de déformations graduelles ou de points pilotes.

8) Méthode selon l'une des revendications 1 à 5, caractérisée en ce que le changement d'échelle est effectué par une méthode analytique de type moyenne de puissance.

9) Méthode selon l'une des revendications 1 à 6, caractérisée en ce que le changement d'échelle est effectué par une méthode numérique par résolution d'un problème d'écoulement local ou global).