FR 2678316 A1 19921231 FR 9108028 A 19910626 La présente invention concerne une méthode pour faciliter la modélisation du processus de formation d'hydrocarbures par craquage dans des roches-mères d'un bassin sédimentaire, dans le but de mieux évaluer la richesse potentielle en hydrocarbures de formations souterraines, de mieux connaître leur localisation, qui permet donc de réduire beaucoup les aléas de l'exploration pétrolière. La méthode selon l'invention concerne plus particulièrement une étape de cette modélisation, celle où l'on cherche à décrire et quantifier le processus de transformation par craquage de la matière organique. La modélisation des processus de formation des hydrocarbures dans des bassins sédimentaires a fait l'objet de nombreuses publications. L'état de l'art dans ce domaine est décrit par exemple dans les publications suivantes - C.R Acad. Sci. Paris, t. 307, Série II, p. 927-934, 1988; The American Association of Petroleum Geologists (A.A.P.G) Bulletin v.74 N,3 p. 309-335; ou . Org. Geochem. vol.16 N, 1-3 p. 1-25. Pour mieux situer Le stade concernant la transformation par craquage où intervient la méthode selon l'invention dans l'ensemble des opérations faites pour modéliser le processus de création des hydrocarbures, on va rappeler succinctement comment s'articulent les étapes successives d'une telle modélisation. L'historique d'un bassin peut être déterminé en tenant compte de données disponibles sur la profondeur et l'âge des différentes limites stratigraphiques, la durée et l'amplitude d'érosions éventuelles, la lithologie de chaque couche sédimentaire, sur des courbes donnant les variations connues ou estimées de la porosité en fonction de la profondeur qui influe sur L'épaisseur des formations et sur des données paléobathymétriques caractérisant l'affaissement du bassin étudié. Toutes ces données permettent de caractériser les modifications au cours du temps de la géométrie d'une colonne sédimentaire, l'historique de L'enfouissement des roches-mères, l'évolution de la porosité et les courbes d'enfouissement. L'étape suivante de caractérisation d'un bassin sédimentaire porte sur l'étude des transferts thermiques. On exploite à cette fin, des données obtenues dans l'étape précédente ainsi que des donnees complémentaires et l'on exprime Le transfert thermique total comme la résultante de deux effets. Le premier est le flux de chaleur par unité de surface produit du gradient géothermique (augmentation de la température avec la profondeur) par la conductivité thermique. Le deuxième terme exprime la quantité de chaleur générée ou emmagasinée dans un volume donné de matériaux et qui dépend de leur capacité thermique. La chaleur créée peut résulter d'une activité radioactive. Elle peut être aussi perdue ou emmagasinée du fait du changement rapide de la température des roches.Par ces calculs thermiques, on arrive à déterminer le profil de température te long de verticales et ses variations durant le processus de formation du bassin, l'historique des paléo-températures pour chaque roche-mère et aussi le profil des flux de chaleur. L'étape suivante de caracterisation d'un bassin sédimentaire où la méthode selon l'invention trouve son application, est l'établissement d'un modèle cinétique rendant compte de La formation des hydrocarbures par craquage primaire, de leur transformation par des réactions secondaires et leur expulsion suivant un critère fondé sur le niveau de saturation en hydrocarbures des micropores dans les roches-mères. On distingue des réactions primaires et secondaires qui sont définies expérimentalement par référence à des expériences de pyrolyse d'un échantillon de roche-mère immature. Sont considérées comme primaires des réactions intervenant dans un dispositif de pyrolyse en milieu ouvert tel qu'un pyrolyseur du type Rock-Eval où l'évacuation des produits de la pyrolyse hors d'un échantillon finement broyé, est facilitée par un gaz vecteur. Elles rendent compte de l'expulsion des hydrocarbures hors d'une roche-mère au fur et à mesure de leur création. Par secondaires, on définit des réactions supplémentaires par lesquelles les hydrocarbures formés se convertissent en fractions plus légères, et que l'on reproduit expérimentalement au moyen de pyrolyses en autoclave, en milieu fermé. Ces réactions rendent compte du processus de conversion sans expulsion hors de la roche-mère. Le schéma de craquage primaire est obtenu en considérant un ensemble de n réactions parallèles suivant lesqueLles à partir d'un potentiel initial xio Ci = 1, 2... n), et d'une énergie d'activation Ei, la masse de matière organique se transforme en une combinaison de m fractions distinctes X1, X2 ... Xm, Les coefficients de cette combinaison étant aij avec j=1, 2 ..., m. La modélisation comporte la résolution d'un système d'équations qui, tant que le seuil d'expulsion n'est pas atteint, exprime les vitesses de formation dq/dt de m fractions différentes X1 à Xm d'hydrocarbures à partir de kérogène, par la relation dq/dt = C.x où : (1) C est une matrice n.m exprimant la formation des m fractions par n réactions parallèles dites primaires; les coefficients de cette matrice sont C.. = A.. exp(-Ej/RT) . a..; (2) 1) ) x, est le vecteur des potentiels résiduels dont les composantes xi (i = 1, 2 ... n) sont des solutions respectives des équations dx./dt = -A. exp (-E./RT). xi (t) (3) avec des valeurs initiales xi(O) dépendant de la teneur en carbone et du potentiel initial xio de la réaction i; Dans la pratique, on peut considérer par exemple la décomposition de la matière organique en quatre fractions C1, C2-5, C6-14 et C15 les paramètres utilisés exprimant la décomposition, sont etalonnés en fonction des données expérimentales obtenues par un pyrolyseur ouvert du type Rock-Eval. Le schéma de craquage secondaire est obtenu en considérant que parmi les fractions X1, X2,... Xm, un certain nombre X1 à Xp sont instables et se décomposent chacune en une combinaison de m fractions distinctes X1 à Xm, les coefficients de cette combinaison étant (i,j) liés entre eux pour respecter la loi de conservation de masse et de conservation de l'hydrogène. Si l'on tient compte des réactions secondaires, les vitesses dq/dt de formation des différentes fractions s'expriment par la relation modifiée dq/dt = C.x + B.q; où (4) B est une matrice m.m qui exprime le craquage secondaire et dont les coefficients sont bij = X(i,j). A'.. exp(-E'j/RT) (5a) pour i = 1 à m, j = 1 à p et i + j ; bjj = -A. . exp(-E'j/RT) (5b) pour j = 1 à p; bij = 0 pour i = 1 à met j = (p+1) à m (5c); et q est le vecteur des concentrations qj(t) de produits résiduels de la fraction i au temps t; Les paramètres du modèle de craquage secondaire, sont étalonnés en fonction de résultats de pyrolyses en autoclave, suivant une méthode d'optimisation telle que les coefficients t(i,j) respectent la conservation du carbone et de l'hydrogène. Les vitesses dq/dt de formations des différentes fractions sont encore modifiées quand on atteint le seuil d'expulsion et deviennent dq/dt = C.x + B.q - k . ff , q (6) T est une matrice diagonale dont les termes T(i,j) représentent Les indices de mobilité des différentes fractions, # est la vitesse d'expulsion et le produit #.#.q est la vitesse de variations des concentrations expulsées.Les paramètres vérifient en outre la relation où (7) V est le volume poreux de La roche ten cmS/g de matrice rocheuse initiale), S est le seuil de saturation nécessaire pour L'expulsion max et p est la masse volumique du mélange d'hydrocarbures en fonction de la température, de la pression et de la composition. Les étapes précédentes de modélisation rendant compte du processus de formation d'hydrocarbures ayant été rappelées brièvement, on va définir à présent la méthode selon l'invention. La modélisation des phénomènes que l'on cherche à réaliser, passe par la résolution d'un système d'équations différentielles suivant le schéma itératif dit de Newton où les inconnues sont la quantité et la composition des hydrocarbures non expulsés, la quantité et la composition des hydrocarbures expulsés et le taux d'expulsion, et L'on utilise une méthode de résolution itérative avec fixation d'un pas de temps. Comme on le verra de façon plus détaillée dans la suite de la description, il s'agit d'une méthode implicite en temps où les paramètres intervenant dans les systèmes d'équation tels que les coefficients des matrices A et B, sont exprimés en fonction des valeurs que les inconnus q et X prennent en fin d'itération et non en début d'itération. On emploie un schéma itératif de Newton pour déterminer les inconnues en fin d'itération car la méthode implicite employée conduit à un système non linéaire à chaque pas de temps. Dans la pratique, la gestion du pas de temps est une source de difficultés. Comme on le verra de façon plus détaillée dans la description qui va suivre, une sélection du pas de temps inappropriée peut conduire à des résultats érronés ou créer des contraintes gênantes pour l'utilisateur qui doit procéder par interpolations en cas de décalages constatés par rapport à la grille usuelle employée pour mesurer les temps géologiques (TPS) ou limiter encore ses options possibles quand, du fait d'un pas de temps imposé, il perd la possibilité de choisir par exemple la précision qu'il recherche ou de fixer une limite en temps aux calculs effectués, quand les ressources de traitement dont il dispose ont une puissance relativement limitée. La méthode selon l'invention permet d'éviter les difficultés évoquées ci-dessus dans la conduite d'opérations de modélisation des processus de transformation de la matière organique au sein de roches-mères, dans le but de déterminer La valeur de paramètres significatifs tels que par exemple la quantité et la composition d'hydrocarbures formées dans un bassin géologique.Elle comporte la réalisation d'une suite d'itérations en découpant la durée du processus en pas de temps successifs et elle est caractérisée en ce que - l'on choisit pour effectuer lesdites itérations, un pas de temps variable teL que la variation maximale (dq) des paramètres significatifs de l'évolution des roches durant ce pas de temps reste inférieure à un certain seuil préétabli (dqm), - on sélectionne la valeur de chaque seuil en fonction des potentiels pétroliers initiaux (xio) des roches-mères et d'un facteur dépendant de la précision recherchée, et - on limite chaque pas de temps incluant un échelon de temps de l'échelle des temps géologiques pour qu'il coïncide avec ledit échelon. On sélectionne par exemple les pas de temps de façon à respecter la relation suivante relation dans laquelle 3q est la variation maximale durant lesdits pas de temps, des concentrations des différentes sections d'hydrocarbures résiduelles, et 3qmax est proportionnel au potentiel pétrolier initial, à la teneur en carbone organique de la roche-mère et dépend d'un terme définissant la précision fixée. De préférence, on limite inférieurement lesdits pas de temps de façon qu'aucun d'eux ne soit inférieur à une fraction définie de la durée d'un échelon local d'une échelle usuelle de temps géologique. La méthode peut comporter par exemple l'utilisation d'un procédé de résolution itératif selon le schéma de Newton avec fixation d'une limitation du nombre d'itérations effectuées quand l'incrément de variation de la concentration d'au moins une fraction des hydrocarbures résiduels dans la roche-mère durant un desdits pas de temps, devient inférieure à une borne déterminée en utilisant un critère d'arrêt défini. La méthode comporte par exemple la détermination des solutions d'une fonction vectorielle G(y) = O où y est un vecteur dont les composantes sont les concentrations q des différentes fractions d'hydrocarbures et la vitesse d'expulsion, cette fonction ayant deux composantes G1(y) et G2(y) dont on connait la valeur maximale, ladite limitation est définie en bornant supérieurement les variations de chacune des deux composantes de ladite fonction vectorielle (G(y)) à chaque pas de temps en tenant compte de la précision recherchée ( ), la borne étant choisie comme une fraction déterminée de ladite valeur maximale. Comme on le verra plus en détail dans la description, le lien établi entre le pas de temps et la variation des paramètres dans cet intervalle et la faculté laissée à l'opérateur de choisir le degré de précision des calculs, a pour conséquence - un gain de temps important sur une bonne part de l'histoire thermique du fait que les pas de temps sont importants tant que la température des roches-mères ne dépasse pas 60 à 70 C; - un traitement des intervalles de temps isothermes de la même façon que les intervalles non isothermes; - un traitement correct des problèmes liés aux érosions; et - une adaptation plus facile de la méthode aux différents moyens de traitement possibles. D'autres caractéristiques et avantages de la méthode selon l'invention apparaîtront mieux à la lecture de la description ci-après en se référant aux dessins annexés où : - la Fig.l est un graphique illustrant la sélection du pas de temps; - la fig.2 montre un exemple schématique de courbe de variation en fonction du temps de la pression P s'exerçant sur les roches sédimentaires au cours de leur enfouissement; - la Fig.3 montre un exemple analogue de variation de leur porosité Pr; - la Fig.4 montre un autre exemple analogue de variation de leur température; - la Fig.5 montre schématiquement au cours du temps La variation du taux de cinq fractions d'hydrocarbures formés et qui restent au sein d'une roche-mère; - la Fig.6 montre schématiquement au cours du temps la variation du taux de cinq fractions d'hydrocarbures formés et qui sont expulsées de la roche; et - la Fig.7 montre dans un cas pratique, les écarts entre les concentrations résiduelles gLobales que l'on peut observer quand on choisit différentes valeurs de la précision. La méthode comporte tout d'abord la sélection d'un pas de temps relié aux variations des paramètres. Le pas de temps est sélectionné pour que la variation maximale 3q d'une itération à la suivante obéisse à la relation On diminue le pas de temps si q est trop proche de qmax et on l'augmente dans le cas contraire, comme le montre la Fig.l. Dans la zone A du graphique, le pas de temps choisi est acceptable. Dans La zone B au contraire, le pas de temps auquel on aboutit n'est pas acceptable car il correspond à un écart trop grand des paramètres, et l'on sélectionne par exemple un pas de temps moitié de pas de temps rejeté. Le pas de temps est ainsi important dans les périodes où la température des roches est encore faible (de l'ordre de 60 à 70 C), ce qui permet des gains de temps importants sur une fraction importante de l'histoire thermique précédant la formation d'hydrocarbures. Pour les intervalles de temps isothermes, on se fixe initialement un pas de temps égal au dixième de leur durée par exemple ensuite de quoi, le mode de gestion défini ci-dessus conduit rapidement à obtenir une valeur optimale. Le critère de sélection choisi permet de tenir compte convenablement des phénomènes d'érosion et aussi indirectement des variations de la vitesse d'expulsion hors de la roche. Dans le premier cas, l'expulsion plus rapide qui se produit, aboutit à faire varier plus vite les paramètres qi ce qui en raison du critère choisi, fait diminuer le pas de temps. Dans le deuxième cas, si la vitesse d'expulsion change rapidement, il se produit corrélativement une diminution rapide des quantités résiduelles qi ce qui fait diminuer aussi le pas de temps. On définit le seuil 3q max à ne pas dépasser comme une fraction du potentiel pétrolier initial, de sorte que les résultats gardent la même précision relative si la teneur en carbone organique est modifiée d'un facteur quelconque : où # (epsilon) est une constante choisie par l'utilisateur et qui définit la précision recherchée. La méthode comporte encore le raccourcissement ou troncation systématique des segments de temps incluant un échelon de l'échelle de temps TPS qui définit lthistoire thermique d'un bassin, de façon à faciliter la comparaison de résultats avec d'autres suivant la même échelle.A cette fin, on impose au pas de temps les contraintes suivantes - si le temps tj+2 marquant un nouveau pas de temps stj+1 contient un échelon TPSk de l'échelle TPS, on le raccourcit de façon que coincide avec cet échelon TPSk; - on impose cependant une limite inférieure à ce pas de temps si le temps tj+2 est en decà de l'échelon TPSk et par rapport à lui dans un intervalle de temps min = 0,1 . e (TPSk - TPSk 1) (10) Dans un tel cas, on impose que 3t = TPSk - t j+l j+l En imposant cette règle, on évite les itérations inutiles dues à des pas de temps exagérément petits. Le pas de temps initial est fixé par exemple au tiers du premier intervalle de temps de l'échelle TPS. L'intervalle de temps convient généralement pour satisfaire le critère 3q La méthode comporte également une étape pour fixer une limite à la série d'itérations du schéma de Newton successivement effectuée et ainsi parvenir plus vite aux résultats recherchés. La discrétisation en temps du système (2a, 2b, 2c) par un schéma implicite aboutit à un système de la forme G(y) = 0 (11) où y est le vecteur (q1(tj+1),.....qm(tj+1), #(tj+1)) des inconnues à la fin du pas de temps et où G est une fonction vectorielle se décomposant en deux sous-systèmes G1 et G2 : G1(y) = (Id - Bj+1 . #t + # . #t) q(tj+1) - C x - q(tj) (12a) G2(y) = V(tj+1) Smax - Id représente la matrice identité, les autres symboles ont la même signification que dans (6) et (7). On utilise le schéma de Newton qui consiste à résoudre ce système par itérations successives, c'est-à-dire à résoudre en #yk une équation où intervient la matrice jacobienne G' : où yk-1 est le vecteur des inconnues (ql,....qm, #) après la k-1ème itération de Newton et #yk est le vecteur incrément.Les inconnues à l'itération suivante s'obtiennent donc par : yk = #yk + yk-1 (14) On constate à l'expérience qu'il faut écarter la solution consistant à interrompre les itérations quand l'incrément sur chacune des inconnues devient inférieur à une certaine fraction de cette inconnue : #yik 10-4 Ce test d'arrêt ne fonctionne pas bien dans les cas où l'un des composants yjk tend vers 0, ce qui est relativement courant.Dans le mode de fractionnement où L'on range les hydrocarbures produits en cinq fractions : C15+, C6-14, C2-5, C1 et coke, par exemple, il est fréquent que la fraction lourde C15+ soit quasiment inexistante au-delà de 180 C par suite du craquage secondaire. Dans ce cas, la condition limitante suivant la relation (15) ne fonctionne pas. La méthode consiste en conséquence à imposer les tests suivants pour décider de l'arrêt des itérations successives. G1i(yk) G2(yk) Les seuils #q'max et #G2max ont été choisis en fonction de l'ordre de grandeur des termes dont on fait la différence pour obtenir les composantes de G1 et G2. Parmi les termes positifs du membre de droite de l'équation le terme prépondérant est (12a) Id.q(tj+1) ce qui amène à définir #q'max par : #q'max = #' . Corg . ou #' est un terme multiplicatif défini par l'opérateur. Le terme positif du membre de droite de 12b est V(tj+1) Smax' dont un majorant est V(t1) Smax puisque le volume poreux (V(t) décroît au cours du temps. On choisit donc : #G2max = #' . V(t1) Smax (17b) Selon la valeur retenue par l'opérateur pour #', l'exigence de convergence est plus ou moins forte, ce qui influence la précision des résultats et les bilans de conservation de la masse. On choisit une précision sur le terme ( ') d'autant plus forte que la précision # requise sur les résultats est élevée : #' = ss . # (18) Dans cette expression, ss doit être très petit devant 1 car la précision sur yk doit être bien plus petite que la variation des inconnues pendant le pas de temps, dont un majorant proche est E.Corg (relation (9)). Dans la pratique on choisit par exemple une valeur ss = 10-3, aboutissant à des bilans de conservation de la masse meilleurs que 10-4 dans les cas testés avec une rapidité de convergence qui limite le nombre d'itérations à effectuer. Les Fig.2 à 4 montrent respectivement des courbes de variation en fonction du temps (en millions d'années ma) de la pression s' exerçant sur les roches sédimentaires au fur et à mesure de la formation, de leur porosité et de leur température dans le cas d'une énergie d'activation relativement basse. La précision 6 (epsilon) imposée par l'opérateur a une incidence importante sur le nombre de pas de temps et donc d'itérations successives. Ainsi, pour une précision epsilon de 10 2, 10 3 et 10 4, les nombre de pas de temps correspondants sont respectivement de 51, 270, 2582 et 25588 dans différents essais pratiques qui ont eté réalisés. Sur les Fig. 5 à 7, on voit également en fonction du temps, des exemples respectivement de l'évolution de la proportion de cinq fractions C15+, C6-14, C2-5, C15+ et coke résiduelles, de ces mêmes fractions expulsées de la roche-mère et la vitesse d'expulsion, obtenus pour une précision epsilon de 10 3 par exemple. La précision choisie a une incidence sur l'écart maximal en valeur absolue, relevé sur les concentrations résiduelles ou expulsées, selon la valeur choisie pour la précision epsilon (10 3 10 2 et 10-1. comparé avec la solution quasi-exacte obtenue avec une précision supérieure, comme le montre la Fig. 8. Cependant on voit que cette incidence est relativement limitée dans la pratique. REVENDICATIONS 1) Méthode pour faciliter la conduite d'opérations de modélisation des processus de transformation de la matière organique au sein de roches-mères, comportant la réalisation d'une suite d'itérations en découpant la durée du processus en pas de temps successifs, dans le but de déterminer la valeur de paramètres significatifs tels que par exemple la quantité et la composition de différentes fractions d'hydrocarbures formées dans un bassin géologique, la méthode comportant la réalisation d'une suite d'itérations en découpant la duree du processus en pas de temps successifs et étant caractérisée en ce que - l'on choisit pour effectuer lesdites itérations, un pas de temps variable tel que la variation maximale (dq) des paramètres significatifs de l'évolution des roches durant ce pas de temps reste inférieure à un certain seuil préétabli (dqm), - on sélectionne la valeur de chaque seuil en fonction des potentiels pétroliers initiaux (x. ) des roches-mères et d'un facteur dépendant io de la précision recherchée, et - on limite chaque pas de temps incluant un échelon de temps de l'échelle des temps géologiques pour qu'il coincide avec ledit échelon. 2) Méthode selon la revendication 1, caractérisée en ce que l'on sélectionne les pas de temps de façon à respecter la relation suivante relation dans laquelle 3q est la variation maximale durant lesdits pas de temps, des concentrations des différentes sections d'hydrocarbures résiduelles, et bqmax est proportionnel au potentiel pétrolier initial, à la max teneur en carbone organique de la roche-mère et dépend d'un terme définissant la précision fixée. 3) Méthode selon la revendication 2, caractérisée en ce que l'on limite inférieurement lesdits pas de temps de façon qu'aucun d'eux ne soit inférieur à une fraction définie de la durée d'un échelon local (TPSk - TPSk 1) de L'échelle de temps géologique (TPS). 4) Méthode selon la revendication 1, 2 ou 3, caractérisée en ce que l'on utilise un procédé de résolution itératif selon Le schéma de Newton avec fixation d'une limitation du nombre d'itérations effectuées quand l'incrément de variation de la concentration d'au moins une fraction des hydrocarbures résiduels dans la roche-mère durant un desdits pas de temps, devient inférieure à une borne déterminée en utilisant un critère d'arrêt défini. 5) Méthode selon la revendication 4, caractérisée en ce qu'elle comporte la détermination des solutions d'une fonction vectorielle G(y) = O où y est un vecteur dont les composantes sont les concentrations q des différentes fractions d'hydrocarbures et la vitesse d'expulsion, cette fonction ayant deux composantes 61 (y) et G2(y) dont on connait la valeur maximale, ladite limitation est définie en bornant supérieurement les variations de chacune des deux composantes de ladite fonction vectorielle (G(y)) à chaque pas de temps en tenant compte de la précision recherchée (epsilon), la borne étant choisie comme une fraction déterminée de ladite valeur maximale.