FR2872584A1 - Methode pour simuler le depot sedimentaire dans un bassin respectant les epaisseurs des sequences sedimentaires - Google Patents

Abstract

- La méthode permet, par un algorithme d'inversion itératif, de prévoir la distribution spatiale de la composition lithologique de sédiments déposés dans un bassin sédimentaire durant un intervalle de temps géologique, ainsi que l'évolution temporelle du profil de dépôt pendant toute la durée du remplissage du bassin, tout en respectant exactement les épaisseurs des séquences sédimentaires mesurées par ailleurs.- Les données d'entrée se composent des cartes d'épaisseur des couches sédimentaires étudiées, de données concernant la localisation et la composition des apports sédimentaires aux frontières du bassin sédimentaire étudié, et de paramètres physiques caractérisant le transport des sédiments pendant la période considérée. Ces données sont déterminées par interprétation de campagnes sismiques, par mesures et observations. Cet ensemble de données est appliqué à une boucle d'inversion itérative initialisée par l'accommodation fournie par un modèle diffusif multi-lithologique stationnaire. Cette boucle fonctionne alors comme un algorithme de point fixe corrigeant les accommodations à l'aide d'un préconditionnement du résidu sur les épaisseurs de la séquence obtenues par un modèle instationnaire. Ce préconditionnement est obtenu en appliquant l'application tangente du modèle multi-lithologique stationnaire au résidu sur les épaisseurs. En sortie, on obtient l'accommodation à la fin de la séquence considérée, ainsi que les solutions en bathymétries et en composition de surface, en tout temps des séquences du modèle instationnaire calé sur les épaisseurs de la couche.- Applications à la localisation de gisements d'hydrocarbures, notamment.

Description

La présente invention concerne une méthode pour établir une carte

tridimensionnelle de la distribution spatiale de la composition lithologique des sédiments déposés dans un bassin sédimentaire pendant un intervalle de temps

géologique donné, ainsi que l'évolution temporelle du profil de dépôt pendant toute la durée du remplissage du bassin, tout en respectant exactement les épaisseurs des séquences sédimentaires mesurées par ailleurs.

En s'appuyant sur la figure 1, qui montre trois étapes de remplissage d'un bassin sédimentaire, on précise tout d'abord quelques termes qui seront utilisés dans la suite. Le remplissage d'un bassin sédimentaire est constitué de couches superposées (Cl, C2 et C3), appelées également séquences sédimentaires, qui correspondent chacune au dépôt des sédiments pendant un sous-intervalle de temps de la période géologique nécessaire au remplissage du bassin. Chaque couche est donc limitée par une surface inférieure (dite "base" de la séquence) et une surface supérieure (dite "toit" de la séquence). Ainsi, le toit d'une séquence coïncide avec la base de la séquence suivante. La base de la séquence initiale correspond au fond du bassin; il s'agit du "socle" du bassin. Ces surfaces évoluent au cours du temps, suite aux déformations mécaniques que subissent certains bassins sédimentaires, pendant ou après leur dépôt. Ainsi, on définit la "subsidence" (S) qui décrit la cote du socle par rapport à un repère absolu lié au globe terrestre (la subsidence varie avec le temps) et l'eustatisme (E) qui décrit les variations de la surface des océans enregistrées simultanément sur l'ensemble de la surface terrestre. On définit également "l'accommodation" (A) comme la distance verticale entre le socle (le fond du bassin) et le niveau marin, celle-ci évoluant également avec le temps. Le "profil de dépôt" représente la surface de dépôt des sédiments à tout instant: il s'agit de la surface topographique. Enfin, la "bathymétrie" (b) décrit la distance verticale entre le profil de dépôt et le niveau de la mer.

Etat de la technique Les progrès récents en géologie, qui ont donné naissance au cours de ces vingt dernières années à la stratigraphie sismique, puis à la stratigraphie génétique, ont profondément modifié la compréhension de l'histoire du remplissage en sédiments des bassins sédimentaires sur de grandes échelles de temps et d'espace, en montrant l'influence primordiale de deux paramètres principaux: l'évolution temporelle de l'accommodation et les apports sédimentaires aux frontières du bassin.

De nombreux modèles et notamment des modèles numériques déterministes, ont été élaborés afin de pouvoir appréhender les implications géométriques et lithologiques de ces nouvelles approches.

Ces modèles numériques simulent le transport et la sédimentation (ou l'érosion) des sédiments dans le bassin, en se basant sur une description plus ou moins complexe de la nature, à partir d'une estimation de l'eustatisme, de la subsidence et de l'apport sédimentaire aux frontières du bassin étudié.

Parmi ces modèles numériques, les modèles diffusifs ont prouvé leur efficacité à travers leur utilisation dans de nombreuses études réalisées notamment pour l'industrie pétrolière, dans le but de mieux, et plus facilement, cerner les zones susceptibles de receler des hydrocarbures. De tels modèles sont décrits, par exemple, dans les documents ci-après: Rivenaes, J.C., 1988, Application of a dual-lithology, depth dependent diffusion equation in stratigraphie simulation, Basin Research, 4, 133146; - dans le brevet suivant du demandeur: FR 2 744 224 (US 5 844 799); et dans la demande de brevet publiée sous le numéro: FR 2 849 211.

Ces deux dernières références portent sur des méthodes de modélisation du remplissage de bassins sédimentaires.

Ces modèles numériques déterministes sont intégrés dans une procédure de 30 calage de leurs paramètres d'entrée, dite "procédure d'inversion". Cette procédure d'inversion vise à ajuster les paramètres du modèle, afin que les résultats fournis par ce dernier s'ajustent au mieux à la réalité observée. Parmi les paramètres à ajuster, on compte l'évolution temporelle de la subsidence et du niveau marin, c'est-à-dire de l'accommodation. Le critère d'ajustement du modèle obtenu repose, entre autre, sur la capacité du modèle à reproduire la géométrie, et notamment la carte des épaisseurs des unités sédimentaires déposées.

Cette procédure d'inversion est la plupart du temps de type "essaierreur", comme décrit par exemple dans la demande précitée de brevet publié sous le numéro FR 2 849 211. Elle peut aussi être automatisée comme décrit par exemple dans le brevet FR 2.776.393 du demandeur, portant sur une méthode de modélisation stratigraphique d'un réservoir, ou bien encore, par exemple, dans la publication suivante: - T.A.Cross and M.A.Lessenger, Construction and application of a stratigraphic inverse model, Numerical experiments in stratigraphy; recent advances in stratigraphic and sedimentologic computer simulations, Special Publication - Society for Sedimentary Geology. 62; Pages 69-83. 1999.

Dans les méthodes sus-citées, la mise en oeuvre de la procédure d'inversion fait appel à des utilisations répétées du modèle direct. De plus, même si elles peiinettent d'améliorer la concordance entre la réalité observée et le modèle obtenu, ces méthodes ne fournissent pas la garantie d'obtenir un accord satisfaisant entre la réalité observée et le résultat du modèle.

Notamment, les cartes d'épaisseur des différentes séquences sédimentaires sont généralement bien contraintes en modélisation stratigraphique, lorsqu'elles proviennent de l'interprétation d'une campagne sismique. Or, l'ajustement exact des cartes d'épaisseur des séquences de dépôt du modèle à celles qui sont données par l'interprétation de la sismique n'est absolument pas garanti par la procédure d'inversion. En pratique, l'ajustement exact n'est quasiment jamais atteint. Dans le cas de la méthode de type essai-erreur, le succès de la procédure dépend entièrement du savoir-faire et de l'intuition de l'utilisateur. Il en est de même si l'on utilise une procédure automatique: celle-ci peut ne jamais aboutir, si les données proposées initialement par l'utilisateur sont trop éloignées de la solution.

En réponse à ces difficultés primordiales d'ajustement aux mesures d'épaisseur, une procédure de modélisation basée sur le principe de diffusion stationnaire est présentée dans le document F'R 03/11 194. Elle propose une approche permettant d'effectuer une simulation basée sur un principe diffusif, dans laquelle les cartes d'épaisseur des séquences sédimentaires deviennent une donnée et non un paramètre d'ajustement. Cependant, cette procédure nécessite une hypothèse selon laquelle les sédiments se déposent à une vitesse constante en toute verticale des séquences sédimentaires modélisées. C'est ce que l'on appelle le modèle de diffusion stationnaire. Cette hypothèse peut se montrer trop réductrice dans certaines situations, notamment lorsque les coefficients de diffusion caractérisant l'efficacité du transport sédimentaire deviennent trop faibles, ce qui est le cas dans les géométries de dépôt présentant des clinoformes, qui traduisent un dépôt du sédiment dès son entrée dans le milieu marin, sans remaniement ultérieur.

La méthode selon l'invention permet, à partir d'un processus diffusif instationnaire, de calculer un modèle du remplissage multi-lithologique d'un bassin, ainsi que l'évolution temporelle du profil de dépôt pendant toute la durée du remplissage de ce bassin, tout en assurant un accord rigoureux entre les cartes d'épaisseur des séquences sédimentaires et celles de la réalité observée.

La méthode selon l'invention L'invention concerne une méthode pour former un modèle de la distribution spatiale de la composition lithologique des sédiments déposés dans un bassin sédimentaire constitué de séquences sédimentaires, ainsi que l'évolution temporelle de la bathymétrie. La méthode comporte les étapes suivantes: É on estime, à partir d'interprétations de mesures issues de campagnes sismiques et de données de puits, des données d'entrée comprenant: - des cartes donnant les épaisseurs actuelles desdites séquences sédimentaires; une carte de la topographie initiale du socle dudit bassin sédimentaire - une distribution de l'apport sédimentaire dans ledit bassin comprenant au moins les descriptions suivantes: sa localisation dans l'espace, sa composition lithologique, l'évolution du volume de l'apport dans le temps au cours du remplissage dudit bassin; - des valeurs de coefficients de diffusion caractérisant le transport desdits sédiments vers ledit bassin; É on applique lesdites données d'entrée à une procédure d'inversion itérative, dans laquelle on calcule les épaisseurs de chacune desdites séquences par résolution d'équations de transport diffusif instationnaire. Dans cette procédure, on met également à jour, de façon itérative, l'accommodation en fin de chacune desdites séquences, à l'aide d'un préconditionnement du résidu entre les épaisseurs actuelles desdites séquences sédimentaires (estimées préalablement) et celles calculées.

Selon l'invention, l'accommodation en fin d'une séquence donnée peut être initialisée par l'accommodation en fin de séquence fournie par un modèle diffusif stationnaire, et le préconditionnement peut s'obtenir par application au résidu sur les épaisseurs, d'un modèle diffusif stationnaire tangent.

L'accommodation initiale d'une séquence peut être déterminée par l'accommodation finale calculée pour la séquence située immédiatement audessous 20 de la séquence considérée.

Selon l'invention, on peut considérer, dans les modèles stationnaire et stationnaire tangent, les valeurs moyennes en temps des apports sédimentaires pour une séquence.

Enfin, selon l'invention, lorsqu'on est proche de la convergence, le 25 préconditionnement peut s'obtenir par un calcul du pas de Newton du problème inverse, à l'aide, par exemple, d'un algorithme itératif de type GMRES préconditionné par un modèle stationnaire tangent.

D'autres caractéristiques et avantages de la méthode selon l'invention, 30 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 trois étapes de remplissage d'un bassin sédimentaire.

- la figure 2 illustre une coupe transversale d'un bassin sédimentaire.

- la figure 3 illustre la succession des étapes et le flux de données de la méthode selon l'invention.

- la figure 4 illustre en détail la seconde boucle itérative de la procédure d'inversion.

les figures 5A à 5D montrent les cartes de références d'un exemple synthétique.

- la figure 5A montre la carte de l'accommodation initiale, ainsi que celle de l'accommodation en fin de séquence; la figure5B représente la carte de solution de référence pour le profil de dépôt en fin de séquence; - la figure 5C représente la concentration moyenne de la couche en lithologie 1; - la figure5D, l'épaisseur de sédiments.

Les figures 6A à 6D montrent les résultats obtenus par la méthode selon l'invention, en termes d'erreur (axe vertical) et de nombres d'itérations (axe horizontal) : - la figure 6A illustre ces résultats pour l'erreur sur la bathymétrie en fin de séquence; - la figure 6B illustre ces résultats pour l'erreur sur l'accommodation en fin de séquence; - la figure 6C illustre ces résultats pour les erreurs sur la concentration moyenne; - la figure 6D illustre ces résultats pour le résidu sur l'épaisseur de sédiments.

Description détaillée de la méthode

La méthode selon l'invention permet, à partir de données d'entrée comprenant des cartes d'épaisseur des couches, d'une carte de la bathymétrie du socle au moment de son dépôt, d'une distribution détaillée de l'apport sédimentaire dans le bassin, et de paramètres physiques caractérisant le transport des sédiments, d'estimer un modèle du remplissage multi-lithologique du bassin ainsi que l'évolution temporelle du profil de dépôt (surface topographique) pendant toute la durée du remplissage du bassin. Un modèle de la distribution spatiale de la composition lithologique des sédiments comporte une estimation de l'accommodation en fin de séquence et de la concentration de surface. La méthode nécessite l'hypothèse selon laquelle l'accommodation varie linéairement avec le temps, pendant chaque intervalle de temps nécessaire au dépôt d'une couche (séquence sédimentaire). Elle permet d'obtenir un accord rigoureux entre les cartes d'épaisseur des séquences sédimentaires ainsi modélisées et celles de la réalité observée. Il est important de noter que le processus diffusif utilisé ici n'est pas limité au cadre de la diffusion stationnaire: la vitesse de sédimentation n'est pas ici supposée constante le long de toute ligne verticale pendant le dépôt d'une séquence.

La méthode peut se décomposer en 2 étapes majeures: É on estime, à partir d'interprétations de mesures de campagnes sismiques et de mesures effectuées aux puits, les paramètres suivants: - les épaisseurs actuelles des couches sédimentaires, en profondeur; - la topographie initiale du bassin; - la distribution détaillée de l'apport sédimentaire dans le bassin, tant sur le plan de sa localisation dans l'espace, de sa composition lithologique, que de l'évolution du volume de l'apport dans le temps au cours du remplissage du bassin; - la valeur de coefficients de diffusion caractérisant le transport des 25 sédiments vers le bassin.

É On applique les données d'entrée à une procédure d'inversion itérative faisant appel à des principes de modélisation avec résolution d'équations de transport diffusif à grande échelle, de façon à déterminer un modèle de remplissage multilithologique ainsi que l'évolution temporelle de la surface topographique au cours du remplissage du bassin, tout en respectant exactement les épaisseurs des séquences sédimentaires.

Chacune de ces deux étapes, l'estimation des données d'entrée de la procédure itérative d'une part, et la détermination des grandeurs recherchées par application de la procédure itérative d'autre part, va être détaillée ci après. Mais auparavant, il est nécessaire de définir, en se référant à la figure 2 qui illustre une coupe transversale d'un bassin sédimentaire, les notations utilisées par la suite: 1 indice des séquences temporelles entre les temps t1 et t1+1 pointées par l'utilisateur.

SZ projection horizontale du bassin.

x coordonnée horizontale (axe X) d'un point dans SZ.

z: coordonnée verticale (axe Z) d'un point du bassin dans un repère fixe.

z = Hm (t) : niveau marin au temps t.

z= H(x,t) : position au temps t du socle (S) = surface de référence soumise aux déplacements verticaux tectoniques.

h(x,t) épaisseur au temps t des sédiments au-dessus du socle avec h(x, to) = 0.

b(x, t) = Hm (t) h(x, t) H (x, t) : bathymétrie au temps t.

coordonnées verticales d'un point du bassin dans le repère ayant pour origine le sommet du bassin et pointant vers le bas 20 ( = z H(x,t) h(x,t) ).

L nombre de lithologies composant les sédiments et indicées par i =1,..., L.

c. (x, e, t) : composition dans la lithologie i des sédiments au point (x, ) du bassin à l'instant t.

c (x, t) composition des sédiments à l'instant t déposés en surface en cas de sédimentation ou passant en surface en cas d'érosion.

nX normale sortante à la frontière an du bassin.

a(x, t) = b(x, t) h(x, t) : accommodation supposée globalement continue en temps et linéaire par morceaux en temps sur chaque séquence, entre les temps tl et t1+1 de sorte que l'on écrit: a(x t) = t - tl Al+1 (x) + t t1+1 Al (x) . t1+1 - tl t1 - tl+l On a par définition Ao (x) = b(x, to).

Estimation des données d'entrée de la procédure itérative La détermination, par la méthode selon l'invention, d'une carte tridimensionnelle de la distribution spatiale de la composition lithologique des sédiments déposés dans un bassin sédimentaire pendant un intervalle de temps géologique donné, nécessite l'estimation préalable des données d'entrée suivantes: 1. les cartes d'épaisseur en profondeur des séquences sédimentaires étudiées, estimées par interprétation des résultats de campagnes sismiques. A chaque séquence, le cumul des valeurs de la carte d'épaisseur permet de calculer le volume total de sédiment introduit dans le bassin pendant la séquence de dépôt. Il faut noter que ces cartes d'épaisseur sont généralement considérées comme les données les mieux contraintes en modélisation stratigraphique. L'épaisseur des sédiments au temps t1 est notée hl (x) (hl (x) = h(x, t1) ), avec ho (x) = 0. Ce sont les données de calage du problème inverse.

2. une carte de la topographie du socle du bassin à l'instant initial du remplissage de ce bassin, estimée par interprétation des données géologiques obtenues par campagnes sismiques et de forage. Plus précisément, cette carte est généralement obtenue en interpolant des valeurs ponctuelles de bathymétrie connues en analysant par exemple l'habitat de la faune fossilisée ou la géométrie des terrains fossilisés (fossile de gouttes de pluie, rides de vagues fossiles, traces de l'affouillement d'organismes marins, sables déposés dans le méandre d'un fleuve, etc...), prélevés en cours de carottage ou observées par des outils de mesures descendus dans les forages. La bathymétrie du bassin à l'instant initial est notée: b (x) = b(x, t ) . Enfin l'utilisateur doit spécifier une valeur de la bathymétrie en un point quelconque xô du bassin et à tous temps t1 que l'on notera bô = b(xô,t1).

3. une distribution détaillée de l'apport sédimentaire (AS) dans le bassin, tant sur le plan de sa localisation dans l'espace, de sa composition lithologique, que de l'évolution du volume de l'apport dans le temps au cours du remplissage du bassin. On décrit la proportion du volume global de sédiments apporté, en tout point de la frontière du bassin, par les fleuves ou les eaux de ruissellement, en fonction du temps. On peut également tenir compte de sédiments directement incorporés à l'intérieur du bassin comme c'est le cas pour la croissance de récifs carbonatés. Les apports sédimentaires sont aussi caractérisés par leur composition lithologique (proportions respectives de sable, d'argile et de carbonate par exemple) également variable dans l'espace et dans le temps. Par exemple, 25% des sédiments peuvent être introduits en une zone centrale du bassin avec une composition de 100% en carbonates, et 75% des sédiments restants introduits de manière uniforme à la frontière du bassin et composés de 50% de sable et de 50% d'argile. Ces données sont estimées par ailleurs par la synthèse des résultats de campagnes sismiques, de mesures géologiques effectuées sur les carottes et/ou les mesures diagraphiques en forage, et d'observations géologiques sur affleurement. Ces données sont donc: g1 (x) : flux total de sédiments (entrant <0 ou sortant >0) défini sur la frontière an.

S1 (x) : terme source de production totale (>0) ou consommation totale (<0) de sédiments défini sur le domaine S2.

aS2é={xxaS2,g1(x)<0}: frontière entrante de sédiments.

aS2s = E aQ, g1 (x) 0}: frontière sortante de sédiments.

np={xES2,S1(x)>0} : sous domaine de S2 de production de sédiments.

52 = {x E S2, S1(x) >_ 0} : sous domaine de 52 de consommation de sédiments.

L

u1 e (x) 0 (avec 1,ui e (x) =1) : fraction de flux entrant en lithologie i définie i=1 sur la portion de frontière aS2é .

L

ut p (x) 0 (avec 1/d,p p (x) =1) : fraction de production en lithologie i définie i=l sur le sous domaine SZp.

On notera S; (x,t) le terme source égal à ,uf p(x) Si(x) sur SZ p et à c; (x,t) Si(x) sur QI.

En cas d'érosion de la première couche, il faut également spécifier la composition initiale à to des sédiments érodés: c; (x, e, to) . 4. des paramètres physiques caractérisant le transport des sédiments pendant la période considérée: les coefficients de diffusion dans le milieu marin et dans le domaine continental. Plus précisément, il s'agit de deux coefficients de diffusion par lithologie: un coefficient de diffusion valable en milieu marin, un coefficient valable en milieu immergé (lacustre ou marin). Ces coefficients quantifient l'efficacité du transport des sédiments en milieux continental et marin à l'époque géologique considérée. Ces coefficients sont estimés préalablement à l'application de la méthode, soit empiriquement, soit par une inversion manuelle d'un modèle direct simplifié, comme par exemple, le modèle stationnaire décrit dans la suite. Les coefficients de diffusion de la lithologie i en milieux marin et continental sont notés: k;,,,,, k/,, . On k.

suppose que le ratio rm = gy est indépendant de la lithologie i. k

Toutes ces données d'entrée sont nécessaires à la description des modèles utilisés dans la procédure itérative.

La figure 3 illustre la succession des étapes ainsi que le flux de données de la méthode selon l'invention. Le profil topographique initial, c'est à dire la surface topographique initiale, est donné en entrée de la procédure globale, c'est à dire en amont du premier niveau d'itération, alors que les autres données (carte d'épaisseur de la couche considérée (hi+l) ainsi que les coefficients de diffusion (k; ,n, k, c) et les apports sédimentaires (AS1) de la couche considérée) sont fournies au premier niveau d'itération.

Application de la procédure itérative Cet ensemble de données est fourni à une procédure itérative d'inversion, comportant deux niveaux d'itérations. Le premier niveau d'itération consiste à traiter les couches successivement dans l'ordre chronologique de leur dépôt (figure 3) . Puis, pour chaque couche, on applique un second niveau d'itérations qui comporte à chaque pas l'appel à un modèle numérique diffusif multi- lithologique instationnaire (figure 3) tel que celui présenté dans Flt 2 744 224, et à des modèles numériques diffusif multi-lithologiques stationnaire et stationnaire tangent (figure 4) décrits dans le document FR 03/11 194. Ce second niveau d'itération a pour objectif d'ajuster l'accommodation à la fin du dépôt de la couche, de manière à respecter strictement l'épaisseur de la couche.

Avant de présenter l'ensemble de la procédure d'inversion itérative, il est nécessaire de décrire les trois modèles utilisés dans la boucle d'itérations.

On considère dans la suite une séquence donnée entre les temps t1 et t1+l sachant que l'accommodation Ai (x) est supposée connue soit par le calcul effectué à la couche précédente soit par la bathymétrie au temps initial si 1=0.

Ym b si b > 0, On définit également la fonction tg (b) = b si b0.

1- Modèle direct diffusif mufti-lithologique instationnaire Ce modèle est celui décrit dans la référence RIVENAES. Étant données g1 (x) , Si (x) , (x) , , m1, , k, c, les variables calculées à la couche précédente (ou données au temps initial) c; (x, e,,t1) , b(x,t), AI (x) et l'inconnue du problème inverse de la séquence Al+l (x) , ce modèle détermine pour tout temps de la séquence l'épaisseur h(x,t) , la bathymétrie b(x,t) la composition de surface c, (x,t) et la composition du bassin et (x, , t) solution du système a, h(x, t)ci (x,0, t) + div( kt, c, (x, t)O (b(x, t)) = St (x, t) sur S2 x

L

E (x, t) = 1 sur SZ x), i=1 L 1 kt cl (x, t)O t/ (b(x, t)).nx = g1(x) sur âSZs x (tl,tl+1), 1=1 kt cc; (x, t)V ip (b(x, t))Énx = ,ut,e (x)g1(x) sur ôS2é x (tl,tl+1), b(x,t = t1) = b(x,tl) sur Q, b(x,t) h(x,t) t tl At+1 (x) + t-t1+1 Al (x) sur SZx(tl,tl+1), t1+1 tl tl tl+l a tci (x, e, t) + a r h(x, t)a eci (x, , t) = 0 sur Q x (0,+ O) x (tl,t1+1), ci(x,0,t) = c1 (x,t) si a,h(x,t) > 0, ci (x, e, t = tl) = ci (x,;, tl).

On notera la solution du modèle instationnaire (INST) sous la forme = INST condensée: b(x,t) h(x,t) c = (x, t) ci (x, e, t)gl (x), Si (x),,ui,e (x),,ui p (x) l ri b(x,tl),ci (x, ,t1),Ai (x) 91+1(x) 2- Modèle diffusif multi-lithologique stationnaire Il s'agit de l'approximation du modèle instationnaire obtenue en considérant que la vitesse de sédimentation est indépendante de t sur la séquence et que la concentration de la couche précédente est indépendante de (en cas d'érosion de la séquence précédente). Le modèle décrit dans la suite, fournit alors directement sans boucle d'inversion une approximation de l'accommodation en fin de séquence Al+1 (x) Étant donnés gl (x) , Si (x) , ,ut e (x) , ,ut p (x) , rm, , kilo, les épaisseurs de sédiments hl (x) et h1 1 (x) , la bathymétrie bô+1 en un point xô+l du bassin au temps t1+1 et la concentration calculée à la couche précédente (ou donnée au temps initial) 2872584 14 ci (x, , tl) , le modèle détermine l'accommodation en fin de séquence Al+1(x) solution du système: 1d (x)O hd+1 (x) hl (x) div( k1 r(x)) = dl (x) + Sl (x) q1 (x) sur Q, t1+1 tl

L

E dis (x) =1 sur Q 1=1

L

E ki cdi (x)V 0(x)Énx = g1(x) sur aQs, i=1 (x)V O(x)Énx =,u e (x)g1(x) sur asze, Ar+1 (x) = W-1 (ç(x)) + h1+1 (x) sur SZ, v-1 (O(xo 1)) = bo+1 avec dis (x) sur (1, iGlil p (x) sur Q1, q1 (x) _ et dis (x) si h1+,(x) ? h1 (x), di (x) = 1 J c1 (X,, )de si h1+, (x) < h1 (x).

L h1+1 (x) h1 (x) h, (x) On notera la solution de ce système stationnaire (STAT) sous la forme condensée: Al+1 (x) = STAT

-

g1 (x),S1 (x),M'i,e (x),,u,,p (x) m1 k1 i,c ci (x, e, tl), h1 (x) bo1+1, xo1+1 h1+1 (x) hl (x) 3- Modèle diffusif multi-lithologique stationnaire tangent: La correction SA de l'accommodation finale est obtenue par résolution de l'application tangente des équations du modèle stationnaire définie au point courant 15 suivant: et (x), b1+1 (x), c s,1+1, c,+1(x)) où et (x) représente l'épaisseur de la séquence, bt+1(x) la bathymétrie en fin de séquence, c, '1+1 la concentration de surface en fin de séquence et où et+1 (x) coïncide avec concentration de surface e: '1+1(x) si et (x) > 0 et représente la concentration moyenne érodée dans la séquence si et (x) < 0.

Étant donné le point courant (et (x), bt+1(x), c;'1+1, 1+1(x , les paramètres de diffusion rm c, kl, le point xô+1, les flux sortants gt (x) sur aSZs, les termes de consommation St (x) sur QI, , et le résidu sur les épaisseurs de sédiments noté R(x) , le modèle tangent calcule la correction de l'accommodation finale SA solution du système linéaire suivant: -ki cdiv[&li (x)V v(bi+i)(x)+c"' (x)VSO(x)] = R(x) c;+ (x) e1 (x)8d, (x) + St (x) (x) sur S , i=] E [&z; (x)0 v(bi+i (x)) + c1 (x)V 8O (x).Jnx = 0 sur aSZ s i=1 k:, [8d: (x)0u(bi+,(x))+c;,i+t(x)V8O(x)Jnx =0 sur aç , SA(x) = SO(x) + R(x) sur SZ, (bi+i (x)) SA(xo+1) =0, db &ï s (x) si e (x) > 0 &l s (x) sa x E SZ t avec &li (x) = t et Sqi (x) = ' t.

0 si et (x) S 0, 0 sa x E SZ, On notera la solution du modèle stationnaire tangent (STATTGT) sous la fouine condensée suivante: ti+1 ti

L

E Sd; (x) = 0 sur Q, i=l L ktc et,i+i (x)V (bc+1(x))Énx = gi (x) sur aSZ S SA(x) = STATTGT e (x), bt+1(x), cs,t+1 (x), e1+1 (x)- t t r k m,c, i,c gt (x), St (x) 1+1 x0 R(x) La figure 3 présente le schéma de fonctionnement de la méthode selon l'invention: Description des deux niveaux d'itérations de la boucle itérative d'inversion 1- Premier niveau d'itération: itérations sur les couches Un premier niveau d'itération (PI) consiste à traiter les couches successivement dans l'ordre chronologiquede leur dépôt, c'est-à-dire du bas vers le haut. La figure 3 illustre ce niveau d'itération. Les données d'entrée pour ce niveau d'itération sont constituées de la carte d'épaisseur de la couche considérée (h1+1), les apports sédimentaires (AS1) et les coefficients de diffusion de la couche considérée 10 également.

2- Deuxième niveau d'itération: itérations sur l'accommodation Comme l'illustre les figures 3 et 4, pour chaque couche, on applique un second niveau d'itérations (DI) qui comporte, à chaque pas, l'appel à un modèle numérique diffusif multi-lithologique instationnaire (MDI) tel que celui présenté dans r'll 2 744 224, et à des modèles numériques diffusifs multi-lithologiques stationnaire et stationnaire tangent présentés précédemment. Ce second niveau d'itération est illustré par la figure 3 et de façon détaillée à l'itération j par la figure 4. Il a pour objectif d'ajuster l'accommodation à la fin du dépôt de la couche (A1+1(x) ), de manière à respecter strictement l'épaisseur de la couche (h1+1). Cette procédure repose principalement sur les hypothèses suivantes, pour une séquence considérée (ter niveau d'itération) : 1. Les coefficients de diffusion sont constants en temps pendant la séquence considérée.

2. L'accommodation varie linéairement avec le temps pendant le dépôt de la séquence étudiée. L'accommodation initiale est détenninée entièrement par l'accommodation finale calculée pour la couche située immédiatement au-dessous de la couche considérée. Si la couche considérée est la couche basale, l'accommodation initiale est donnée par le profil de dépôt à l'instant initial qui est une donnée d'entrée de la procédure. L'accommodation finale est ajustée par la procédure itérative relative à la couche considérée.

Les données d'entrée pour ce niveau d'itération sont constituées de g1 (x) , Si (x) , /4,e (x) , u p (x) , rm, , k!, , les variables calculées à la couche précédente (ou données au temps initial) et (x, e,t1) , b(x, t1) et Ai (x). Dans le sens direct, le modèle instationnaire suppose connu en plus des données précédemment définies, 5 l'accommodation à la fin de chaque séquence Al+1(x) . C'est cette accommodation qui est inversée à chaque séquence par la boucle d'inversion. Al+1(x) est donc l'inconnue du problème inverse de la séquence. La boucle d'inversion a donc pour but de résoudre le problème de calage sur les épaisseurs données par le modèle instationnaire: trouver l'accommodation en fin de séquence Al+1 (x) telle que h(x, t1+1) = h1+1(x), b(x,t) g1 (x), Si (x),,u e (x),,ti, p (x) h( 1 1 x, t) = INST rm'' ki,c c,S (x, t) b(x, t1),ci (x, e, t1),Ai (x) _c1 (x, , t)_ Al+1(x) La boucle d'inversion est initialisée par l'accommodation A 1(x) fournie par le modèle diffusif multi-lithologique stationnaire (MDS, figure 3). Les données d'entrées décrites précédemment permettent en effet de calculer directement sans boucle d'inversion, via le modèle diffusif multi-lithologique stationnaire, une approximation de l'accommodation en fin de séquence A 1(x) , qui sera ensuite mise à jour dans la boucle itérative:

-

g1 (x), S1 (x),,ui,e (x), ii, p (x) 1 rm,c, i,c ci (x, e, t1), h1 (x) 1+1 xo1+1 bo, h1+1(x) h1 (x) Cette boucle fonctionne, après initialisation, comme un algorithme de point fixe corrigeant les accommodations successivement calculées (A1+1(x)) par un 20 préconditionnement du résidu (R) sur les épaisseurs de la séquence (h1+1) jusqu'à obtention d'un résidu suffisamment faible (inférieur à e donné) sur les épaisseurs.

1+1 1+1 b(xo,t1+1) bo A 1(x) = Al+1 (x) = STAT Le préconditionnement s'obtient en appliquant l'application tangente du modèle multilithologique stationnaire (MDST) au résidu sur les épaisseurs.

En sortie de la boucle d'inversion, on obtient l'accommodation à la fin de la séquence considérée Al+1 (x) , ainsi que les solutions en bathymétries b(x,t) et en concentrations (composition de surface) cf (x,t) , en tout temps de la séquence du modèle instationnaire calé sur l'épaisseur de la couche h1+1.

L'algorithme proposé dans ce brevet est un point fixe initialisé par le modèle stationnaire et préconditionné par le modèle stationnaire tangent: Initialisation:

-

g1 (x), Si (x),,u,e (x), 1 p (x) AI+1 (x) = STAT = Al0+1(x) R = oo (R est le résidu sur les épaisseurs des séquences) est donné. e est une faible valeur acceptable pour le résidu R Boucle: TANT QUE IIRII e Calcul du résidu R et du point de calcul du modèle stationnaire tangent: g1 (x), S1 (x),, i,e (x),,u: p (x) 1 1 rm,c k1,c b(x, t1),ci (x, , t1), 41 (x) AI+1 (x) = INST R(x) = h1+1 (x) h(x, t1+1), b1+1(x) = b(x, t1+1), e1+1(x) = h(x, t1+1) h1 (x), c s,1+1(x) = c i (x, t1+1), Ci '1+1 (x) si e1+1 (x) % 0, 1 h(x,t, 1) et (x, e, t1)de si e1+1 (x) < 0,.

l+1 (x) hr (x) Préconditionnement:

-

et (x) br+1(x), c= ,r+1 (x) c? (x) r k m,c, c g1 (x), Sr (x) 1+1 x0 R(x) 5 Correction de l'accommodation: AI 1(x) F Ar+1 (x) + SA(x)

FIN TANT QUE

Le schéma numérique utilisé pour cet algorithme est le schéma en espace (schéma de discrétisation des équations sur le domaine défini par le bassin) décrit par exemple dans la référence Eymard et aI déjà citée. Il est similaire à celui introduit dans le papier de Rivanaes déjà cité également. Il s'adapte à des maillages du domaine SZ Cartésien, Corner Point (CPG), triangulaire, avec ou sans raffinement local.

Remarques a. Sans préconditionnement la correction s'écrirait simplement: Ar+1(x) Ar+1 (x) + R(x) R(x0i+1) . On peut montrer, dans le cas monolithologique (L=1) linéaire (rr, =1), que le 20 taux de convergence de l'algorithme préconditionné est au moins égal à 0.5 C1+1(x) _ 5A(x) = STATTGT quels que soient les paramètres du problème (coefficient de diffusion, domaine 52, intervalle de temps de la séquence, ...), alors que le taux de convergence de l'algorithme sans préconditionnement tend vers 1 lorsque k (t1+1 t1) tend Diamètre(n) vers zéro et n'est donc pas robuste. De plus, on montre que si l'on choisit un coefficient de diffusion du modèle tangent égal à 0.55 fois celui du modèle instationnaire pour l'étape de préconditionnement, le taux de convergence de l'algorithme préconditionné est toujours au moins égal à 0.1.

b. On peut montrer, dans le cas monolithologique (L=1) linéaire (rn, =1), que le problème inverse de calage exact des épaisseurs de la couche est bien posé au sens où le bruit haute fréquence sur les données n'est pas amplifié.

c. L'algorithme s'étend facilement dans le cas où les apports dépendraient aussi du temps au sein de la séquence. Il suffit de considérer les valeurs moyennes en temps dans la séquence dans les modèles approchés stationnaire et stationnaire tangent.

d. Pour les très petits ratios rm = k'jm (inférieurs ou égaux à 0.01), l'algorithme de l,, point fixe préconditionné par le modèle stationnaire tangent peut être amélioré en terme de robustesse en lui ajoutant une étape de recherche linéaire selon la direction SA. Précisément, on cherche, par recherche linéaire, une valeur de a E]0,1] , telle que la norme, notée f(a), du résidu sur les épaisseurs obtenu avec la mise à jour AI+, (x) + aSA(x) de l'accommodation en fin de séquence décroisse.

e. L'algorithme de point fixe, préconditionné par le modèle stationnaire tangent, peut être accéléré lorsque l'on est suffisamment proche de la convergence. Typiquement, on estime être proche de la convergence lorsque le résidu sur les épaisseurs ou la norme du pas SAa décru d'un facteur de l'ordre de 10 à 100. Cette accélération est obtenue en calculant le pas de Newton du problème inverse défini par:

- --

SANewton âA1+, - R Ce calcul peut être effectué par un algorithme itératif GMRES préconditionné par le modèle stationnaire tangent, tel que décrit dans le document suivant: - Y. SAAD, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Statist. Comput., 7 5 (1986), pp. 856-869.

On obtient ainsi une convergence quadratique du résidu des épaisseurs, lorsque l'on est assez proche de la solution: asymptotiquement, la norme du résidu sur les épaisseurs, à une itération donnée, est de l'ordre du carré de celle de l'itération précédente.

Exemple d'application Il s'agit tout d'abord de construire un exemple synthétique sur une seule couche d'un bassin 3D, permettant d'apprécier les résultats de l'invention. Le domaine S2, est un rectangle (0, Lx) x (0, Ly) avec Lx =100 km et Ly = 50 km. Les flux entrants aux frontières sont égaux à -1,5 m2/an sur le coté x=0, à -0,75 m2/an sur le coté y=50 km et les flux sortants sont nuls sur les deux autres cotés. Les sédiments sont constitués de deux lithologies (L=2) de coefficients de diffusion kt, = 500 m2/an, k2, = 50 m2/an, kl m = 50 m2/an et kz m = 5 m2/an. Les flux fractionnaires entrants en x=0 et en y=50 km sont égaux à 0,5. La composition initiale du bassin est égale à 0,5. La durée de la séquence est de 10 Ma. Le maillage de calcul est Cartésien uniforme de taille 40x20 mailles et le pas de temps est unifoime égal à 1 Ma. L'exemple synthétique est construit en considérant une bathymétrie initiale définie, en km, par bo (x) = Ao (x) =100(1 cos z(x y)) . L'épaisseur de la première couche (l = 0) à caler est obtenue par un calcul du modèle instationnaire pour une accommodation en fin de séquence égale à Al (x) = 400(1 cos r(x + y)) . Outre les données en épaisseur, on obtient la solution, qui servira de référence, pour la concentration moyenne dans la couche et pour la bathymétrie en fin de séquence. Cette solution de référence nous servira à calculer les erreurs en concentration moyenne et en bathymétrie à chaque itération de l'algorithme d'inversion.

Les figures 5A à 5D montrent les cartes de références, c'est à dire celles constituant l'exemple synthétique et que l'on cherche à estimer par la méthode selon l'invention. La figure 5A montre la carte de l'accommodation initiale (AI), ainsi que celle de l'accommodation en fin de séquence (AFS). La figure 5B représente la carte de solution de référence pour le profil de dépôt (i.e. l'opposé de la bathymétrie) en fin de séquence. La figure 5C représente la concentration moyenne de la couche en lithologie 1 et la figure 5D, l'épaisseur de sédiments.

On applique ensuite la méthode selon l'invention à ce modèle synthétique. Les épaisseurs obtenues sont utilisées comme données du problème inverse. On impose également la valeur bô de la bathymétrie obtenue en fin de séquence au point de contrôle xô du maillage. Ce point de contrôle correspond au point où la bathymétrie est minimale en valeur absolue. On peut repérer sur la figure 5B, par le trait en pointillé vertical, le point xô où la bathymétrie (ou de façon équivalente l'accommodation en fin de séquence) est imposée.

Les figures 6A à 6D montrent les résultats obtenus par la méthode selon l'invention, en termes d'erreur (axe vertical) et de nombres d'itérations (axe horizontal). Ces figures montrent les courbes de convergence de la boucle d'inversion avec (PNL1, PNLi) et sans (NPNL1, NPNLi) préconditionnement, et pour les deux normes Ll (PNL1, NPNL1) et L (PNLi, NPNLi) définies respectivement par: 1 Jlu(x)ldx et IluL =sup SZlu(x)I. mesure(Q) La figure 6A illustre ces résultats pour l'erreur sur la bathymétrie en fin de séquence, la figure 6B pour l'erreur sur l'accommodation en fin de séquence, la figure 6C pour les erreurs sur la concentration moyenne et enfin, la figure 6D illustre ces résultats pour le résidu sur l'épaisseur de sédiments.

Les deux algorithmes avec et sans préconditionnement sont initialisés avec le modèle stationnaire. On constate clairement le gain très important en vitesse de 30 convergence pour l'algorithme avec préconditionnement.

llull1 =

Claims (7)

REVENDICATIONS

1) Méthode pour former un modèle de la distribution spatiale de la composition lithologique des sédiments déposés dans un bassin sédimentaire constitué de séquences sédimentaires, ainsi que l'évolution temporelle de la bathymétrie, caractérisée en ce que la méthode comporte les étapes suivantes: É on estime, à partir d'interprétations de mesures issues de campagnes sismiques et de données de puits, des données d'entrée comprenant: - des cartes donnant les épaisseurs actuelles desdites séquences 10 sédimentaires; - une carte de la topographie initiale du socle dudit bassin sédimentaire; - une distribution de l'apport sédimentaire dans ledit bassin comprenant au moins les descriptions suivantes: sa localisation dans l'espace, sa 15 composition lithologique, l'évolution du volume de l'apport dans le temps au cours du remplissage dudit bassin; - des valeurs de coefficients de diffusion caractérisant le transport desdits sédiments vers ledit bassin; É on applique lesdites données d'entrée à une procédure d'inversion itérative, dans laquelle on calcule les épaisseurs de chacune desdites séquences par résolution d'équations de transport diffusif instationnaire, et on met à jour de façon itérative l'accommodation en fin de chacune desdites séquences, à l'aide d'un préconditionnement du résidu entre les épaisseurs actuelles desdites séquences sédimentaires et celles calculées.

2) Méthode selon la revendication 1, dans laquelle l'accommodation en fin d'une séquence donnée est initialisée par l'accommodation en fin de séquence fournie par un modèle diffusif stationnaire.

3) Méthode selon l'une des revendications précédentes, dans laquelle le préconditionnement s'obtient par application au résidu sur les épaisseurs, d'un 30 modèle diffusif stationnaire tangent.

4) Méthode selon l'une des revendications précédentes, dans laquelle l'accommodation initiale d'une séquence est déterminée par l'accommodation finale calculée pour la séquence située immédiatement au- dessous de la séquence considérée.

5) Méthode selon l'une des revendications précédentes, dans laquelle on considère, dans les modèles stationnaire et stationnaire tangent, les valeurs moyennes en temps des apports sédimentaires pour une séquence.

6) Méthode selon l'une des revendications précédentes, dans laquelle, lorsqu'on est proche de la convergence, la mise à jour de l'accommodation en fin de chacune desdites séquences s'obtient par un calcul du pas de Newton du problème inverse.

7) Méthode selon la revendication 6, dans laquelle on calcule le pas de Newton du problème inverse par un algorithme itératif de type GMRES préconditionné par un modèle stationnaire tangent.