Justin Herbert

Merci Charline

Denis-Jo Audrey

P Durand-Riard

Gestion

J H Shaw

A Plesch

éléments finis . . . . . .

 3.8Volumetric mismatch between both restorations at restoration step 1 . . . . . 3.9 Volumetric mismatch between both restorations at restoration step 2 . . . . . 3.10 Volumetric mismatch between both restorations at restoration step 3 . . . . . 3.11 Scheme of bad property transfer between two unrestored models . . . . . . . . 3.12 Example of bad property transfers in the analog model . . . . . . . . . . . . . 3.13 Norm of the difference between restoration displacements on faults (step 1) . . 3.14 Norm of the difference between restoration displacements on faults (step 2) . . 3.15 Norm of the difference between restoration displacements on faults (step 3) . . pour ton encadrement quand je n'étais qu'un petit master, et pour l'aide et le soutien que tu as pu m'apporter pendant ma thèse. Mais surtout merci pour les moments à écouter Waka Waka (au grand désespoir de Fifi). Fifi les patates c'est lourd ! Merci pour toute l'aide que tu as pu m'apporter pendant la thèse et avant quand j'étais un petit 2A/3A. Désolé si mes goûts musicaux t'ont cassé les oreilles ou si la vue de mon petit chien t'a donné envie de vomir... Ta vitalité est très appréciée au labo. Attention à l'orthographe, toi qui travailles sur les réseaux de factures... euh de fractures. Ce fut un plaisir de passer des soirées à jouer à des jeux de société autour d'un verre avec toi, Julie et vos deux beaux enfants. Force, honneur, respect, robustesse ! Marion, ma directrice de thèse qui fait aussi de la restauration, ce fut un plaisir d'être dans le même bureau. Merci pour ton aide. Bravo de m'avoir supporté, surtout lors de la fin de ma thèse (rédaction, préparation de la soutenance). Un grand merci d'avoir relu et corrigé (plusieurs fois) mon manuscrit. Évite de trop roupiller pendant ton temps de travail :) ! J'espère que mes bons conseils, que je te donne depuis ta 2A, ne vont pas te manquer (I am useful !). Au fait, il faut que tu programmes des équations :). Et faudra que tu me fasses d'autres muffins aux fruits rouges (c'est une tuerie) ! Jon ! On a partagé tellement de choses ensemble : Kaamelott, la paléontologie, la vodka zoolou... Ce fut cool de passer ces 3 ans (euh non, presque 4 !) dans le même labo. Elle est où la poulette ? Un jour on ira voir Guillaume pour qu'il nous reconnaisse en tant que tels ! On en a gros ! C'est super Kaamelott ! Désolé si j'ai jeté les cassettes vidéos de Jurassic Park (ça sert à quoi les VHS quand on a des DVD en 2017 ?). Gaby (mon petit Master !), merci pour tes blagues si originales ! Bravo de m'avoir supporté en tant qu'encadrant. Heureusement que tu as été là avec Faustine lors de la thèse d'Arnaud ! Antoine La Maz ! My FEM's brother :) ! On a été assez complémentaires : moi j'étais davantage sur la théorie et la biblio, et toi sur l'implémentation et le test des différentes librairies numériques. Chapeau pour tout ton travail sur RINGMecha. Mais bon qui utilise RINGMecha en 2017 ? C'est quand qu'on se fait une zienkiewicz ? Team Blandan ! Merci à toi et à Mathilde pour les moments passés ensemble, notamment au cinéma. Anthony ! Tu sais que je suis allergique aux genêts ? :) Tu es un aventurier géologue. Balaise sur le terrain (ça m'arrangeait puisque ce n'est pas mon cas...). Tu fus un super pote pendant l'école d'ingé et puis après. Tu ne disais pas que tu ne ferais pas de thèse ? Ça va me manquer les petits restos, les cinés, les parties de wii... Attention avec l'électricité ! Pierre A. (Président de la République), celui qui veut réparer les modèles. Si un jour ça marche, non quand ça marchera, ce sera super ! Tu ne m'en veux pas si je t'appelle Emmanuel ? Tu viendras faire un marathon à Boston ? Margaux ! On va manger japonais ? La serveuse était inquiète en pensant que tu étais ma copine... c'était drôle. Faut vraiment qu'on se motive à aller à la salle de sport ! Attention à ta jambe :). Jeanne, merci pour ton aide, en particulier quand j'étais en master. Avec toi c'est toujours précis, carré et efficace. J'ai beaucoup appris à tes côtés, notamment en programmation. Guillaume R. le mec qui a révolutionné la modélisation des chenaux avec des arbres. Super thèse ! Tu es quelqu'un de très efficace. Merci de m'avoir aidé pour mon déménagement. À bientôt à Boston ! Modeste qui porte bien son nom. Tu es très efficace et tu feras 3 thèses en 1. C'est très impressionnant tout ce que tu arrives à faire et c'est dommage qu'on ait pas collaboré davantage sur un projet. C'est quand que tu viens dans RINGMesh/RINGMecha ? Nicolas C., un grand merci d'avoir racheté une bonne partie de mes meubles... quand ça fait plaisir et que ça débarasse comme disait Pierre. N'oublie pas d'arroser les plantes. Julien, tu as le futur de la restauration entre tes mains. Avec toi tout est less ! J'ai hâte de voir la fin de tes travaux sur le meshless. Anne-Laure ! Yo ! On va manger chinois ? Ou japonais ? La réponse est toujours oui ! Merci pour ta bonne humeur au quotidien. Ce fut un plaisir que tu sois au labo plusieurs mois. S'il te plaît, ne me détruis pas la main :(. Bonne continuation en thèse. Paul ou Paulot ! Master of FEM. Un grand merci pour ton aide dans les équations. Je dois t'avouer que depuis que je te connais j'ai souvent une chanson de Stromae en tête : "Paulot aime les moules frites. Sans frite et sans mayo. Allez Paulot !" Laurent d'Australie qui n'aime pas mes cocottes cassoulet ! Tes "ça gazouille" vont me manquer. Jérémy, un autre fan de Kaamelott ! Et grand gamer. Il est malheureusexiv ment loin le temps où on était dans le même bureau. Tu me parlais des NURBS et des livres fantastiques que tu lisais. Pablo : cómo ça va ? Tu n'as jamais réussi à prononcer mon nom correctement :) ! Ce qui est bien avec toi c'est que tu as toujours la banane ! Pauline, la pro des karst ! Je suis toujours impressionné que tu aies fait eau puis géologie numérique. Tu as un beau parcours. Tu es toujours appliquée dans ce que tu fais et engagée dans l'enseignement. Tu apportes beaucoup de joie autour de toi. Il ne faut pas que tu aies peur des équations, elles sont tes amies ! Le labo et l'école ont de la chance de t'avoir. Dommage que je ne me souvienne pas que tu fus mon ouvreuse à la séance de Titanic. Gautier ou Gogo. Smart guy. Toujours prêt à aider. Tu m'as appris avec Fifi les bases de la programmation dans Gocad et Linux. Merci pour ton aide. Romain mon encadrant en 2A. J'ai beaucoup appris, surtout en Qt. Ce fut cool de recollaborer avec RINGMesh. Maxence, expert en éléments finis et membre de la team Blandan, ce fut chouette de parler programmation et de méthodes numériques autour d'un (plusieurs) verre(s). Si seulement tu avais été dans la team RINGMecha, on aurait pu faire de grandes choses ! Tony ! C'est super gentil d'être venu à ma soutenance ! Dommage que tu ne sois pas resté très longtemps au labo. C'est quand qu'on fait une partie de tarot ? Sophie R. : Kasteel ? Rouge ! Tu as réussi le défi d'être à la hauteur de Fatima. Je ne dirai pas que tu étais mieux ou moins bien pour éviter de me mettre une de vous deux à dos. En tout cas, merci pour ta bonne humeur quotidienne, les rigolades et les soirées au Milton... et bien sûr pour tout le travail que tu as fait. Fatima, la super secrétaire ! Toujours pas si vieille que ça :). Merci beaucoup pour ton aide et tout ce que tu fais pour le labo. Tu es toujours efficace et agréable. J'ai eu peur quand j'ai su que tu allais partir temporairement (heureusement, ça a été). Christine ! Un grand merci car sans tes cours de Python je ne me serais pas découvert une passion pour la programmation. Tu es toujours souriante et gentille. Super prof à faire des blagues du type "Les passages pythons" :). Tu apportes beaucoup de bonne humeur au labo et tu vas me manquer. Bruno, le mec le plus balaise en programmation que je connaisse. Il code plus vite que son ombre, sans EDI, sans auto-complétion, sans rien... juste Emacs et il modifie les binaires directement ! Merci pour ton aide quand j'avais des questions sur Geogram. Merci Frantz pour les discussions sur la restauration. C'est dommage qu'on ne se soit pas rencontrés plus tôt. Merci Emmanuel pour tout l'aide et le temps que tu m'as donnés quand j'ai fait le portage de RestorationLab. Je remercie Bruno Leflon et Fabien Bosquet pour le très bon stage qu'ils m'ont proposé. J'ai beaucoup appris auprès de vous. Christophe, boss en programmation. Merci pour ton aide quand j'avais des problèmes, notamment de programmation sous Windows et de portage. Tu m'as aidé au lieu de juste dire RTFM. Jean-Jacques Royer, merci pour l'aide apportée pour mon projet 2A et tous les conseils lors de ma thèse.

Vous êtes toujours là pour apporter une connaissance très bénéfique à un sujet de recherche. Pierre Jacquemin, toujours discret mais efficace. Bravo pour tout ton travail dans le projet Gocad et merci pour ton aide dans l'équipe. Jean-Laurent Mallet, merci pour votre travail en Géosciences. Vous avez apporté une contribution significative à la recherche en géologie numérique. Je me suis épanoui intellectuellement dans le laboratoire que vous avez créé.

Je remercie plusieurs autres personnes qui ne sont pas intervenues dans ma thèse mais dans ma vie et qui ont apporté leur soutien jusqu'à ce que je devienne Docteur. Rémi, Clara, Violette, Ahn, Cédric, Gaëlle, Zoya et plein d'autres de l'ENSG ! Amélie, Oriane, Katia, Cécilia, Axel, Kenneth, Pauline, Roxanne et d'autres de la prépa ! Rémy de la Science Académie ! Wilfred au lycée ! Désolé pour ceux que j'oublie. Merci à Ludovic T. pour m'avoir fait adorer la physique chimie en prépa. Un grand merci à David G. qui a amplifié mon intérêt pour les sciences pour le transformer en passion. Après tes cours, je savais que je voulais faire un métier scientifique. Merci à Paris Montagne, à la Science Académie et à Marie G. de m'avoir fait découvrir le fabuleux monde de la recherche quand j'étais au lycée. Je fais un coucou à mes amis du tir à l'arc : Georges, Daniel G., François L., Christine, Bob, Daniel et Jeanette, Bernard et Jeanine, Albert et Jeannine, Thierry, et bien d'autres. Merci au petit groupe de calligraphie chinoise. Que de bons moments passés les samedis avec Li, Fengyu, Marie-Jo, xv

Introduction

La restauration structurale correspond à un ensemble de méthodes visant à retrouver la configuration passée des structures géologiques. Dans cette thèse, nous nous concentrons aux méthodes qui annulent l'action des failles et qui déplissent les couches géologiques. Les méthodes traditionnelles sont géométriques et/ou cinématiques [e.g., Chamberlin, 1910, Dahlstrom, 1969, Gratier, 1988, Rouby et al., 2000, Massot, 2002, Groshong, 2006, Fossen, 2016]. Elles simplifient les principes fondamentaux de la mécanique et imposent un style de déformation. La restauration géomécanique, en intégrant les lois de la mécanique des milieux continus et le comportement élastique des roches, a pour objectif de s'affranchir des limites des approches classiques de restauration [e.g., Maerten et Maerten, 2001, Santi et al., 2003, Muron, 2005, Maerten et Maerten, 2006, Moretti et al., 2006, Guzofski et al., 2009, Durand-Riard, 2010, Maerten et Maerten, 2015, Al-Fahmi et al., 2016, Tang et al., 2016]. Néanmoins, plusieurs études ont remis en question l'applicabilité des approches mécaniques pour obtenir un état restauré cohérent et physique. Les conditions aux limites de la restauration géomécanique sont notamment remises en cause [Durand-Riard, 2010, Lovely et al., 2012, Durand-Riard et al., 2013b]. Par ailleurs, il existe plusieurs méthodes de restauration géomécanique. À ce jour, il n'y a pas d'inventaire exhaustif des différentes approches qui établit les avantages et inconvénients de chacune. Cela rend difficile le choix d'une méthode en particulier.

Cette thèse présente d'abord une synthèse bibliographique des différentes méthodes géomécaniques de restauration. Puis, une étude des conditions aux limites est réalisée sur un modèle analogique structural. La déformation au cours du temps de ce modèle analogique est observée sur une coupe qui sert de référence pour les restaurations. Enfin, nous comparons nos résultats de restauration sur le modèle analogique avec ceux obtenus par une méthode de restauration géométrique basée sur le modèle Géo-Chronologique (GeoChron) [Mallet, 2014, Medwedeff et al., 2016].

1 La restauration géomécanique : revue des différentes méthodes La restauration géomécanique, comme chaque méthode de restauration structurale, annule les déformations qu'ont subies les roches afin de retrouver leur état passé. Typiquement, l'horizon supérieur est remis à plat, et le rejet défini par les deux traces de cet horizon avec chaque faille est contraint à être nul. Les horizons sous-jacents suivent la déformation induite par celle imposée par l'horizon supérieur. En cas de déformations entre deux phases de dépôts sédimentaires, ou de déformations syn-sédimentaires, une restauration séquentielle est effectuée. Après la remise à plat de l'horizon supérieur, la couche supérieure est ôtée et une nouvelle restauration est faite sur l'horizon sous-jacent. Ce processus itératif se termine quand il ne reste plus de couche géologique à restaurer. La revue proposée dans cette thèse traite du déroulement d'une étape élémentaire de la restauration séquentielle.

xvii Résumé étendu

Définition de la restauration géomécanique : principes et équations fondamentales

La restauration géomécanique correspond à une simulation basée sur la mécanique des milieux continus, déformant un modèle géologique actuel afin de retrouver sa géométrie passée. Deux lois fondamentales en mécanique définissent cette déformation : la conservation de la masse et la conservation du moment linéaire. Ces deux lois ensembles constituent l'équation classique du mouvement des solides en mécanique [e.g., Muron, 2005, Maerten et Maerten, 2006, Belytschko et al., 2013] :

∇ • σ + ρb = ρü, (1)
avec ∇• la divergence, σ le tenseur des contraintes de Cauchy, ρ la densité des roches, ü le vecteur accélération et b les forces volumiques (gravité et potentiel amortissement). Pour des raisons de simplicité et de réversibilité, le comportement mécanique des roches utilisé dans la restauration géomécanique est élastique. Au sein du modèle géologique, différentes propriétés élastiques sont définies généralement de manière isotrope (module de Young et coefficient de Poisson ou équivalents). Pour prendre en compte les phénomènes de flexion-glissement (flexural slip), Durand-Riard [2010], Durand-Riard et al. [2013a] utilisent des propriétés élastiques isotropes transverses. Une loi linéaire ou non-linéaire est utilisée pour définir le tenseur des contraintes de Cauchy σ dans l' Équation (1). La loi linéaire correspond à la loi de Hooke applicable uniquement en petites déformations. Pour les grandes déformations, les lois nonlinéaires issues de l'hyperélasticité sont utilisées. Muron [2005] utilise la loi Néo-Hookéenne qui est l'extension de la loi de Hooke aux grandes déformations [Belytschko et al., 2013], et Moretti et al. [2006] et Moretti et Titeux [2007] utilisent la loi de Simo-Miehe.

Les conditions aux limites

Toute simulation mécanique a besoin de conditions aux limites afin de solliciter et de faire évoluer le domaine étudié jusqu'à atteindre un nouvel état d'équilibre ( Équation (1)). En restauration géomécanique, les conditions aux limites sont principalement des conditions en déplacement (conditions de Dirichlet). Les conditions en force de type gravité ou en traction (conditions de Neumann) ne semblent pas être utilisées en restauration même si elles sont mentionnées par Maerten et Maerten [2006]. Les conditions aux limites en général permettent de contraindre l'état restauré à partir de connaissances a priori des géométries passées. Typiquement, l'horizon supérieur est remis à plat et le jeu des failles est annulé pour cet horizon. D'autres conditions peuvent être imposées pour orienter le déplacement le long d'une direction (e.g., murs fixes, raccourcissement, extension). Par ailleurs, la plupart des méthodes de résolution numérique nécessitent de bloquer un nombre suffisant de degrés de liberté pour assurer l'unicité de la solution.

Différentes méthodes de résolution d'un problème de restauration géomécanique

Plusieurs méthodes numériques sont utilisées dans le cadre de la restauration géomécanique. La plus utilisée est la méthode des éléments finis [e.g., Zienkiewicz et Taylor, 2000a, Muron, 2005, Moretti et al., 2006, Maerten et Maerten, 2006, Tang et al., 2016] qui résout l' Équation (1) sur un maillage volumique (Section 1.4). Cette méthode a l'avantage de prendre en compte les hétérogénéités mécaniques et les géométries complexes. En revanche le coût calculatoire peut être important. Les contraintes de maillages sont aussi importantes même si les approches implicites ou sans maillage tentent de les réduire [e.g., Frank et al., 2007, Durand-Riard, 2010, Durand-Riard et al., 2010, Maerten, 2014, Zehner et al., 2015]. La résolution par éléments frontières est intéressante du fait que le modèle géologique n'a pas xviii Section 1. La restauration géomécanique : revue des différentes méthodes besoin d'être maillé par des éléments volumiques (e.g., tétraèdres, hexaèdres). Ceci a pour avantage de réduire les contraintes de maillages et le coût de calcul. En revanche les hétérogénéités mécaniques sont plus difficiles à prendre en compte. La dernière méthode est celle de la masse-ressort qui peut reposer sur une représentation par frontières ou un maillage volumique [e.g., Terzopoulos et al., 1987, Provot, 1995, Bourguignon et Cani, 2000, Macaulay et al., 2015, Midland Valley, 2017a,b]. Au lieu de résoudre l' Équation (1) sur chaque élément, elle définit un système de ressorts (par exemple correspondant aux segments des éléments) et la résolution se fait sur cet ensemble d'éléments 1D au lieu d'une formulation 3D comme dans la méthode des éléments finis. Cette formulation est beaucoup plus simple à mettre en place mais est très sensible à la discrétisation du maillage qui définit l'orientation des ressorts.

Résolution de la restauration géomécanique par éléments finis

La méthode des éléments finis est une méthode numérique permettant de résoudre une équation différentielle dans un domaine donné [e.g., Zienkiewicz et Taylor, 2000a,b, Hughes, 2012, Belytschko et al., 2013, Bathe, 2014]. Le domaine en question est discrétisé en une multitude d'éléments définissant un maillage. L'équation différentielle est résolue sur chacun des éléments. Puis la solution globale est obtenue en sommant la contribution de chaque élément. Après application du principe des travaux virtuels et plusieurs opérations mathématiques sur l' Équation (1), la méthode des éléments finis aboutit au système matriciel :

M • ü + C • u + F int -F ext = 0, (2)
avec M la matrice de masse, C la matrice d'amortissement, ü le vecteur accélération, u le vecteur vitesse, F int le vecteur des forces internes et F ext le vecteur des forces externes. Les méthodes actuelles de restauration géomécanique ne prennent pas en compte le temps, seul l'état stationnaire est important. Muron [2005], Maerten et Maerten [2006] et Moretti et al. [2006] utilisent donc la forme statique de l' Équation (2) :

F int = F ext .
(3) Muron [2005] et Moretti et al. [2006] résolvent ce système d'équations de façon globale (tous les éléments sont pris en compte dans ce système en même temps). Maerten et Maerten [2006] le résolvent de façon locale en chaque noeud avec un processus itératif de type Gauss-Seidel [Golub et Van Loan, 1996]. Muron [2005] utilise une autre méthode nommée la relaxation dynamique [e.g., Papadrakakis, 1981, Underwood, 1983, Oakley et Knight, 1995a,b]. Cette méthode garde les termes temporels de l' Équation (2) et les définit avec une intégration par différences finies. Ceci permet d'aboutir à une relation de récurrence définissant le déplacement u n+1 au temps t n+1 en fonction uniquement des grandeurs aux temps précédents :

u n+1 = u n + ∆t un+ 1 2 with un+ 1 2 = 2 -c∆t 2 + c∆t un-1 2 + 2∆t 2 + c∆t M -1 • F ext,n -F int,n , (4)
avec c un coefficient d'amortissement, et ∆t le pas de temps entre les temps t n et t n+1 . Cette formulation ne nécessite pas de résolution matricielle comme l'approche statique qui est coûteuse en mémoire. Cependant elle présente des problèmes de stabilité [Oakley et Knight, 1995a, Muron, 2005, Belytschko et al., 2013].

Gestions des failles dans les approches éléments finis

En restauration la définition de conditions aux limites adaptées aux failles est fondamentale pour assurer une bonne cohérence géologique du modèle restauré. Ces conditions assurent xix qu'il n'y a ni espace ni pénétration entre deux blocs de failles. Elles permettent aussi de déjouer l'action des failles et de remettre l'horizon supérieur dans son état topologique avant déformation tectonique, i.e., continu. Dans la littérature, deux méthodes existent pour gérer les contacts de failles. La première est purement géométrique et basée sur des conditions de Dirichlet (déplacement) [Santi et al., 2003, Muron, 2005, Moretti et al., 2006]. La deuxième repose sur la théorie des contacts mécaniques [Muron, 2005, Maerten et Maerten, 2006, Wriggers et Laursen, 2006].

1.5.1 Approche géométrique des contacts de failles Santi et al. [2003], Lepage et al. [2004], Muron [2005] et Moretti et al. [2006] utilisent une approche maître/esclave géométrique pour gérer les contacts de failles. Cette méthode consiste en un ensemble de conditions de déplacement imposé sur les noeuds de la surface de faille du bloc esclave pour le mettre en contact avec la surface de faille du bloc maître. Dans cette approche le maître reste relativement fixe par rapport à l'esclave, i.e., la condition de contact n'impacte pas le maître et le mouvement de ce dernier est dicté par les autres conditions aux limites et les lois mécaniques de conservation. Si on note A s un noeud de la surface de faille esclave à projeter, A m i les noeuds de l'élément (e.g., pour un triangle, i ∈ 0; 2 ) de la surface de faille maître contenant le projeté de A s , et a i les coordonnées barycentriques du projeté de A s dans l'élément maître, le vecteur déplacement P appliqué au noeud esclave A s est :

P = n i a i A m i -A s .
(5) Muron [2005] utilise cette approche géométrique dans le cadre de la restauration géomécanique résolue par la relaxation dynamique en se basant sur les travaux de Hallquist [1998]. À chaque pas de temps, il corrige les vitesses nodales pour prendre en compte les conditions de contacts. Il définit un coefficient de partitionnenment cinématique pour partager le déplacement entre le maître et l'esclave, permettant ainsi un mouvement bilatéral.

1.5.2 Approche mécanique des contacts de failles L'autre approche pour gérer les contacts de failles est basée sur la mécanique des contacts [Wriggers et Laursen, 2006] et est utilisée par Muron [2005], Maerten et Maerten [2006] et Tang et al. [2016]. Cette méthode repose sur la troisième loi de Newton et implique donc un équilibre des contraintes de part et d'autre de la faille. En restauration, le glissement le long de la faille se fait sans friction. Ce type de contact repose aussi sur une approche maître/esclave, mais contrairement à la méthode de contact géométrique, le bloc maître n'est pas relativement fixe. Le mouvement entre le maître et l'esclave est partagé et est une conséquence de la minimisation de l'énergie. La différence entre le maître et l'esclave se fait sur la pénétrabilité. Les noeuds du bloc esclave doivent être sur le maître (ni espacement ni pénétration) mais le contraire est possible. Pour éviter ce dernier problème, Maerten et Maerten [2006] alternent le rôle maître/esclave entre les blocs.

2 Validation de nouvelles conditions aux limites dans la restauration géomécanique 3D : application sur un modèle analogique extensif

Les conditions aux limites, utilisées dans la restauration géomécanique, sont controversées voire qualifiées de non-physiques [Durand-Riard, 2010, Lovely et al., 2012, Durand-Riard et al., 2013b] sont faibles et dont la paléo-géométrie est connue sur une coupe au cours de la déformation du modèle analogique. Cette étude permet de considérer les conditions aux limites comme responsables de la qualité de l'état restauré [Gratier et al., 1991, Schultz-Ela, 1992].

Modèle analogique structural

Le modèle analogique est un ensemble de matériaux (silicone, et grains de sable et de pyrex) déformé en laboratoire [e.g., McClay, 1990]. L'expérience a été réalisée par l'IFPEN1 et C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System2 . Ce modèle est composé initialement d'une couche de silicone à sa base et d'une couche de sable au dessus. La déformation est gravitaire et extensive, induite par l'inclinaison du modèle analogique structural de 1.5 • . Au cours de la déformation, des couches de sable et des couches de pyrex sont alternativement déposées pour imiter des dépôts syn-sédimentaires. Cette expérience permet de représenter des structures similaires à celles observées en domaines extensifs salifères comme dans le Golfe du Mexique, en Angola et au Maroc [e.g., Weijermars et al., 1993, Moretti et Callot, 2012]. Une tomographie 3D a permis d'obtenir une imagerie 3D des structures du modèle analogique dans son état déformé [e.g., Colletta et al., 1991, Callot et al., 2012]. À l'aide du géomodeleur SKUA-GOCAD [Paradigm, 2015], nous avons construit un modèle structural numérique représentant les structures observées dans la tomographie (Figure 1). Ce modèle comporte 8 horizons et un réseau de failles complexe (22 failles normales dont 5 failles coupées et déplacées par des failles plus récentes, aucune faille isolée et 27 branchements de failles). Nous avons généré un maillage volumique composé de tétraèdres, à l'aide des outils dévelopés par Lévy [2015], Pellerin et al. [2015], Si [2015a,b], Botella et al. [2016], Botella [2016] et Pellerin et al. [2017], afin de pouvoir réaliser des restaurations.

Conditions aux limites classiques et non-classiques

Afin de mener à bien la restauration séquentielle, nous avons testé différentes conditions aux limites. L'objectif est de définir les conditions aux limites permettant d'obtenir des paléogéométries proches de celles fournies par les images tomographiques.

Conditions aux limites classiques

Dans nos restaurations, nous avons défini les conditions aux limites suivantes. Premièrement, l'horizon supérieur est remis à plat. Deuxièment, le mur ouest est fixe en X et Y (Figure 1). Troisièment, les murs nord et sud sont fixes en Y. Enfin, des conditions de contacts xxi Résumé étendu de faille annulent les rejets de l'horizon supérieur et empêchent les trous et les pénétrations entre les blocs de failles [Muron, 2005, Wriggers et Laursen, 2006, Maerten et Maerten, 2006].

Conditions aux limites non-classiques

Nous avons testé une condition optionnelle de raccourcissement au niveau du mur est, condition en déplacement similaire à celles proposées par Durand-Riard [2010], Lovely et al. [2012] et Durand-Riard et al. [2013b]. Afin de pouvoir restaurer le réseau de failles complexe, nous avons défini deux nouvelles conditions de contacts. La première connecte les bords internes des surfaces de failles. La deuxième connecte les composantes connexes d'une même faille, coupée et déplacée par une faille plus récente. Avant d'appliquer cette condition aux limites, nous restaurons sans cette condition et si la continuité de la faille coupée et déplacée est apparente dans l'état restauré, nous refaisons une restauration avec cette condition. Ces nouvelles conditions de contact reposent sur le même principe que celui utilisé pour les conditions de contact classiques, à savoir la mécanique du contact [Muron, 2005, Maerten et Maerten, 2006, Wriggers et Laursen, 2006].

Restauration du modèle analogique : une condition aux limites de raccourcissement nécessaire

Nous avons restauré séquentiellement le modèle analogique en utilisant la librairie RING-Mecha 3 [Chauvin et Mazuyer, 2016] qui est basée sur les travaux de Muron [2005] et Durand-Riard [2010]. Comme le sable et le pyrex sont mécaniquement équivalents [Panien et al., 2006], nous avons défini un matériau élastique homogène pour l'ensemble du modèle à restaurer, avec un module de Young de 70 GPa et un coefficient de Poisson de 0.2 [Holtzman et al., 2009]. Pour vérifier la qualité de nos restaurations, nous avons des images tomographiques (au niveau du mur nord) de l'évolution de la déformation du modèle analogique.

Nous avons restauré en appliquant pour chaque étape de la restauration séquentielle un raccourcissement issu des images tomographiques. Nous montrons, grâce aux images tomographiques, que qualitativement et quantitativement les modèles restaurés sont géométriquement acceptables. Nous avons aussi essayé de restaurer chaque étape sans imposer de raccourcissement. En revanche, nous utilisons comme modèle non restauré à chaque étape celui obtenu lors de l'étape précédente avec le raccourcissement observé sur les images tomographiques afin d'éviter de cumuler des erreurs au cours de la restauration séquentielle. Nous montrons que sans cette condition aux limites de raccourcissement, la restauration géomécanique échoue dans la capture de l'extension du modèle analogique structural.

Évaluation du raccourcissement

Nous avons montré qu'un raccourcissement est nécessaire pour restaurer proprement le modèle analogique. Nous avons estimé le raccourcissement grâce aux images tomographiques. Dans le cadre de la restauration de modèles géologiques, il est nécessaire d'utiliser des méthodes pour estimer le raccourcissement à appliquer. Nous avons essayé différentes méthodes géométriques dont la méthode de la surface transférée. Cette méthode s'applique en coupe et permet, en plus de tester la validité d'une interprétation structurale, d'estimer l'extension ou le raccourcissement qu'a subi un domaine géologique en fonction de l'aire qui a été déplacée le long d'un niveau de détachement et application sur le modèle analogique structural la profondeur de ce niveau [e.g., [START_REF] Groshong | Area balance, depth to detachment, and strain in extension[END_REF], Groshong et al., 2003, Groshong, 2006].

Les estimations obtenues par cette méthode sont semblables à celles fournies par les images tomographiques avec moins de 15% d'erreur.

Estimation du raccourcissement par une analyse de la dilatation

Comme le modèle analogique est extensif, dans le sens de la restauration le volume du modèle devrait diminuer ou demeurer constant (il n'est pas attendu que le volume augmente). Nous avons refait la restauration séquentielle en faisant varier pour chaque étape la valeur de la condition aux limites de raccourcissement. L'objectif est de déterminer quelle valeur de raccourcissement minimise le nombre de tétraèdres avec une dilatation positive (augmentation du volume). Nous montrons sur le modèle analogique que plus le raccourcissement augmente, plus le nombre de tétraèdres incohérents diminue. Par ailleurs, ce nombre atteint un plateau lorsque l'on se rapproche du raccourcissement attendu.

3 Comparaison entre la restauration géomécanique et la restauration basée GeoChron : application sur le modèle analogique structural

La restauration géomécanique a montré dans plusieurs études ses avantages sur les méthodes classiques de restauration [Maerten et Maerten, 2006, Plesch et al., 2007, Guzofski et al., 2009]. Néanmoins, comparer deux méthodes de restauration n'est pas forcément évident puisqu'en restauration, par définition, il n'y a pas de solution de référence. De plus, de nombreuses incertitudes, notamment structurales, peuvent biaiser l'analyse [e.g., Bond et al., 2007, Bond, 2015, Cherpeau et Caumon, 2015]. Le même modèle analogique structural a été restauré avec une approche géométrique fondée sur le calcul des coordonnées chronostratigraphiques [Mallet, 2014, Medwedeff et al., 2016]. La géométrie du modèle analogique étant observable sur une coupe au cours du temps, nous proposons une comparaison objective entre les deux méthodes de restauration.

Restauration géométrique basée sur le modèle Géo-Chronologique

Le modèle analogique a été restauré en utilisant une méthode de restauration [Medwedeff et al., 2016] basée sur le modèle Géo-Chronologique GeoChron [Mallet, 2014] qui établit une bijection entre chaque point du sous-sol avec son équivalent dans l'espace de dépôt. Ainsi, GeoChron correspond à la formulation mathématique de l'espace de Wheeler. Mallet [2014] construit un champ orthogonal uvt qui correspond à l'espace chronostratigraphique [e.g., de Groot et al., 2006, Monsen et al., 2007, Wu et Hale, 2015, Labrunye et Carn, 2015, Karimi et Fomel, 2015] avec u et v les deux coordonnées qui suivent la stratigraphie (plan horizontal dans l'espace de Wheeler), et t la coordonnée perpendiculaire à la stratigraphie. Pour construire ce champ, deux styles de déformation sont disponibles : la déformation par flexion-glissement (flexural slip) ou la minimisation de la déformation [Mallet, 2014]. Medwedeff et al. [2016] utilisent la minimisation de la déformation qui correspond le plus au style de déformation du modèle analogique. Plusieurs contraintes peuvent être définies lors du calcul de l'espace uvt. Medwedeff et al. [2016] n'utilisent que l'horizon à remettre à plat comme contrainte sur le champ t, et ajoutent des contraintes pour assurer de bons contacts entre les blocs de failles. Enfin, Medwedeff et al. [2016] convertissent la dimension temporelle t de l'espace de Wheeler en une dimension métrique afin d'être cohérent avec un espace restauré. La conversion respecte les épaisseurs verticales du modèle analogique.

xxiii Résumé étendu

Comparaison géométrique entre la restauration géomécanique et la restauration GeoChron

Nous avons comparé les trois premières étapes de la restauration séquentielle réalisées avec la restauration géomécanique, basée sur les travaux de Muron [2005] et Chauvin et al. [2017], avec la restauration faite avec la méthode géométrique proposée par Medwedeff et al. [2016]. Nous montrons que, pour chaque étape de restauration, les deux modèles restaurés sont géométriquement très similaires, aussi bien qualitativement que quantitativement. Ainsi les deux méthodes de restauration fournissent des résultats équivalents. Néanmoins, tout comme la restauration géomécanique sans condition aux limites de raccourcissement (Section 2), la restauration basée GeoChron seule n'est pas capable de capturer l'extension qu'a subi le modèle analogique lors de sa déformation (le long de l'axe X). Afin de comparer la restauration géomécanique, avec des conditions de raccourcissements, avec la restauration basée GeoChron nous avons multiplié la composante X de toutes les coordonnées du solide restauré par GeoChron par un facteur afin de faire correspondre les murs ouest et est avec ceux du modèle restauré issu de la restauration géomécanique.

Impact des paramètres mécaniques dans la restauration géomécanique : cas homogène

Bien que les deux méthodes de restauration fournissent des modèles restaurés similaires géométriquement, elles ne sont pas équivalentes notamment du fait que la méthode basée sur le modèle Géo-Chronologique impose un style de déformation sans intégrer le comportement mécanique des roches. Par ailleurs, il y a des incertitudes sur le choix des paramètres élastiques dans la restauration géomécanique. C'est pourquoi nous avons réalisé plusieurs restaurations géomécaniques en variant le module de Young et le coefficient de Poisson dans le modèle analogique structural tout en restant dans un cas mécaniquement homogène.

Impact du module de Young

Afin d'étudier l'influence du module de Young nous avons testé six nouvelles restaurations avec pour coefficient de Poisson 0.2 et pour module de Young respectivement (en Pascal Pa) : 10 6 , 10 8 , 2.5 × 10 10 , 4.5 × 10 10 , 10 11 et 10 12 . Le champ de déplacement de ces restaurations est presque identique à celui du modèle restauré dans la Section 2.3. Par ailleurs, le module de Young n'a aucun impact sur le changement global et local de volume après restauration, ce qui est logique au vu de la définition du module de Young.

Impact du coefficient de Poisson

Nous avons réalisé cinq restaurations en modifiant le coefficient de Poisson (le module de Young est maintenu constant à 70 GPa). Les différents tests du coefficient de Poisson sont : 0.01, 0.1, 0.3, 0.4 et 0.45. Le coefficient de Poisson a un impact beaucoup plus significatif que le module de Young sur le champ de déplacement, même si l'écart avec le modèle restauré de référence reste de l'ordre de la limite des incertitudes structurales définie par [Chauvin et al., 2017]. En revanche, le volume global du modèle change de manière significative en fonction du choix du coefficient de Poisson. Ainsi, les coefficients de Poisson impliquant une dilatation irréaliste, d'un point de vue mécanique et géologique, peuvent être rejetés.

Impact des paramètres mécaniques dans la restauration géomécanique : cas hétérogène

La méthode de restauration basée sur GeoChron impose un style de déformation global à tout le modèle sans prendre en compte les hétérogénéités mécaniques locales, contrairement

Conclusions

Les différentes méthodes de restauration géomécanique sont toutes basées sur la conservation de la masse et la conservation des moments linéaires. Elles utilisent des propriétés élastiques pour prendre en compte la réaction des roches face à des contraintes. Les conditions aux limites sont principalement en déplacement et constituent des hypothèses sur la paléo-géométrie. La méthode numérique la plus utilisée est celle des éléments finis, elle-même subdivisée en plusieurs méthodes : statique locale, statique globale et relaxation dynamique. Ces trois méthodes sont très proches. La relaxation dynamique est moins coûteuse en mémoire et plus simple à implémenter, mais elle a des problèmes de stabilité numérique et est très lente. Les deux autres méthodes sont plus stables et très proches numériquement. L'approche masse-ressort simplifie les équations mais permet de résoudre plus rapidement un problème de restauration. L'approche par éléments frontières a quant à elle pour principal avantage de ne nécessiter qu'une représentation par frontière du modèle géologique et non d'un maillage volumique. Ceci permet de limiter les contraintes de maillage, et de réduire le coût calculatoire. En revanche les hétérogénéités sont moins bien prises en compte. Enfin, il existe deux manières de gérer les contacts de failles : géométriquement ou mécaniquement.

La restauration d'un modèle analogique obtenu en laboratoire a permis de tester l'influence des conditions aux limites. Nous avons montré qu'un déplacement latéral était nécessaire pour capturer l'extension qu'avait subi le modèle analogique structural. La méthode de la surface transférée a fourni de bonnes estimations des raccourcissements à prescrire au cours de la restauration séquentielle. Nous avons également proposé une approche empirique basée sur la dilatation pour estimer le raccourcissement. Enfin, le jeu de conditions aux limites proposé en contexte extensif pourrait être utilisé dans d'autres contextes et devrait améliorer les restaurations des cas d'études géologiques.

Nous avons montré que la restauration géomécanique et celle basée sur GeoChron, appliquées sur le modèle analogique extensif, fournissent des champs de déplacement relativement similaires. Toutefois, la méthode basée GeoChron nécessite une opération après restauration pour prendre en compte la dilatation du modèle analogique structural lors de sa déformation, ainsi que l'extension due aux failles non-représentées [e.g., Kautz et Sclater, 1988, Marrett et Allmendinger, 1992, Baxter, 1998, Groshong et al., 2003]. Cette opération est similaire à la condition de raccourcissement nécessaire dans la restauration géomécanique. Le module de Young n'a pas d'impact sur la restauration dans le cas homogène. Le coefficient de Poisson a davantage d'impact en fonction de l'amplitude de la déformation imposée par les conditions aux limites. Ce paramètre mécanique est très important pour les changements de volumes et la localisation de la déformation. Finalement les deux méthodes de restauration sont assez équivalentes dans le cadre de la restauration du modèle analogique. Néanmoins, l'approche basée GeoChron est plus simple à utiliser et très rapide. En revanche elle est peu flexible contrairement à la restauration géomécanique qui permet de définir des conditions aux limites personalisées et de jouer sur le comportement mécanique. Ces avantages de la méthode mécanique permettent une meilleure analyse des domaines géologiques complexes.

Introduction Importance of structural restoration in geomodeling

Necessity of understanding the past The Earth's crust is in constant evolution. Understanding its changes through time is fundamental in several domains. In petroleum geosciences, the economical interest of a basin is directly related to its history. Several steps are important for the formation of oil and gas. First, organic matter and sediments must be deposited in enough quantity, and be transformed in kerogen by bacteria under specific conditions [e.g., Hatcher et al., 1983, Derenne et al., 1992, Gupta et al., 2007]. Second, these sediments have to be buried. According to the burial depth and the subsequent temperature and pressure, the kerogen may produce oil and/or gas [Suggate, 1998], generating a source rock. If the gas pressure is high enough, hydrocarbons leave the source rock. For petroleum industry purposes, these fluids must be concentrated and localized to form conventional reservoirs. The presence of a permeable reservoir rock and of natural barriers (impermeable level, crystallized faults, etc.) is essential for trapping. To produce these fluids, rock permeability is also an important factor, which is function of the rock nature, its compaction and its deformation through time. All the previous steps are tracked by petroleum geologists, but a major component is the timing and the succession of these steps. If the sediments are not buried deep enough for enough time, oil/gas windows are not reached [Landais et al., 1994]. A petroleum trap must be formed prior to fluid expulsion. The knowledge of tectonic deformation through time provides valuable insight on the stress within the geological domain and therefore on the potential fracturation impacting the permeability of rocks [Macé, 2006, Maerten and Maerten, 2006, Mejía-Herrera et al., 2014, Stockmeyer et al., 2017]. In mining, the history of rock deformation, fracturing and fluid circulation is also crucial to understand and locate crystallization areas which may present economical minerals [e.g., Tripp and Vearncombe, 2004]. Finally, a pure academic understanding of a geological domain is an important topic [e.g., Muñoz, 1992, Mouthereau et al., 2001, Zhou et al., 2006]. Broadly speaking, restoration corresponds to an ensemble of methods permitting to quantitatively assess the tectonic history of a geological domain.

A variety of restoration processes

Restoration in the large sense covers several complementary processes. The classical methods of unfolding and unfaulting remove the action of tectonic forces. In addition, several techniques can be performed such as the reconstruction of eroded units [e.g., Dimakis et al., 1998, Moucha and Forte, 2011, Boukare et al., 2012, Godefroy et al., 2014], the isostasy compensation [e.g., Allen andAllen, 2013, Lovely et al., 2015], the thermal subsidence due to mantle thermal effect [Royden andKeen, 1980, Allen andAllen, 2013], the rock decompaction due to a change of load [e.g., Athy, 1930[START_REF] Sclater | Continental stretching; an explanation of the post-Mid-Cretaceous subsidence of the central North Sea basin[END_REF], Durand-Riard et al., 2011, Allen and Allen, 2013], and the reverse migration of channelized systems [Rongier et al., 2015, Parquer et al., 2016]. Most of these restoration methods are valuable to assess the evolution of structures and the consistency of a structural interpretation. They have implications in both academic and industrial domains. This thesis focuses on the methods to unfold and unfault a geological model, in particular on the geomechanical approaches [e.g., Maerten and Maerten, 2001, Muron, 2005, Maerten and Maerten, 2006, Moretti et al., 2006, Guzofski et al., 2009, Durand-Riard et al., 2013].

Structural restoration: definition and objectives

Structural restoration is the process which consists in unfaulting and unfolding a geological model. It aims to recover the paleo-geometry of geological units through time [e.g., Chamberlin, 1910, Dahlstrom, 1969, Gratier, 1988, Rouby et al., 2002, Maerten and Maerten, 2006, Moretti et al., 2006]. A myriad of methods exists, based on several assumptions such as the deformation style and the partial knowledge a priori of the paleo-geometries. In addition to recovering the history of a geological domain through time, structural restoration has additional goals. A primary objective of structural restoration is the validation of structural interpretations. Geological data are from outcrops, seismic recordings and exploration wells. They are sparse, limited and uncertain. Moreover, the data analysis is equivocal due to the subjectivity of the interpretation made by a geologist and to the current limit of the knowledge in geosciences [e.g., [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF], Bond et al., 2007, Wellmann et al., 2010, Bond, 2015, Cherpeau and Caumon, 2015]. Restoration of a geological model permits to test the structural interpretations. A restored model presenting inconsistencies discredits the initial interpretation which needs to be reworked [e.g., Dahlstrom, 1969, Moretti and Larrère, 1989, Groshong, 2006]. Another objective of structural restoration is to assess and understand tectonic deformations, for instance to understand the folding process [Plesch et al., 2007, Guzofski et al., 2009, Li et al., 2013], to detect uplift events as in the Annot reservoir analog by Durand-Riard et al. [2011], to analyze faulted salt dome [START_REF] Al-Fahmi | Restorations of faulted domes[END_REF], or to understand the collision stages of the Taiwan orogen [Mouthereau et al., 2001]. Finally, assessing the deformation, and consequently the stress, through time permits to localize fracture areas which are prone to fluid circulation and mineral crystallization [e.g., Macé, 2006, Maerten and Maerten, 2006, Mejía-Herrera et al., 2014, Stockmeyer et al., 2017].

Geometric/kinematic restoration methods: definition and limits

Classical structural restoration methods are based on geometric and kinematic assumptions. They have been applied from the beginning of the last century to nowadays [e.g., Chamberlin, 1910, Dahlstrom, 1969, Gratier, 1988, Moretti et al., 1990, Rouby et al., 2000, Groshong, 2006, Medwedeff et al., 2016, Fossen, 2016]. A multitude of methods working on cross sections, horizon surfaces and volumetric models exists.

Cross-section balancing

The restoration of cross sections, called balancing method, was the first ensemble of restoration methods. The first balanced cross sections were on compressive contexts and based on geometric assumptions of area conservation and potentially conservation of bed lengths and/or bed thicknesses [e.g., Chamberlin, 1910, Dennison and Woodward, 1963, Dahlstrom, 1969, Kiefer and Dennison, 1972, Hossack, 1979, Mitra and Namson, 1989, Groshong, 2006]. Gibbs [1983] extended the balancing of cross sections to extensive contexts, followed by several authors [Davison, 1986, Rowan andKligfield, 1989]. The major drawback of cross-section restoration methods is that they assume that the deformation is within the plane of the cross section, i.e., the out-of-plane deformation is considered as minimal. Typically, strike-slip deformation cannot be handled by such methods.

Map/multi-map restoration

Map restoration has been developed to study horizontal deformations [e.g., Cobbold and Percevault, 1983, Gratier et al., 1991, Gratier and Guillier, 1993, Rouby, 1994, Samson, 1996, Massot, 2002, Dunbar and Cook, 2003, Ramón et al., 2015]. A horizon is represented by a folded and faulted surface. Map restoration aims to flatten this surface, based on geometrical conservations, and to cancel the action of faults. The displacement provided by such restoration permits to analyze the strain and the stress along the horizon surface. Multi-map restoration has been developed to restore a stack of surface horizons [Samson, 1996, Léger et al., 1997, Williams et al., 1997, Rouby et al., 2000, Griffiths et al., 2002]. From the restoration of the uppermost horizon, the restored state of the underlying horizons is determined with the help of a kinematic law such as flexural slip or simple shear. Although the horizon surfaces (2D) are represented in a 3D space, such methods are qualified as 2.5D methods since the restoration mechanism is not a real 3D deformation applied on the entire volume.

3D restoration methods

Several works extended the geometrical restoration in 3D [Massot, 2002, Muron, 2005, Medwedeff et al., 2016], accounting for internal deformation within geological units. All these methods are based on the minimization of the deformation and on volume conservation. Massot [2002] restores the uppermost horizon and then interpolates the displacement field within a regular grid which represents the geological volume. Such an interpolation is performed by Discrete Smooth Interpolation DSI [Mallet, 1989[Mallet, , 1992[Mallet, , 1997]]. Muron [2005] extends this approach in tetrahedral meshes to respect more complex topologies. Another volumetric method is the one developed by Medwedeff et al. [2016]. This method relies on the Geo-Chronological model, abbreviated GeoChron, which only permits two styles of defomation (minimal deformation and flexural slip) [Moyen et al., 2004, Mallet, 2004, Moyen, 2005, Mallet, 2014].

The avenue of geomechanical restoration methods: definition and limits

The necessity of geomechanical approaches

The previously described restoration processes simplify the real rock mechanics by an idealized one. For example, kinematic rules such as simple shear or flexural slip are used. This imposes a style of deformation whereas it is not perfectly known in nature. Thus, one of the drawbacks of classical restoration methods is that internal rock deformation must be known by the geologist and simplified. In addition, the fundamental principle of mass conservation is reduced to surface or volume conservation whereas it is recognized that rock density may change due to geological processes [Ramsay andWood, 1973, Hossack, 1979] even if decompaction is applied in particular to extensive contexts [e.g., Athy, 1930, Allen andAllen, 2013]. Moreover, these restoration methods are largely used in two dimensions (cross-section and map/multi-map restorations) although 3D methods are more and more developed [Massot, 2002, Muron, 2005, Medwedeff et al., 2016]. Several authors highlighted the necessity of incorporating mechanics into geological modeling [e.g., Fletcher and Pollard, 1999, Gjerde et al., 2002, Muron, 2005, Maerten and Maerten, 2006, Guzofski et al., 2009, Al-Fahmi et al., 2016].

The mechanics-based restoration

Mechanics-based restoration corresponds to a boundary value geomechanical simulation. The internal deformation is not known a priori . The resulting strain is a consequence of both rock mechanical behaviors and applied boundary conditions. This approach has been being developed since the 2000's [e.g., Maerten and Maerten, 2001, Santi et al., 2002, Muron, 2005, Moretti et al., 2006, Maerten and Maerten, 2006, Plesch et al., 2007, Guzofski et al., 2009, Durand-Riard et al., 2010, 2013, Tang et al., 2016]. It provides for instance valuable insight on the deformation mechanisms, the paleo-geometry and the potential fracture areas [Macé, 2006, Maerten and Maerten, 2006, Mejía-Herrera et al., 2014, Stockmeyer et al., 2017]. The geological model is filled by elastic rock mechanical properties simulating the response of rocks due to a mechanical constraint. The restoration displacement is not governed by any geometric or kinematic rule, but by fundamental laws of motion: mass conservation and linear momentum conservation. Boundary conditions, generally displacement conditions, are defined to constrain the restored geometry from several assumptions. The common hypotheses are (1) the uppermost horizon was flat and horizontal at its deposition time, and (2) it was not faulted (continuous horizon). Other constraints can be defined to respect a specific geological knowledge such as a direction of deformation.

Physical issues

Although mechanics-based restoration enabled several advances in structural geology, the geomechanical restoration presents many issues. The boundary conditions, mostly set to unfold and unfault a model, are often qualified as unphysical [Lovely et al., 2012]. Moreover the choice of the boundary conditions is not unique and is at the discretion of the geologist. This problem of boundary conditions is recurrent in geomechanical restoration and often addressed [Durand-Riard, 2010, Lovely et al., 2012, Durand-Riard et al., 2013]. Durand-Riard et al. [2013] suggest that classical boundary conditions may fail to properly restore strike-slip and highly oblique-slip faults. They propose to add piercing points, i.e., known fault slip, to improve restoration results. However, that assumes the knowledge of such fault constraints (channel offsets for instance) which may not be always available. In addition, such local constraints must be numerous enough and may provide unrealistic strain within the neighborhood of the piercing points as pointed by Durand-Riard et al. [2013]. Durand-Riard et al. [2013] also show that constraining restoration by a displacement condition on a wall corresponding to the opposite forward displacement provides a more accurate recovery of the expected strike-/oblique-slip for the entire structure. Similarly, Lovely et al. [2012] found that applying a lateral displacement, which is equal in magnitude but of opposite sign to the forward tectonic load, yields a strain field consistent with the forward strain patterns. The classical boundary conditions do not permit such results. Moreover, the geomechanical restoration only uses elasticity and neglects other mechanical behaviors such as plasticity or visco-elasticity [Gerbault et al., 1998]. Due to the physical issues faced by the geomechanical restoration methods, the question of their ability to properly recover the paleo-deformation remains. Even if several studies show the benefits of mechanics over the geometric and kinematic methods [Maerten and Maerten, 2006, Plesch et al., 2007, Guzofski et al., 2009], there is no clear guidelines on the choice between these two restoration families.

Interactivity problem

Most of the mechanics-based restoration methods need a volumetric mesh of a structural model, i.e., a boundary representation of the geological domain with horizon and fault surfaces as boundaries [e.g., Muron, 2005]. The generation of such a model is complicated and is an active research topic [Frank et al., 2007, Pellerin et al., 2014, Zehner et al., 2015, Anquez et al., 2016, Botella, 2016]. Although Durand-Riard [2010] and Durand-Riard et al. [2010] relax the meshing constraints by defining horizons by isovalues of a scalar field, the meshing is still problematic for complex fault networks [Vidal-Royo et al., 2012, Pellerin et al., 2015]. This limits the applicability of the geomechanical restoration to be used as a validation tool due to the time required to edit an initial interpretation and to remesh it.

Diversity of restoration methods

There are several geomechanical restoration methods. All of them rely on fundamental equations of continuum mechanics and on elasticity. The divergences concern the numerical approach. The three main approaches are based on the boundary element method [Gjerde, 2002], on the mass-spring method [e.g., Terzopoulos et al., 1987, Provot, 1995, Gibson and Mirtich, 1997, Bianchi et al., 2003, Shackleton et al., 2008, Midland Valley, 2017a,b], or on the finite element method [e.g., Zienkiewicz and Taylor, 2000a,b, Muron, 2005, Maerten and Maerten, 2006, Moretti et al., 2006, Belytschko et al., 2013, Tang et al., 2016]. Concerning the finite element method, it is subdivided into several approaches: the dynamic relaxation [e.g., Underwood, 1983, Oakley and Knight, 1995, Muron, 2005], the global static approach [Muron, 2005, Moretti et al., 2006], and the local static approach [Maerten and Maerten, 2006]. All the different approaches are described in Chapter 1. The choice of a specific numerical approach has numerous implications. It influences for instance: the type of mesh to represent a geological model and the subsequent meshing constraints, the geomechanical simplifications, the implementation ease, the computer memory requirement, and the computational time of a restoration process. Concerning the fault contact boundary conditions, there are two primary methods. One is purely geometric [Muron, 2005, Moretti et al., 2006], and the other relies on the mechanical contact theory [Muron, 2005, Wriggers and Laursen, 2006, Maerten and Maerten, 2006, Tang et al., 2016]. Considering the multiplicity of the geomechanical restoration methods, it may be difficult for a geologist to choose the most appropriate method for a specific case study [Moretti and Titeux, 2007].

Problematic, approach and contributions of this thesis Problematic and manuscript organization

We show that the mechanics-based restoration has several limits. In this thesis, we provide a variety of guidelines for the effective and robust use of geomechanical restoration. In particular we focus on the following questions:

• What are the limits and the advantages of each geomechanical restoration method?

• Which boundary conditions provide a restored state that is physically and geologically acceptable?

• What is the place of the geomechanical approaches in comparison to the traditional geometrical approaches?

• What is the impact of elastic mechanical parameters on the restoration?

This manuscript is composed of three chapters. Chapter 1 is a review of the different geomechanical methods to unfold and unfault a geological model. The main objectives are: (1) to describe how these methods work and (2) to highlight the differences between them and what they imply for a structural geologist. Chapter 2 deals with the restoration of a structural sandbox model analogous to supra-salt extensional structures. The aim of this study is to define proper boundary conditions for the restoration of the analog model knowing its paleo-geometries on a cross section. Chapter 3 is a comparison of two restoration methods: mechanics-based and GeoChron-based. This comparison is performed on the structural sandbox model. The knowledge of the paleo-geometries leads to an objective comparison between both methods. These chapters aim to (1) ease the choice of a restoration method, and (2) define proper boundary conditions in mechanics-based restoration.

The presented work on the mechanics-based restoration is the continuity of the work of Muron [2005] and Durand-Riard [2010]. The geomechanical tool used in the PhD has been developed in a C++ mechanical library, based on the finite element method, by Antoine Mazuyer and Benjamin Chauvin. This external library called RINGMecha4 [Chauvin and Mazuyer, 2016] is the successor of RestorationLab5 which is a SKUA-GOCAD plugin [Paradigm, 2015]. A significant part of this thesis consisted in rewriting the code of RestorationLab into RING-Mecha to:

• Update and clean the code after several years of development.

• Check all parts of the code and fix existing issues.

• Convert the code into a standalone program without graphical interface. In this way, all the source code is controlled and easy to script.

• Adapt the code to current programming (e.g., parallel solvers and mesh libraries

Introduction

Several methods have been developed for the mechanics-based restoration. The variants may differ in the used numerical procedure, in the mesh type, and in the way to handle boundary conditions. These differences may impact the restored state, geometrically or physically, and computational time or computer memory requirements. The choice of a specific method is not trivial. This chapter is a review of the different geomechanical restoration approaches to unfold and unfault a geological structural model. This review is theoretical, i.e., there is no comparison between two methods on a natural case study. We base this review on the pusblished mechanics-based restoration methods and on the used mathematical founding principles. We hope that this review will ease the understanding of the geomechanical restoration in general and of the subsequent methods. For simplicity, we have placed most mathematical details in Appendix 1.A.

1.1 Mechanics-based restoration: toward a restoration process including rock mechanical behavior

Traditional restoration approaches are based on geometric and kinematic assumptions [e.g., Chamberlin, 1910, Dennison and Woodward, 1963, Dahlstrom, 1969, Kiefer and Dennison, 1972, Hossack, 1979, Cobbold and Percevault, 1983, Gratier, 1988, Rowan and Kligfield, 1989, Gratier et al., 1991, Rouby, 1994, Samson, 1996, Léger et al., 1997, Williams et al., 1997, Rouby et al., 2000, Griffiths et al., 2002, Rouby et al., 2002, Massot, 2002, Muron, 2005, Groshong, 2006, Mallet, 2014, Ramón et al., 2015, Medwedeff et al., 2016, Fossen, 2016]. Although these methods have provided valuable insight on paleogeometries, and the subsequent interpretations, on many field cases, several authors argue that mechanics must be taken into account to properly handle geological processes, in particular in restoration [e.g., Fletcher and Pollard, 1999, Maerten and Maerten, 2001, Gjerde, 2002, Gjerde et al., 2002, Muron, 2005, Moretti et al., 2006, Maerten and Maerten, 2006, Guzofski et al., 2009]. This section presents the physical and geometrical bases of such a mechanics-based restoration.

A restoration based on continuum mechanics

Instead of defining the motion by purely kinematic or geometric laws, as in classical restoration methods, geomechanical restoration uses continuum mechanics and material constitutive laws to evaluate structural deformation. In mechanics-based restoration, the geological domain is represented as a solid which is mechanically deformed backward to its paleo (restored) geometry. Two fundamental laws, which stem from continuum mechanics, govern the motion. Let Ω and m respectively denote the geological domain to restore and its mass. The commonly accepted assumption in restoration and more globally in rock mechanics is that the domain keeps a constant mass through time. Thus in continuum mechanics, and so in geomechanical restoration, the first law is the mass conservation [e.g., Malvern, 1969[START_REF] Marsden | Mathematical foundations of elasticity[END_REF], Wriggers and Laursen, 2006, Belytschko et al., 2013]:

dm dt = 0, (1.1)
where dm dt denotes the derivative of the mass m by the time t. In continuum mechanics, solid deformation and motion are characterized by the conservation of linear momentum, which corresponds to Newton's second law. The linear momentum p of a domain Ω is defined by its mass m multiplied by its velocity u. Let f denote the sum of the forces applied on the domain. The conservation of linear momentum is described by:

dp dt = f . (1.2)
Both laws are combined into a single equation which forms the basis of all the geomechanical restoration methods [Belytschko et al., 2013, p. 115-118]:

∇ • σ + ρb = ρü, (1.3)
with ∇• the mathematical divergence, σ the Cauchy stress tensor, and b the body forces. The body forces concretely correspond to the gravity force in solid mechanics. In the dynamic relaxation method, which is a numerical method to solve a restoration problem by Muron [2005] (Section 1.3.2), a damping term is added to the body forces. Equation (1.3) defines the mechanical equilibrium of the restored state.

A restoration based on elasticity assumption

Restoration aims to obtain a reasonable and consistent geometry of a geological model after removal of tectonic deformations. By definition it is the reverse of a forward process. As forward deformation is irreversible (plasticity) and the paleo-stresses is poorly known, it is not possible to directly reverse time. Elasticity, which is reversible, is used in mechanics-based restoration by simplicity and as a simple assumption [Muron, 2005, Maerten and Maerten, 2006, Moretti et al., 2006, Guzofski et al., 2009]. A justification is that elastic deformation is a reasonable approximation for incremental deformations recorded in growth strata, and restored sequentially. According to Maerten and Maerten [2006], even if rock reaction is not completely elastic and that a part of the information is potentially lost, it can be compensated by some mechanical assumptions. In addition, during a sequential restoration the stress is not accumulated, i.e., the stress is reset to zero before each restoration step. This relaxation of the stress avoids reaching high and unphysical strain [e.g., Maerten andMaerten, 2006, Guzofski et al., 2009]. Indeed, in nature the stress decreases by plasticity, fracturing, and ductile deformation. Thus the stress relaxation in restoration aims to mimic the nature stress dissipation through time even if it cannot be perfectly equivalent.

To be able to solve Equation (1.3), the Cauchy stress σ must be defined, here by an elastic law. Small deformation assumption is often used for simplicity. It simplifies the equations by a linear relationship between the Cauchy stress σ and the linear Green-Lagrange strain (which reduces the computational time). This relation is given by Hooke's law [e.g., Ramsay and Huber, 2000]:

σ = D : , (1.4)
with D the elastic fourth-order tensor and : the double dot product (Appendix 1.B.1). The small deformation assumption is valid as long as the quadratic terms in the strain equation are negligible. Salençon [2005] defines that the norm of the gradient of the displacement u, here the restoration displacement, relatively to the unrestored configuration ∇ 0 must be far lower than 1:

∇ 0 u 1. (1.5)
Concretely, if we consider a 1D solid, the deformation must be far smaller than the size of the original solid. In practice, 10% of deformation as a maximum should be considered. It is important to note that the small deformation assumption can be used even if the displacement is large. In that case, most of the displacement corresponds to a translation and/or rotation which has no effect on the strain, i.e., a solid under a pure rigid motion has a null strain [Fossen, 2016, p. 24].

As tectonic deformations are not small in general, the small deformation assumption may be invalid in the mechanics-based restoration [Moretti and Titeux, 2007]. Muron [2005], Moretti et al. [2006] and Moretti and Titeux [2007] use hyperelastic laws to relate the stress to the strain, which means that the work does not depend on the path of deformation [Belytschko et al., 2013]. Such laws are not linear, making the solving of the restoration more complex and more time consuming. Muron [2005] uses the Neo-Hookean law (Appendix 1.B.3) which is an extension of Hooke's law for large deformations [e.g., Belytschko et al., 2013, Zienkiewicz and Taylor, 2000b, Muron, 2005]. Moretti and Titeux [2007] use the Simo-Miehe law [Lorentz, 2013].

Constitutive laws incorporate mechanical rock properties which are defined within the geological model to represent rock behavior. Mechanics-based restoration generally uses isotropic elastic parameters which can be described with Young's modulus E and Poisson's ratio ν [e.g., Sokolnikoff, 1956, Gjerde, 2002, Muron, 2005, Moretti et al., 2006, Maerten and Maerten, 2006, Gercek, 2007]. Other equivalent parameters can be used such as the first Lamé parameter λ, the second Lamé parameter µ (also called shear modulus) or the bulk modulus K (Appendix 1.C). These parameters may vary within the domain and are isotropic (no direction dependency). It is usual to define different elastic parameters per geological layer, to reflect the nature of rocks. It is important to note that the change of volume from the unrestored state to the restored state is controlled by Poisson's ratio and the applied boundary conditions. Indeed, the mass is conserved and the mechanical deformation is minimized but the volume is not necessary maintained.

It is well-known that sedimentary layers may slide along bedding limits (flexural slip) corresponding to thin weak beds in particular in compressive contexts [e.g., Donath and Parker, 1964, Dubey and Cobbold, 1977, Tanner, 1989, Fossen, 2016]. To emulate this effect in restoration, Durand-Riard [2010] and Durand-Riard et al. [2013a] upscale a stack of geological layers, under flexural slip mode, into a single layer defined by a transverse isotropic material [Backus, 1962, Salamon, 1968, Watkinson and Cobbold, 1981, Crea et al., 1981, Graham and Houlsby, 1983, Chalon et al., 2004, Titeux, 2009, Proix, 2010]. This avoids the explicit introduction of several sliding interfaces or the definition of multiple thin weak layers, which would be a meshing constraint and imply irrealistic mesh refinement [Durand-Riard et al., 2013a]. Consequently, this avoids a significant increase of the computer memory cost and of the computational restoration cost. A transverse isotropic mechanical behavior is modeled by five elastic parameters (Appendix 1.D.1). The upscaling proposed by Durand-Riard et al. [2013a] for restoration purposes stems from the approach of Salamon [1968]. The five elastic parameters derive from the elastic parameters of each layer in the stack under flexural slip deformation, including the thin weak beds and their proportion (Appendix 1.D.2). A difficulty of this approach is to evaluate the volumetric proportion of the weak beds relatively to the entire stack. In addition, the use of pure transverse isotropy implies that the stratigraphy is perfectly horizontal, i.e., a bias is introduced in the case of non-horizontal layers. Tilted transverse isotropy enables to consider the stratigraphy normal, and not the depth axis, as the anisotropy axis [Gornet, 2008, Danek et al., 2010, Kostecki, 2010, 2011].

Boundary conditions to ensure a geologically consistent restored state

The equation of motion (1.3) defines the evolution of the solid as a result of mechanical constraints. These constraints are the boundary conditions of the geomechanical restoration. They ensure that the restored state is geologically consistent with external knowledge, and that the restored state is unique. The classical boundary conditions in restoration aim to reverse the action of tectonic forces: unfolding and unfaulting (Figure 1.1). They are mainly displacement (Dirichlet) conditions. Indeed, traction (Neumann) boundary conditions which would enable to unfold and unfault are unknown and difficult to assess [Lepage et al., 2004, Muron, 2005, Moretti et al., 2006]. A null Neumann boundary condition is defined on the boundaries in which no Dirichlet or contact condition is defined (free boundaries).

Displacement boundary conditions: unfolding and numerical stability

Displacement boundary conditions are largely used in mechanics-based restoration. Their main use is to flatten the uppermost horizon, i.e., the depth component of the top horizon is imposed at a specific value. Another use of Dirichlet conditions is to restrict the motion along a specific direction. For instance, in the case of evidence of a direction without significant deformation, the geological model is imposed to not move along that direction. Durand-Riard [2010] prevents two lateral walls from moving along the horizontal direction orthogonal to the main horizontal direction of folding in a fault-bend fold model. Similarly, Chauvin et al. [2017] apply such conditions on an extensional model. In some cases, a part of the model is constrained to move with a specific displacement. Durand-Riard [2010], Durand-Riard et al.

[2013b] and Chauvin et al. [2017] apply a displacement condition on a wall respectively in a compressive context, a strike-slip context and an extensional context. Finally, some numerical methods need to have enough constraints to converge into a unique solution, by avoiding a multitude of equivalent restoration solutions which would only differ by a global translation or a global rotation [e.g., Lepage et al., 2004, Moretti et al., 2006, Durand-Riard et al., 2013a,b].

Unfaulting based on geometrical and/or physical boundary conditions

Fault contact boundary conditions are defined to remove the action of faults. There are three main kinds of fault contacts (represented in Figure 1.1). The first one consists in avoiding gap and overlap between the fault blocks. This condition ties the two mirrors of a fault. The second condition ties both cutoff lines of a horizon. This condition is generally only applied on the uppermost horizon to handle synsedimentary deformation. Indeed, in the case of growth stratigraphy, all the horizons did not undergo the same amount of slip. For the horizons below the uppermost horizon, the cutoff line geometry in the restored state is a consequence of the mechanical simulation. Note that in the case of pre-growth stratigraphy, all the horizon cutoff lines for a same fault are a priori removed at the same time. Finally, when two clear markers are visible on each side of a fault, for instance two markers of a same channel path, a specific contact condition can be defined to tie both points. Durand-Riard et al. [2013b] call such contact points piercing points. These points represent a meshing constraint since two nodes must exist on these locations. Contact conditions are performed by Dirichlet conditions [Muron, 2005, Moretti et al., 2006] or by contact mechanics [Muron, 2005, Maerten and Maerten, 2006, Wriggers and Laursen, 2006]. The latter is a mix between Dirichlet and Neumann conditions.

Numerical representation of a geological model in restoration

Equation ( 1.3) has no analytical solution in the general case. Thus to solve it and get the restoration displacement, it is necessary to use a numerical method. Several numerical methods have been used in restoration. They need a numerical representation of the geological model. We develop here only the geological representation used in the published restoration methods. There are two main representations of a geological model. The first one is the boundary representation (Figure 1.3a) [Caumon et al., 2004[Caumon et al., , 2009]]. It is a representation by surfaces. Each surface represents a geological surface, generally a fault or a horizon. The ensemble of these surfaces constitutes a sealed volume called a structural model. The second representation is a 3D mesh composed by 3D volumetric elements (Figure 1.3b). In this representation the volumetric elements are generally conformal to horizons and faults. The volumetric elements are most often tetrahedra since these elements are able to fit complex geometries [Muron, 2005, Zehner et al., 2015]. When the volumetric elements are not conformal to the horizons, the latter can be represented by isovalues of a scalar field within the mesh [e.g., Frank et al., 2007, Durand-Riard, 2010, Durand-Riard et al., 2010]. Thus, the volumetric elements are crossed by the horizons, leading to an approximation of the elastic parameters along the horizons since these parameters are stored within each element (Figure 1.2). For a tetrahedron, Durand-Riard et al. [2010] define average elastic parameters by a weighted mean. Let E r and ν r denote Young's modulus and Poisson's ratio of the red layer in Figure 1.2, and E b and ν b denote Young's modulus and Poisson's ratio of the blue layer. The cut of a tetrahedron by the implicit horizon defines a volume V r of the tetrahedron within the red layer, and a volume V b of the tetrahedron within the blue layer. Average Young's modulus E av and Poisson's ratio ν av in the tetrahedron are [Durand-Riard et al., 2010]:

E av = V r × E r + V b × E b V r + V b ; ν av = V r × ν r + V b × ν b V r + V b . (1.6)
Finally, above the implicit uppermost horizon to restore, the mesh has no geological meaning and is only present to store the scalar field which defines the stratigraphy. Durand-Riard et al. [2010] define a rubber-like material to avoid a mechanical artifact of this part of the mesh (E = 0.2 GPa and ν = 0.5).

The typical workflow (Figure 1.3) in geomechanical restoration which includes the geological representation construction is: (1) the construction of a structural model (for instance from seismic data) (Figure 1.3a), ( 2) the meshing of the structural model by volumetric elements if necessary (Figure 1.3b), (3) the assignment of mechanical properties within the geological model (Young modulus and Poisson's ratio or equivalent, Figure 1.3c), ( 4) the definition of the boundary conditions to unfold and unfault (Figure 1.3d), ( 5) the resolution of Equation ( 1.3) by a numerical method (Figure 1.3e), i.e., the restoration displacement field is computed, (6) the application of the restoration displacement field to get the restored model (Figure 1.3f), (7) the removal of the uppermost layer (Figure 1.3g). In the case of poly-stage deformations or of growth strata, if it remains at least a geological layer, the workflow runs until no layer remains from the step ( 4). This workflow may vary for certain restoration methods.

Numerical methods for solving a geomechanical problem

In this section we present the different numerical methods which are used in geomechanical restorations to solve Equation ( 1.3). The main families of numerical methods are the boundary element method [Beer and Poulsen, 1994, Gjerde, 2002, Maerten et al., 2006, Marshall et al., 2008, Maerten et al., 2010], the mass-spring method [Terzopoulos et al., 1987, Provot, 1995, Shukla and Jayakumar, 2011, Liu et al., 2013, Macaulay et al., 2015], and the finite element method [e.g., Zienkiewicz and Taylor, 2000a,b, Maerten and Maerten, 2001, Santi et al., 2003, Muron, 2005, Moretti et al., 2006, Maerten and Maerten, 2006, Hughes, 2012, Belytschko et al., 2013, Bathe, 2014, Tang et al., 2016]. The latter is the most common method in mechanics-based restorations and is itself subdivided into several methods.

Geomechanical restoration solved by the finite element method

In this section we describe the principle of the finite element method. Santi et al. [2003], Muron [2005], Maerten and Maerten [2006], Moretti et al. [2006], Durand-Riard [2010], Tang et al. [2016] and Chauvin et al. [2017] developed geomechanical restoration tools using this numerical method.

Principle

The finite element method is a mathematical method to approximate the solution of differential equations. This numerical approach discretizes the studied medium into elements. A 3D mesh composed by volumetric elements is constructed (Figure 1.3b). The differential equations are then solved on each element. The global solution is obtained by the concatenation of all the elementary solutions.

The finite element method corresponds to a family of methods. The classical steps of the finite element method are:

• The conversion of the strong form into a weak form: Equation (1.3) is transformed into an equation which can be solved.

• The discretization of the weak form: the medium is split into elementary media. • The passage to reference elements: equations defined on the elementary media are transformed to be defined into reference elements in which interpolation functions are well known.

a) b) c) e) f) g) d) E1, ν1 E2, ν2
• The integrals in the equations are approximated by sums (Gauss quadrature).

• All the matrices, defined on each element, are concatenated into global matrices representing the entire geological domain. This step is called assembly.

• The resolution of the problem: the displacement is computed using the static method (Section 1.3.1) or the dynamic relaxation method (Section 1.3.2).

Application of the finite element method to the mechanics-based restoration

Equation ( 1.3) is defined on the restored domain Ω. This space is not known since it is the purpose of the restoration. In the following development, we choose to convert Equation ( 1.3) into the unrestored space as done by Muron [2005] since such a space geometrically remains the same through the restoration simulation. This formulation is called the total Langrangian formulation [e.g., Muron, 2005, Belytschko et al., 2013, Bathe, 2014]. Other formulations, which should provide the same solution, are possible (e.g., the updated Lagragian formulation) but we just develop this one. The formulation of Equation ( 1.3) on the unrestored geometry is (Appendix 1.E):

∇ 0 • P + ρ 0 b = ρ 0 ü, (1.7)
with P the nominal stress (Appendix 1.B.2), ρ 0 the unrestored rock density, ∇ 0 • the divergence in the unrestored space, b the body forces, and ü the acceleration vector. Note that this new equation is not the equilibrium of the unrestored state but the equilibrium of the restored state expressed in the unrestored space.

Equation (1.7) is called the strong form. The first step of the finite element method consists in transforming this strong form into a weak form which can be solved. This weak form is mathematically and physically equivalent to the strong form. This transformation is performed by the virtual work principle [e.g., Dym and Shames, 1973, Zienkiewicz and Taylor, 2000a,b, Muron, 2005, Wriggers and Laursen, 2006, Belytschko et al., 2013, Bathe, 2014]. The virtual work principle introduces a virtual displacement δu which is applied on the solid at equilibrium state. A light deformation of a solid enables to get its characteristics such as its stress. This permits to get the real displacement from the unrestored state to the restored state. In the case of elasticity in solid mechanics, the principle of virtual work implies the principle of least action [Dym and Shames, 1973, p. 124-127]. Thus, it corresponds to a minimization of the deformation, from the unrestored state to the restored state. Equation (1.7) is homogeneous to a force. A simple way to understand how to get the equivalent weak form is to multiply Equation (1.7) by the virtual displacement (a work correspond to the product between a force and a displacement) and then to integrate over the entire body. After some basic mathematical developments the final shape of the weak form is (Appendix 1.F):

Ω 0 δu • ρ 0 ü dΩ 0 + Ω 0 (∇ 0 (δu)) T : P dΩ 0 = Γ 0 δu • t 0 dΓ 0 + Ω 0 δu • ρ 0 b dΩ 0 .
(1.8) At this stage, there is no approximation yet. The weak form needs to be discretized because of three main reasons. Firstly, it is defined continuously on a domain Ω 0 . So to get the solution at several locations a discrete domain must be considered. Secondly, the integrations are on a domain whose geometry is arbitrary. A simplification of this geometry Ω0 is needed to compute integrals in (1.8). Thirdly, the different mechanical parameters such as the density, Young's modulus and Poisson's ratio, are not necessary constant in the whole domain. Classically, these parameters are different for each geological layer and potentially within them. Therefore the medium Ω 0 is split into several media Ωe 0 such as (Figure 1.4):

Ω 0 ≈ Ω0 = ne e Ωe 0 ; ∀ (i, j) ∈ 1; n e 2 /i = j | Ωe 0 | e=i ∩ Ωe 0 | e=j = ∅. (1.9)
Such a discretization is called a mesh. Here we can see that due to the dicretization, the geological model boundaries (faults and horizons) are approximated and that the conformity to the structure is dependent on the mesh size. Thus, the finite element method, by the discretization, approximates the solution of a differential equation [e.g., Lo, 2002]. Moreover, the shape of the elements may have an important impact on the result. For a good quality, equilaterality in solid mechanics is generally wished [e.g., Parthasarathy et al., 1994, Shewchuk, 2002, Munson, 2007]. The weak form is discretized using: (1.10)

Similar relations for surface integrals may be established. Thus, instead of solving the weak form on the whole domain, it is locally solved on each element. For each element a local system of equations is defined, such a system is represented by elementary matrices. Several finite element solvers assemble the matrices of all the elements into global matrices to get a single system of equations to solve: .11) with M the mass matrix, F int the vector of the internal forces, F ext the vector of the external forces, C the damping matrix, and u the velocity. Note that the two latter terms stem from the body forces b and are used in the dynamic relaxation method (Section 1.3.2). Geomechanical restoration approaches relying on a finite element procedure solve the above equation, in this form or another. As the number of methods is high, Section 1.3 is dedicated for their description.

M • ü + C • u + F int -F ext = 0, ( 1

Geomechanical restoration solved by the boundary element method

The boundary element method has seldom been used to perform a geomechanical restoration problem [Gjerde, 2002], even if this approach is very common in geomechanical applications [e.g., Beer and Poulsen, 1994, Maerten et al., 2006, Marshall et al., 2008, Maerten et al., 2010]. This numerical method just requires a boundary representation, i.e., a mesh composed by volumetric elements (e.g., tetrahedra) is not required. Thus, for geological applications, just the structural model is needed (Figure 1.3a). Equation (1.3) is thus solved only on the surfaces and then the solution can be computed within the volume on specific points. In the boundary element method the surfaces are discretized, introducing a geometrical approximation. As in the finite element method, the equation of motion is solved in each (surface) element. The global solution is obtained by assembling all the local ones.

The use of just a boundary representation instead of a full volumetric mesh is a huge advantage of the boundary element method over the finite element method. Indeed, it reduces the meshing constraints. Moreover, the finite element method solves the equation of motion everywhere (dependent of the mesh discretization) even in areas in which the solution is not needed. Conversely, the solution within the volume is performed only on necessary points in the boundary element method. Furthermore, as only the surfaces are discretized, there are far less mesh nodes. Thus, the boundary element method requires less computer memory and less computational time than the finite element methods. However, the resulting system of equations implies matrices which are not sparse and not symmetrical [e.g., [START_REF] Gangming | A new boundary element method coupled with FEM packages[END_REF]. The fact that the matrices are not sparse is not really a problem, since the number of degrees of freedom are farly reduced in comparison to the finite element method. Finally, the boundary element method is more appropriate for linear problems with homogeneous mechanical properties within a close volume. Mechanical heterogeneities and non-linearities are more difficult to handle than in the finite element method. As a consequence, a volumetric discretization by volumetric elements may be necessary.

Geomechanical restoration solved by the mass-spring method

The mass-spring method [e.g., Terzopoulos et al., 1987] is used by Midland Valley [2017b] for restoration purposes [Shackleton et al., 2008, 2011, Ghail, 2014, Macaulay et al., 2015, Verdon et al., 2015]. There are few information on this method. In this section, we try to provide the bases of this mass-spring approach from the literature. We hope that the general understanding of this numerical approach will ease the comprehension of the specific geomechanical restoration tool of Midland Valley [2017b].

The mass-spring method solves Equation (1.3) using a set of mass points connected by springs (Figure 1.5), in a boundary representation or in a volumetric model composed by volumetric elements [Shackleton et al., 2008, Midland Valley, 2017a]. Thus this approach approximates a physical 3D problem into an ensemble of 1D problems, making this approach faster than the finite element method but less realistic [Gibson and Mirtich, 1997, Bourguignon and Cani, 2000, Bianchi et al., 2003]. The boundary conditions initiate the deformation on some points (Figure 1.5c), for instance to flatten the uppermost horizon or to remove fault slip [Midland Valley, 2017a]. This generates forces on the springs sharing the impacted nodes. The internal force

F int 0 , contained in ∇ • σ in (1.
3), acting on a node N 0 corresponds to the sum of the action of the n ∈ 1; n springs connected to this node [e.g., Provot, 1995, Louchet et al., 1995, Bianchi et al., 2003]:

F int 0 = - n i=1 k 0i N 0 N i -N 0 0 N 0 i N 0 N i N 0 N i , (1.12)
with N i a node connected to N 0 by a spring N 0 N i (vector) of stiffness k 0i , and N 0 0 N 0 i the natural length of the spring (when unloaded) N 0 N i (Figure 1.6). These forces are propagated through the mesh, spring by spring (Figure 1.5d). This propagation is controlled by the motion equation and the stiffness of each spring. The resolution of Equation ( 1.3) is pseudo-temporal, i.e., the deformation is propagated through time within the mesh and a damping enables to reduce more and more the energy of the system until equilibrium. In a way, this approach is similar to the dynamic relaxation which is a specific finite element method (Section 1.3.2).

By approximating a 3D problem by a set of springs, the mass-spring method is very sensitive to the choice of the spring system. The points and their masses must be appropriately chosen. Points generally correspond to mesh nodes. The mass of a point can be determinated with the volume of the Voronoi cell of the point in the Voronoi diagram [e.g., Aurenhammer, 1991] of the entire mesh [Bourguignon and Cani, 2000], using a user-defined density. For the definition of the springs, a basic method is to choose the mesh element edges. For integrating mechanical anisotropy, another choice of springs can be defined, leading to springs which do not necessary correspond to mesh element edges [Bourguignon andCani, 2000, Midland Valley, 2017a]. In all cases, the mechanical response of a mass-spring system is very sensitive to the choice of the spring system, and so of the mesh, and of the spring stiffnesses [Bourguignon and Cani, 2000, Bianchi et al., 2003, Midland Valley, 2017a]. This issue is for instance addressed by Bourguignon and Cani [2000] and Bianchi et al. [2003], among other authors. Moreover, the mass-spring method has issues to respect some rules such as the volume conservation [Gibson andMirtich, 1997, Bourguignon andCani, 2000] which is a common rule in restoration in general (even if the geomechanical restoration preserves mass and not volume). Several authors developed approaches to limit the changes of volume after a mass-spring simulation [Lee et al., 1995, Promayon et al., 1996, Bourguignon and Cani, 2000]. For the numerical resolution, several methods exist such as the explicit Euler integration [Provot, 1995] and the implicit integration [Baraff andWitkin, 1998, Liu et al., 2013]. Here the mass points are on mesh nodes, and the springs correspond to the mesh edges. c) Displacement of a node of the mesh (green node) due to boundary conditions. The springs connected to this node (green springs) are in elongation. Blue nodes and blue springs do not move yet. d) The elongation of the springs in c) generates forces, controlled by the spring stiffnesses, acting on the nodes connected to them (green nodes). The motion of these nodes implies a length change of the springs connected to them (green springs). Thus the motion propagates, spring by spring, with a damping which progressively dissipates the energy, until the system reaches an equilibrium.

Figure 1.6: Nodal internal forces in a mass-spring system. The internal forces of a node N 0 correspond to the sum of the action of each spring connected to this node. N i , with i ∈ 1; n , is one of the n nodes connected to the node N 0 by a spring of stiffness k 0i and of natural length N 0 0 N 0 i .

Chapter

Comparison between the different numerical methods

The different fundamental similarities and differences between the numerical methods used in geomechanical restoration are summarized in Table 1.1.

Geomechanical restoration solved by different finite element methods

Many restoration codes are based on a finite element method such as Dynel software developed by Schlumberger [Maerten and Maerten, 2006, Maerten, 2010, Vidal-Royo et al., 2012, Maerten and Maerten, 2015, Schlumberger, 2017], Kine3D sofware which is a GOCAD-SKUA plugin developed by IFPEN [Moretti et al., 2006, Moretti, 2008, Paradigm, 2015], and RING-Mecha1 library developed by the RING team in the laboratory GeoRessources [Muron, 2005, Durand-Riard, 2010, Chauvin and Mazuyer, 2016]. In this section we explain the different methods used in the literature to solve Equation (1.11) for restoration purposes. As far as we know, all the solvers described in the literature to solve a restoration problem do not consider the time. Thus, the transient evolution from the unrestored to the restored states does not matter, only the steady-state is important. There are two methods to get the steady-state while neglecting the transient part. The first method is the static method. That simplifies Equation (1.11) since the temporal terms disappear. The second method is the dynamic relaxation method. This technique keeps the temporal terms to transform Equation (1.11) into a recurrence relation, reducing the computer memory requirements.

Static finite element method

In the static form, mathematical terms in Equation (1.11) are not functions of time. Hence u and ü are null and Equation (1.11) becomes:

F int = F ext .
(1.13) Muron [2005] and Moretti et al. [2006], Moretti and Titeux [2007] deal with Equation 1.13 in a global way, i.e., a unique system of equations for the whole mesh is solved. Thus, the internal forces F int and the external forces F ext are computed for the entire domain. Their definitions depend on the mechanical problem. Fault contact and large deformations are sources of non-linearities. If there is no non-linearity, the internal forces F int are defined by a linear relation according to the displacement (thanks to the Hooke's law and the linear Green-Lagrange strain, Appendix 1.B.1):

1.3.1.1 Global approach
F int = K • u, (1.14)
with K the stiffness matrix representing the mechanical behavior of the entire model. From Equation (1.13), the product between the stiffness matrix K and the displacement vector u is equal to the external force vector F ext :

K • u = F ext .
(1.15)

u can be assessed as it is the unique unknown. The intuitive way to do that is to invert the stiffness matrix K: u = K -1 • F ext . In practice this method cannot be realized in a reasonable time since the stiffness matrix can have a dimension greater than one million in restoration applications. As the stiffness matrix K is sparse and symmetric, the system 1.15 can be efficiently solved by numerical methods such as the Cholesky decomposition or the conjugate gradient [e.g., M.R and Stiefel, 1952, Davis and Duff, 1993, Demmel et al., 1995, Golub and Van Loan, 1996, Schenk and Klaus, 2004, Muron, 2005, Gould et al., 2007, Dureisseix, 2008, Bathe, 2014].

In the case of non-linearities, i.e., if there are fault contacts or if the small deformation assumption is not appropriate, the internal forces are not linear with the displacement u. Thus it is not possible to solve directly Equation (1.13). A classical method to handle that in solid mechanics, and used by Muron [2005], Moretti et al. [2006], and Moretti and Titeux [2007], is to linearize Equation (1.13) [e.g., Zienkiewicz and Taylor, 2000b, Wriggers and Laursen, 2006, Belytschko et al., 2013, Bathe, 2014]. In the literature, the Newton-Raphson algorithm is often used (Figure 1.7). Let f (u) denote the substraction between the internal and the external forces:

f (u) = F int -F ext .
(1.16)

The Newton-Raphson algorithm solves f (u) = 0 iteratively. At iteration k + 1, a Taylor expansion of f is computed from the previous step k:

f (u k+1 ) f (u k ) + ∂f (u) ∂u u=u k • (u k+1 -u k ) = 0.
(1.17)

Then the new solution u k+1 is found. Then on this new solution a new Taylor expansion of f is done and a new solution of u is determined and so on until the difference of the solutions u from one step to another becomes negligible. The derivative of (1.16) is:

A = ∂F int ∂u - ∂F ext ∂u = K int -K ext . (1.18)
A is the Jacobian matrix2 [Belytschko et al., 2013, p. 341-343]. A is a function of the derivatives according to the displacement u of the internal forces F int , called tangent stiffness matrix K int , and of the external forces F ext , called load stiffness matrix K ext . The expression of these derivatives is obtained by linearization [Belytschko et al., 2013, Section 6.4]. The algorithm to solve non-linear system is in Algorithm 1. An initial solution u 0 is necessary to begin the Newton-Raphson procedure. To ensure the convergence, it is necessary to define the initial solution as the best assumption of the final solution. Convergence will depend on this initial solution. The linear solution u linear = K -1 • F ext may be chosen [Muron, 2005]. In addition, other iterative methods than Newton-Raphson may be used such as the modified Newton-Raphson method or the quasi-Newton methods which use second order Taylor expansions [Zienkiewicz and Taylor, 2000b].

Algorithm 1 Static non-linear algorithm. Global approach as developed by Muron [2005] and Moretti et al. [2006] Displacement Maerten and Maerten [2006] solve Equation (1.13) locally instead of globally. They only use the small deformation assumption so they also have to solve Equation (1.15). However, instead of solving it for the entire system at once, they use a Gauss-Seidel algorithm which is an iterative method [e.g., Golub andVan Loan, 1996, Dureisseix, 2008]. For each iteration all the nodes are mechanically analyzed to determine their displacements. The order of analysis is not important [Maerten and Maerten, 2006]. For a node, Equation (1.15) is locally built. Only the surrounding elements of the current node are considered. The external forces applied to the current node correspond to the deformation stress of the surrounding elements. The stiffness matrix is defined on the current node and corresponds to the addition of surrounding element stiffness matrices. The other nodes which belong to the surrounding elements are defined as fixed [Maerten and Maerten, 2006]. Boundary conditions applied to the current node, if any, are taken into account. Finally the displacement of the current node is computed and the node is displaced. Then the algorithm considers another node and the same process is repeated. As pointed by Maerten and Maerten [2006], an important fact is that the current node uses the new position of the moved nodes at the current iteration and not just the nodal positions at the previous iteration. The iterative solver ends when the nodal displacements are negligible.

u ← u linear ; Newton-Raphson step k ← 0 repeat k ← k + 1 Compute r ← F int k-1 -F ext k-1 Compute A ← K int -K ext Integrate Dirichlet boundary conditions into A Solve A • ∆u = -r u ← u + ∆u until ∆u < or k = maximum number of iterations return u 1.3.1.2 Local approach
The method developed by Maerten and Maerten [2006] has several advantages. First, the iterative process enables a better handling of the non-linear constraints such as fault contacts. Second, the Gauss-Seidel algorithm provides an under-constrained system. Thus it is not necessary to define a set of Dirichlet boundary conditions to guarantee the uniqueness of the solution (pin nodes, lines or walls to avoid translation and rotation of the entire model). Finally, as the algorithm is iterative and just a local stiffness matrix is built at each step, the allocated memory is far smaller than the allocated memory of a global stiffness matrix [Maerten and Maerten, 2006]. 1.3.1.3 Application of boundary conditions: displacement, traction and body forces Boundary conditions are defined to get a consistent restored model and to ensure a unique solution of the differential equations. In restoration, there are two kinds of boundary conditions. The Dirichlet conditions, also called essential boundary conditions, consist in imposing the solution, i.e., a specific displacement, on some boundaries. The Neumann conditions, also called natural boundary conditions, consist in setting a gradient of the displacement, i.e., a positive or negative traction, on some boundaries. The boundary conditions used in mechanics-based restoration are mainly the Dirichlet ones. Traction and gravity conditions are seldom used even if proposed by several authors [Maerten andMaerten, 2006, Titeux andRoyer, 2006].

For a node, a Dirichlet condition fixes at least one component of the nodal displacement. The three components are not necessary all fixed. If a node is not fixed on a specific component, this node is free to move on this component, i.e., the motion on this component is a result of the mechanical resolution. When a nodal component is fixed, the linear or tangent stiffness matrix and the external forces are modified to take into account this Dirichlet condition. Several methods exist such as the penalty method, the Lagrange multiplier method, the augmented Lagrangian method and the perturbed Lagrangian method [e.g., Belytschko et al., 2013]. Muron [2005] uses a penalty method. This technique permits to impose a specific displacement by multiplying the corresponding elements in the matricial system by a high value. Thus, the imposed displacement has an important weight in the system resolution. This method is quite easy to implement but is very controversial. Indeed, the system resolution is sensitive to the choice of the penalty factor [Nour-Omid and Wriggers, 1987].

In the case of the implicit restoration, in which the horizons are defined by isovalues of a scalar field [e.g., Chilès et al., 2004, Cowan et al., 2004, Moyen et al., 2004, Frank et al., 2007, Calcagno et al., 2008, Durand-Riard, 2010, Durand-Riard et al., 2010], two approaches are possible. If the uppermost horizon is explicit, as the boundary conditions on horizons (datuming) are generally applied only on that horizon, so explicit boundary conditions are used as defined above. In that case, the scalar field defining the implicit horizons below follows the deformation of the volumetric mesh. If the uppermost horizon is implicit, implicit boundary conditions are used as proposed by Durand-Riard [2010] and Durand-Riard et al. [2010]. Durand-Riard et al. [2010] gather the tetrahedral nodes which are the closest to the implicit horizon. Then, on each of these nodes, a Dirichlet boundary condition disp is applied to define the restored depth of the implicit horizon:

disp = Z r + d (n) -Z u (n) , (1.19)
with Z r the desired depth of the implicit horizon after restoration, d (n) the distance of the node n to the implicit horizon, and Z u (n) the depth of the node n in the unrestored state. An equivalent condition is used to define the position of points, lines and surfaces along an axis. For more details about this algorithm, see Appendices 3-6 in Durand-Riard et al. [2010]. Durand-Riard et al. [2010] show with an anticline model that implicit restoration provides a very similar restored model as explicit restoration. The use of implicit horizons is a huge progress to reduce meshing constraints and enables the restoration of complex geological models. However, a concern can be pointed out in the use of implicit boundary conditions as defined by Durand-Riard et al. [2010]. Indeed, they define boundary conditions within a volumetric mesh. However, by definition a boundary condition is set on the mesh boundary and the classical finite element method does not consider implicit conditions. More studies and mathematical evidence on this approach should be performed. The extended finite element method, which mathematically defines implicit boundary conditions, may be considered [e.g., Moës et al., 1999, Moës and Belytschko, 2002, Siavelis et al., 2010, 2011, Siavelis, 2011, Siavelis et al., 2013, Belytschko et al., 2013].

1.3.2 Dynamic relaxation method 1.3.2.1 Principle
The dynamic relaxation method [e.g., Day, 1965, Otter, 1965, Underwood, 1979, 1983], also called pseudo-transient analysis, is another finite element method to solve Equation (1.11). Santi et al. [2002] and Muron [2005] developed this method for geomechanical restoration purposes. This method has been used in several studies by Plesch et al. [2007], Guzofski et al. [2009], Durand-Riard [2010], Durand-Riard et al. [2013b] and [START_REF] Stockmeyer | Interplay Between Extension, Salt and Pre-Existing Structure, Offshore Angola[END_REF]. Contrary to the static methods, the dynamic relaxation keeps the temporal terms of Equation (1.11). However, the time is still not taken into account. Temporal terms are kept to transform Equation (1.11) into a recurrence relation. This avoids the matrix system of equations present in the static method. A damping term reduces the transient part of a time dependent simulation in order to quickly reach the static solution which is the objective [e.g., Shizhong, 1988, Muron, 2005]. Thus, the transient steps have no physical meaning and do not matter.

Equation (1.11) is function of the time t. Time is by definition a continuous space. To get a solution at a specific time, it is necessary to convert this continuous space into a discrete space. Let {t 0 , t 1 , . . . , t n-1 , t n , t n+1 , . . . , t ∞ } with n ∈ N denote the time discretization. Let ün , un and u n respectively denote the acceleration, the velocity and the displacement at time t n . Equation (1.11) can be rewritten as:

M • ün + C • un + F int,n -F ext,n = 0.
(1.20) Time integration is used to solve such a problem. Until the end of this chapter, we assume a constant step between two consecutive times: ∆t = t n+1 -t n ∀n ∈ N. The dynamic relaxation uses an explicit time integration. In the finite element method, explicit means that the solution at the time t n only depends on the previous time steps and not on other values of the time t n (implicit time integration). Here explicit/implicit has no link at all with the structural modeling. The explicit time integration in question is the central difference method [Belytschko et al., 2013] which belongs to the finite difference method. For a better precision, half-time integration is used in central difference [e.g., Oakley and Knight, 1995a, Muron, 2005, Belytschko et al., 2013]: .21) In the central difference scheme, it is common to use a diagonalized version of the mass matrix and of the damping matrix. This process is called lumping [Oakley and Knight, 1995a, Zienkiewicz and Taylor, 2000a, Belytschko et al., 2013]. Moreover the damping matrix is often assumed to be linearly proportional to the mass matrix [Muron, 2005, Oakley andKnight, 1995a]:

un+ 1 2 = u n+1 -u n ∆t ; un = un+ 1 2 + un-1 2 2 ; ün = un+ 1 2 -un-1 2 ∆t . ( 1
C = cM , (1.22)
with c a damping factor. Inserting (1.21) and (1.22) into (1.20) it follows [Oakley andKnight, 1995a, Muron, 2005] (Appendix 1.G):

u n+1 = u n + ∆t un+ 1 2 with un+ 1 2 = 2 -c∆t 2 + c∆t un-1 2 + 2∆t 2 + c∆t M -1 • F ext,n -F int,n . (1.23)
This equation has the advantage of being simple to solve since step n+1 can be determined by previous steps. Contrary to the static method (Section 1.3.1), no system of equations needs to be solved [e.g., Muron, 2005]. As the resolution is node by node, no global assembly is required, reducing the computer memory cost. In the linear case, Muron [2005] assembles even so the stiffness matrix once and for all at the beginning of the restoration process to limit the computational time, to the detriment of the computer memory. Furthermore, nonlinearities are easier to handle since no linearization is necessary. Equation (1.23) is always linear even if there are non-linearities such as large deformations or contact conditions [Oakley and Knight, 1995a]. Moreover, as the mass matrix M is often chosen lumped [e.g., Oakley and Knight, 1995a, Zienkiewicz and Taylor, 2000a, Muron, 2005, Belytschko et al., 2013], that is to say diagonalized, computing its inverse is straightforward (equal to the inverse of the diagonal components). Dynamic relaxation, by the direct computation of the next time step and the use of a one-dimensional table to store the data, is a straightforward and parallelizable algorithm [Oakley and Knight, 1995a,b, Oakley et al., 1995, Topping and Khan, 1994]. Parameters are chosen to accelerate the resolution. For instance, as the mass matrix M and the damping matrix C are not physical, they are defined to quickly reach the steadystate [Oakley and Knight, 1995a]. The major drawback of the dynamic relaxation method is that it lacks stability and it can be time-consuming due to the small time step required to reach convergence [Oakley and Knight, 1995a, Muron, 2005, Belytschko et al., 2013]. To select the parameters which will help to a better convergence, the reader can refer to Cassell and Hobbs [1976], Papadrakakis [1981], Oakley and Knight [1995a], and Muron [2005].

Algorithm 2 Dynamic relaxation algorithm. Method developed by Muron [2005] for restoration purposes.

u 0 ← 0 ; u0 ← 0 ; u 1 2 ← ∆t 2 M -1 • F ext,0 ; u 1 ← ∆t u 1 2 ; n ← 1 repeat Compute external forces F ext,n Compute internal forces F int,n un+ 1 2 ← 2-c∆t 2+c∆t un-1 2 + 2∆t 2+c∆t M -1 • F ext,n -F int,n u n+1 ← u n + ∆t un+ 1 2 n ← n + 1 until energy < tolerance 1.3.2.

Definition of the Dirichlet boundary conditions

The definition of the displacement boundary condition consists in a correction of un+ 1 2 in Equation (1.23) [Muron, 2005]. As in the static approach, a displacement condition corresponds to a vector V disp . As the dynamic relaxation is a pseudo-temporal process, this vector is applied incrementally on the node [Belytschko et al., 2013, p. 337]. At the time t n , V disp is recomputed in function of the current geometry and a very small fraction of this vector is applied. As the boundary conditions are applied by very small increments, the dynamic relaxation method is very slow [Belytschko et al., 2013, p. 337]. This portion of the projection vector replaces un+ 1 2 prior to the calculation of un+1 (Equation (1.23)). To be homogeneous to a velocity, V disp is divided by the time step ∆t. Thus the corrected un+ 1 2 corrected is:

un+ 1 2 corrected = f actor × V disp ∆t , (1.24)
with f actor ∈ 0; 1 the factor which defines the small increment of V disp .

Conclusions on the various finite element approaches

Table 1.2 provides a summary of the different finite element methods used in the mechanicsbased restoration. Concerning the choice between a static method and the dynamic relaxation, the main difference between these two methods is on the numerical stability and the computational time. Implicit methods, in the sense of numerical resolution and not structural modeling, are more stable than explicit methods which can easily diverge. Indeed, explicit methods are stable under specific conditions [Oakley andKnight, 1995a, Belytschko et al., 2013]. The dynamic relaxation method is an explicit method and therefore has convergence issues [e.g., Papadrakakis, 1981, Muron, 2005, Moretti and Titeux, 2007]. A small time step is generally required to ensure convergence, but that increases the computational time [Oakley andKnight, 1995a, Belytschko et al., 2013]. The main interest of such an unstable method is that it significantly reduces the computer memory requirements [e.g., Oakley and Knight, 1995a]. Static methods which use a global system of equations for the entire mesh need much computer memory and the matricial resolution may be heavy. However, nowadays, memory is not a major problem in modern computers anymore. For a more detailed comparison between the two methods, the reader can refer to Pica and Hinton [1981].

Handling faults in geomechanical restoration solved by the finite element method

Handling faults in restoration is fundamental. Indeed, their geometry, their mechanical behavior and their evolution through time have a huge impact for instance on reservoir characterization. As we have seen, mechanics-based restoration is based on continuum mechanics. Dealing with faults is peculiar since by nature a fault represents a discontinuity in rock. Although a fault corresponds to a volume in nature, it is generally represented by a surface for simplicity in geomodeling [Fossen, 2016, p. 178]. In restauration, faults correspond to surfaces. Another hypothesis is that fault contact is compliant: between two fault mirrors neither holes nor gaps are wished. In the literature there are two main methods to deal with fault contact. Moretti et al. [2006] and Muron [2005] develop fault contact boundary conditions based only on geometrical assumptions, whereas Maerten and Maerten [2006], Muron [2005] and Tang et al. [2016] use a method based on contact mechanics [Wriggers and Laursen, 2006]. All the fault contact methods use a master/slave approach. However, the definition of master and slave may vary. All the contact methods aim to connect the nodes of the slave fault onto the master fault. There are three kinds of contacts. The first one connects a slave node onto a master node (node to node boundary condition). This contact removes a fault dip slip. The second one connects the nodes of a slave line (continuous set of edges) onto a master line (node to edge boundary condition). This condition is used to tie horizon cutoff lines, generally of the uppermost horizon. The last contact is the connection of the nodes of a slave surface onto a master surface (node to surface boundary condition). This condition is used to avoid gap and penetration between fault mirrors. This fault contact method is based on a slave/master approach. The slave is constrained to move relatively toward the master. The latter remains fixed relatively to the slave [Moretti et al., 2006, Muron, 2005]. The master motion is not constrained by the slave but by other potential boundary conditions and by the mechanical simulation. Thus the restored state is strongly dependent on the choice of the slave and the master [Moretti et al., 2006]. As this approach is purely geometric, there is no notion of frictionless gliding which is a mechanical concept. However, in the case of absence of tangential constraint, such geometric contacts may mimic a frictionless behavior.

1.4.1.2 Geometric fault contact in the static method Lepage et al. [2004] and Moretti et al. [2006] use a purely geometric fault contact method in a static finite element solver. Each node of the slave fault is projected onto the master fault. The point on the master fault is not necessary a mesh node. Let A s denote the slave node to project. Let A m denote the master point which corresponds to the projection of A s on the master fault. A m is a point on a mesh edge or on a mesh facet. The mesh edge or the mesh facet is defined by the nodes A m i (Figures 1.8-1.9) with i ∈ 0; n -1 , n being the number of nodes in the mesh element (e.g., 2 for an edge, 3 for a triangle, 4 for a quadrangle). The position of A m within the mesh element is defined by its barycentric coordinates a i [e.g., Meyer et al., 2002, Moretti et al., 2006] such as: The projection vector P applied on the node A s is the vector linking A s to A m . Using the barycentric coordinate system it comes [Moretti et al., 2006]:

A m = n i a i A m i . ( 1
P = n i a i A m i -A s .
(1.26)

The way to project the slave nodes is a decision which may have important consequences on the restored state [Moretti et al., 2006]. Moretti et al. [2006] advise to use the same curvilinear abscissas on the uppermost horizon cutoff lines to tie the uppermost horizon parts. At the year of their paper, there was no rule to ensure neither gap nor penetration along fault mirrors. A projection, orthogonal to the master surface, may be used as a simple guess [Muron, 2005]. Furthermore, as the master block moves independently of the slave block, the contact process is non-linear (Figure 1.10). Indeed, contact conditions change according to the relative position of the slave to the master through the mechanical processes until convergence. Finally, Moretti et al. [2006] use geometric contact conditions to tie geological layers along horizon interfaces. Indeed in their approach the mesh nodes are duplicated along the horizons. Contrary to contact conditions for faults, they prevent any sliding along horizons, i.e., horizons are "sticking interfaces". 1.4.1.3 Geometric fault contact in the dynamic relaxation method Santi et al. [2003] and Muron [2005] also use a geometric method for the fault contacts. As for the Dirichlet boundary conditions, the method consists in a correction of un+ 1 2 in (1.23) to take into account the contact condition [Hallquist, 1998, Muron, 2005]. As in the static approach, slave nodes are projected onto the master and the projection vector P is defined in a same manner. As for the displacement boundary conditions, the projection vector is not entirely applied on the slave node, but is done incrementally [Belytschko et al., 2013, p. 337]. At the time step n, P is recomputed and a very small fraction of this vector is applied. This portion of the projection vector is added to un+ 1 2 prior to the calculation of un+1 (Equation (1.23)). P is divided by the time step ∆t to be homogeneous to a velocity.

Thus the corrected

un+ 1 2 corrected is:

un+ 1 2 corrected = un+ 1 2 + f actor × P ∆t , (1.27)
with f actor ∈ 0; 1 the factor which defines the small increment of P . This velocity correction enables to project the slave onto the master. As said previously the fault contact is unilateral, i.e., the master is not impacted by the fault contact condition. Muron [2005] extended the velocity correction to the master to enable the sharing of the motion between the master and the slave. He defines a kinematic partitioning factor β ∈ [0; 1], and the slave node projection vector becomes βP instead of P . The master nodes, which in theory are not constrained by the fault contact condition in the geometrical approach, are constrained to move. For each node of the master fault surface, a projection vector is defined. This projection vector is multiplied by 1 -β. This approach of splitting the displacement between the hanging wall and the footwall by a kinematic partitioning factor is similar to the approach of Midland Valley [2017a] to handle fault contacts. In frictionless contact, t m t is null. nm and τ m are respectively the normal and the tangential vectors at the point xm on the master boundary. nm is oriented toward the exterior of the master domain Ω m . These relations are the same for the slave.

Geomechanical fault contact methods

Principle

The second approach of fault contact boundary condition is based on the theory of contact mechanics [Zienkiewicz and Taylor, 2000b, Wriggers and Laursen, 2006, Belytschko et al., 2013] and used in geomechanical restorations by Muron [2005], Maerten and Maerten [2006], Guiton and Zammali [2007] and Tang et al. [2016]. In this formulation, the starting point is the equilibrium of the stresses between both fault blocks, which corresponds to the Newton's third law. In other terms, the stress on a point of the master fault is equal to the opposite stress on the same location on the slave fault [Wriggers and Laursen, 2006, p. 70]: .28) with t m the stress on the point of the master fault, and t s the stress on the same point of the slave fault. The stress on a surface can be decomposed into a normal stress t n and a tangential stress t t (Figure 1.11). Tangential stresses correspond to the forces which restrain the motion along the (fault) surface, i.e., they correspond to frictions. An important assumption in mechanics-based restoration is that faults slide freely, i.e., there is no tangential stress. Under this frictionless assumption, the tangential stresses are null [Wriggers and Laursen, 2006, p. 70]:

t m = -t s , ( 1
t s t = t m t = 0. (1.29)
Furthermore, in contact mechanics, formulation of contacts is stated by the non-penetration of bodies. Indeed, several bodies move due to physical processes, and contact conditions are defined to avoid the penetration of a body into another. The gap between the slave and the master g is composed by a normal gap g n oriented perpendicular to the master, and by a tangential gap g t defined by two coordinates in the tangential plane of the master (Figure 1.12). As a penetration is not physical, the normal gap must be positive or null [Wriggers and Laursen, 2006, p. 60]. In this formulation, the constraint of non-penetration is only on the slave, i.e., the slave cannot penetrate the master. In restoration, to ensure geological consistency between fault mirrors, a no gap condition is defined which has a formulation equivalent to no-penetration. Thus, the goal of fault contact is to reach a null normal gap:

g n = 0.
(1.30)

Figure 1.12: Gap vector decomposition. The gap vector g is defined by the slave node x s and the projection of this node on the master boundary xm . The gap vector is decomposed into a normal component g n and a tangential component g t . In frictionless contact, g t is not explicitly defined.

For restoration purposes, g n must be (close to) null. nm and τ m are respectively the normal and the tangential vectors at the point xm on the master boundary. nm is oriented toward the exterior of the master domain Ω m . Ω s is the slave domain.

As the motion is frictionless, there is no constraint on the tangential gap g t . All the previously defined considerations are integrated into the weak form of the motion equation (1.8) (included in the surface stress term):

δW contact = Γ contact 0 -t n δg n dΓ 0 .
(1.31)

1.4.2.

Resolution

The goal of the fault contact is to reach a null normal gap (Equation (1.30)). Equation (1.31) forces the contact between fault mirrors and such a constraint is function of the normal. As the definition of the normal gap is function of the solid geometry, the contact conditions are necessarily non-linear (Figure 1.13). Indeed, as the tangential gap is not imposed, and will be a consequence of all the boundary conditions and the mechanical settings, after a linear step of mechanical simulation, the "residual" normal gap may not be null (Figure 1.13d). Thus, as for the geometrical contact approach, only an iterative process can step by step reduce the absolute magnitude of the normal gap to reach a null value.

As the other terms of the weak form (1.8), Equation (1.31) is discretized by the finite element method, which defines a contact force. The expression of this force depends on the chosen numerical method. Muron [2005] uses the penalty method and Tang et al. [2016] the augmented Lagrangian method. In the approach of Muron [2005], a Newton-Raphson procedure transforms the non-linear equation of motion which includes the contact and potentially large deformations. Thus, the derivative of the contact force by the displacement is needed. This derivative is then added to the tangent stiffness matrix.

It is important to note that, contrary to the geometrical approach, both master and slave move, i.e., the master is not relatively fixed and the slave does not move relatively to the master. The displacement is bilateral and the magnitude of displacement for each side is function of the strain minimization. The only difference between the master and the slave is that the slave can neither penetrate nor have a gap with the master but the contrary is possible. To avoid the master penetration/gap into the slave, Maerten and Maerten [2006] iteratively alternate the master and the slave definition. Thus, the slave (master) side at a specific step of the non-linear resolution loop becomes the master (slave) side at the next step. After convergence, there is no gap/penetration along the fault.

1.4.3 Some considerations on the fault contacts

1.4.3.1 Master/slave choice
The previous sections developed the two contact methods used in geomechanical restoration to handle faults. Both methods are based on a master/slave approach which is different for both methods. In the case of contact mechanics, the master/slave scheme refers to the penetrability of the slave into the master. The slave cannot penetrate the master but the contrary is possible. In the case of a slave mesh fine enough to fit the geometrical irregularities of the master fault surface, the master should not penetrate too much the slave. In addition, even if there is penetration, this one can be judged acceptable if it is small enough in comparison to the size of the fault. Thus, the choice of the master and of the slave is less critical in contact mechanics than in the geometrical contact approach in which the master/slave definition relates to a relative displacement between fault blocks. Furthermore, other considerations should be taken into account to properly choose the master and the slave. The slave surface is preferentially the surface with a finer mesh, a lower stiffness, a smaller extension or a higher curvature [De Soza, 2015].

Physical consistency of the fault contacts

The fault contact boundary conditions tie fault cutoff lines and avoid any gap/penetration between fault blocks. They just ensure a geometrical consistency in the restored state. As no real temporal method is currently used in mechanics-based restoration, it is not possible to say that these contacts enable a sliding along the faults. In the dynamic relaxation method, the contact conditions are incremental. As a result, the transient part of the resolution looks like a sliding, even if by definition the transient part is not physical in this pseudo-temporal method.

In addition, the geometrical contact approach, by defining Dirichlet conditions on the three space directions, constrains a lot the slave surface. This is in contradiction with the purpose of the geomechanical restoration to have the fault slip as a result of a mechanical process. Thus, the geomechanical contact method, by only constraining the motion along the master normals, allows more freedom in the motion, although both slave and master are constrained. An interesting approach for the geometrical method would be to constrain the motion along a direction and leave free the two other directions.

Furthermore, the mechanical contact approach was initially developed to avoid solid penetrations [Wriggers and Laursen, 2006]. In a mechanical simulation, solids move and deform due to forces, and potentially enter in contact. Contact mechanics is used to physically represent this contact which is a consequence of the motion of the solids by mechanical constraints. In restoration, the contacts are mechanical constraints and not a consequence of solid motion. This point may be philosophical but it could be related to the unphysical boundary conditions in the geomechanical restoration [Lovely et al., 2012].

Finally, contact mechanics assumes a stress equality on each side of a fault [Wriggers and Laursen, 2006, p. 70]. An important question is the geological validity of this statement. In the geometrical contact method, there is no rule which implies this stress equilibrium [Muron, 2005].

Comparison between the different methods

We try to sum-up the similarities and the divergences of the different fault contact methods in Table 1

Conclusions

Since the beginning of the century, geomechanical restoration methods have been developed to overcome issues of traditional restoration methods and to propose a real 3D mechanical approximation of the past structures. Several geomechanical restoration approaches exist, all based on mass and linear momentum conservations, and on elasticity. They use boundary conditions to constrain the paleo-geometries (e.g., datuming, fault contacts) and to ensure the uniqueness of the restored model after a numerical simulation. The divergences between the different methods are on several aspects.

The mesh, i.e., the numerical representation to solve differential equations of mechanics, may be a boundary representation (boundary element and mass-spring methods) or a volumetric mesh composed by volumetric elements (finite element and mass-spring methods).

In the case of a volumetric mesh, the volumetric elements can be conformal to the horizons (explicit modeling) or not (implicit modeling, horizons are defined by isovalues of a scalar field). Three main numerical methods are used to solve a restoration problem. The boundary element method has the advantage of lower meshing constraints but it has difficulty to handle mechanical heterogeneities. The mass-spring method solves a 3D problem by a multitude of 1D problems. Therefore this approach is fast but is not a real 3D resolution. The finite element method is the most largely used and has the main advantage to handle complex geometry and complex mechanical heterogeneities.

The finite element method is declined into various approaches. All of them neglect the temporal part of the motion equation and focus on the steady-state. The static approach completely removes the temporal terms of the motion equation. The local method solves the restoration problem at the node scale whereas the global method solves it globally in the entire mesh. The static method is more stable and faster than the dynamic relaxation method which keeps the temporal terms in order to avoid a matricial system of equations. The main advantage of the dynamic relaxation is the reduction of the computer memory cost.

Concerning the boundary conditions, the main differences concern the fault contacts. Two fault contact methods exist, each based on its own master/slave approach which defines an asymmetrical constraint between the footwall and the hanging wall. The geometrical method imposes that the slave moves toward the master which is relatively fixed. It is interesting to point that some authors extended this method to share the displacement between the master and the slave. The mechanical method is based on the stress equilibrium across the fault. Moreover, there is no imposed relative displacement between the slave and the master. The only constraint is that the slave cannot penetrate the master. This review on mechanics-based restoration was performed to provide insight about all the physical hypotheses and mathematical resolution methods. This is a mandatory step to address the challenges of the geomechanical restoration approach in an effective way. The general conclusions of this thesis (p. 139) detail these challenges.

Appendices 1.A Notations

The notations used in this paper are mainly inspired from those used by Belytschko et al. [2013].

1.A.1 Mathematical symbols

The mathematical symbols presented in this chapter are in Table 1.4.

1.A.2 Einstein notation

Einstein notation consists in implicitly represent a summation by two repeated indexes [Einstein, 1916]:

a i b i = i a i b i .
(1.32)

1.A.3 Tensors

A tensor is a mathematical object defined by its order and its components. Scalars corresponds to tensors of order 0, vectors to tensors of order 1 and matrices to tensors of order 2. Tensors with an order strictly superior to 2 may be defined. Tensor order o corresponds to the number of base vectors which define the tensor. Let e i denote a base vector with i ∈ 1; o . A vector v is defined by v = v i e i , a matrix by M = M ij e i ⊗ e j , a three-order tensor by C = C ijk e i ⊗ e j ⊗ e k . . . with ⊗ the tensor product. In this paper, as in Belytschko et al. [2013], scalar has no particular font and may be lowercase or uppercase, vectors are in bold and lowercase, and matrices and higher order tensors are in bold and uppercase. There are Section 1.B. Elastic stress/strain laws some exceptions which correspond to classical notations found in the literature such as σ for the Cauchy stress tensor (order 2).

1.A.4 Transpose

B = A T with b ij = a ji .
(1.33)

1.A.5 Matrix product C = AB with C ij = A ik B kj .
(1.34)

AB = BA.
(1.35)

1.A.6 Scalar product

The scalar product between two vectors a and b is generally noted as:

a • b = a i b i .
(1.36)

1.A.7 Product between a vector and a matrix

c T = a T B with c i = a j B ji .
(1.37)

1.A.8 Double dot product

A : B = A ij B ij .
(1.38)

1.A.9 Gradient of a vector

F = ∇f with F ij = ∂f i ∂x j .
(1.39) (1.42) is generally written in the Voigt form by simplicity [e.g., Belytschko et al., 2013]:

1.A.10 Gradient of a scalar f = ∇f with f i = ∂f ∂x i . (1.40) 1.A.11 Divergence of a matrix a = ∇ • A =    ∂A 11 ∂x 1 + ∂A 21 ∂x 2 + ∂A 31 ∂x 3 ∂A 12 ∂x 1 + ∂A 22 ∂x 2 + ∂A 32 ∂x 3 ∂A 13 ∂x 1 + ∂A 23 ∂x 2 + ∂A 33 ∂x 3    =    a 1 a 2 a 3    with a j = ∂A ij ∂x i . ( 1
{σ} = [D] • { } (1.43)        σ xx σ yy σ zz σ yz σ xz σ xy        =        C xxxx C xxyy C xxzz C xxyz C xxxz C xxxy C yyxx C yyyy C yyzz C yyyz C yyxz C yyxy C zzxx C zzyy C zzzz C zzyz C zzxz C zzxy C yzxx C yzyy C yzzz C yzyz C yzxz C yzxy C xzxx C xzyy C xzzz C xzyz C xzxz C xzxy C xyxx C xyyy C xyzz C xyyz C xyxz C xyxy        •        xx yy zz 2 yz 2 xz 2 xy       
.

(1.44)

The elasticity tensor in the general (anisotropic) case is function of 21 coefficients [Bower, 2010, p. 77]. Among these coefficients there are Young's moduli E i with i ∈ {x, y, z}, Poisson's ratios ν ij and shear moduli µ ij with (i, j) ∈ {x, y, z} 2 | i = j. E i represents the theoretical constraint along the axis i to reduce by two the size of the solid along this axis. ν ij defines the strain along the axis j due to a strain along the axis i. In the general case ν ij = ν ji . µ ij is defined by the division of the shear stress in the plane (i, j) by the strain in this plane due to a force along the axis i. In the general case µ ij = µ ji .

1.B.2 Stress tensor conversion

In mechanics, several stress tensors exist such as the Cauchy stress tensor σ which is defined on the deformed (restored) domain Ω [Belytschko et al., 2013, Section 3.4]. This stress tensor is symmetric. The equivalent of the Cauchy stress in the initial (unrestored) domain Ω 0 is the nominal stress P . The nominal stress is not symmetric and its transpose is the first Piola-Kirchhoff stress tensor P T . The second Piola-Kirchhoff stress tensor S is often used since it is symmetric. The conversions between these stress tensors are [Belytschko et al., 2013, p. 106]:

P = JF -1 • σ (1.45) S = P • F -T , (1.46)
with F the deformation gradient, J its determinant, F -1 its inverse, and F -T the inverse of the transpose of F .

1.B.3 Neo-Hookean law

Neo-Hookean law is an extension of Hooke's law for large deformations [e.g., Treloar, 1948, Ogden, 1997, Muron, 2005, Belytschko et al., 2013]. It is a non-linear isotropic hyperelastic law which means that the work does not depend on the path of deformation [Belytschko et al., 2013]. The second Piola-Kirchhoff stress S of such a material is defined by a potential energy ψ such as [Belytschko et al., 2013, p. 248]: .47) with C the right Cauchy-Green tensor. According to Malvern [1969] and Belytschko et al. [2013], the potential energy ψ can be expressed by the invariants of C. Let I 1 , I 2 and I 3 denote the three invariants of C. The expression of the invariants and their derivatives are [Carlson and Hoger, 1986, Oakley and Knight, 1995a, Zienkiewicz and Taylor, 2000b, Belytschko et al., Section 1.C. Isotropic materials 2013]:

S = 2 ∂ψ (C) ∂C , ( 1
I 1 (C) = trace (C) ; I 2 (C) = 1 2 trace (C) 2 -trace C 2 ; I 3 (C) = det (C) (1.48) ∂I 1 ∂C = I ; ∂I 2 ∂C = I 1 I -C T ; ∂I 3 ∂C = I 3 C -T . (1.49)
In addition, C is symmetric, so:

C T = F T • F T = F T • F T T = F T • F = C.
(1.50) Hence the expression of the second Piola-Kirchhoff stress S:

S = 2 ∂ψ (C) ∂C = 2 ∂ψ (I 1 , I 2 , I 3 ) ∂C = 2 ∂ψ ∂I 1 ∂I 1 ∂C + ∂ψ ∂I 2 ∂I 2 ∂C + ∂ψ ∂I 3 ∂I 3 ∂C = 2 ∂ψ ∂I 1 I + ∂ψ ∂I 2 I 1 I -C T + ∂ψ ∂I 3 I 3 C -T = 2 ∂ψ ∂I 1 I + ∂ψ ∂I 2 (I 1 I -C) + ∂ψ ∂I 3 I 3 C -1 S = 2 ∂ψ ∂I 1 + I 1 ∂ψ ∂I 2 I -2 ∂ψ ∂I 2 C + 2I 3 ∂ψ ∂I 3 C -1 .
(1.51)

A Neo-Hookean material is defined for the potential energy [Belytschko et al., 2013, p. 252]:

ψ = 1 2 λ (ln J) 2 -µ ln J + 1 2 µ (I 1 -3) , (1.52)
with ln the natural logarithm. In addition, J = det (F ) = det F T and the determinant of the product of two matrices is equal to the product of the determinants. Thus:

I 3 = det (C) = det F T • F = det F T × det (F ) = J 2 .
(1.53) By injecting this relation into (1.52) it comes:

ψ = 1 2 λ ln I 3 2 -µ ln I 3 + 1 2 µ (I 1 -3) . (1.54)
The derivative of ψ by the invariants of C are:

∂ψ ∂I 1 = µ 2 ; ∂ψ ∂I 2 = 0 ; ∂ψ ∂I 3 = λ ln J -µ 2J 2 .
(1.55)

Finally [Belytschko et al., 2013, p. 252]:

S = λ ln (J) C -1 + µ I -C -1 .
(1.56)

1.C Isotropic materials 1.C.1 Elastic parameter simplification in the isotropic case

In an isotropic material, the elastic behavior is the same whatever the direction, hence:

∀i ∈ {x, y, z} , E i = E ; ∀ (i, j) ∈ {x, y, z} 2 | i = j, ν ij = ν. ( 1

1.C.2 Conversion of the isotropic elastic parameters

The conversion from Young's modulus E and Poisson's ratio ν to Lamé's parameters λ µ is [Sokolnikoff, 1956, p. 68]:

λ = Eν (1 + ν) (1 -2ν) ; µ = E 2 (1 + ν)
.

(1.58)

The second Lamé parameter µ is the shear modulus. The conversion from Young's modulus E and Poisson's ratio ν to bulk modulus K is [Sokolnikoff, 1956, p. 70]:

K = E 3 (1 -2ν)
.

(1.59)

1.C.3 Isotropic elasticity tensor

For a linear material, the Voigt definition of D is:

[D] =        2µ + λ λ λ 0 0 0 λ 2µ + λ λ 0 0 0 λ λ 2µ + λ 0 0 0 0 0 0 µ 0 0 0 0 0 0 µ 0 0 0 0 0 0 µ       
.

(1.60)

1.D Transverse isotropic materials used to reach flexural slip mode

1.D.1 Transverse isotropic materials

A transverse isotropic material has an isotropic elastic behavior along parallel planes and an anisotropic elastic behavior along the orthogonal direction to these planes. We consider here the the anisotropic axis is depth axis z, i.e., x and y define isotropic planes. (Figure 1.14).

x and y are completely equivalent, therefore [Bower, 2010, p. 84]:

E x = E y = E p ; ν xy = ν yx = ν p ; ν xz = ν yz = ν pt ; ν zx = ν zy = ν tp ; µ xz = µ yz = µ t , (1.61)
with p and t which respectively mean plane and transverse. To be consistent with Bower's notation, we note E z = E t , and µ xy = µ p . E p is Young's modulus within the isotropic plane. E t is Young's modulus along the anisotropic axis (orthogonal to the isotropic plane). ν p is Poisson's ratio within the isotropic plane. ν pt represents the strain along the anisotropic axis due to a strain within isotropic planes. ν tp represents the strain along isotropic planes due to a strain along the anisotropic axis. µ p is the shear modulus in the isotropic planes. µ t is the Section 1.D. Transverse isotropic materials used to reach flexural slip mode shear modulus of the planes containing the anisotropic axis. ν pt and ν tp are not equal and are related by:

ν pt E p = ν tp E t .
(1.62)

A transverse isotropic material is defined by the five following elastic parameters: E p , E t , ν p , ν pt (or ν tp using (1.62)), and µ t (µ p is deduced from E p and ν p ). In the case of a transverse isotropic behavior, the elastic tensor in the Voigt form is [Bower, 2010, p. 84]: .63) with:

[D] =         d 11 d 12 d 13 0 0 0 d 12 d 11 d 13 0 0 0 d 13 d 13 d 33 0 0 0 0 0 0 d 44 0 0 0 0 0 0 d 44 0 0 0 0 0 0 (d 11 -d 12 ) 2         , ( 1
d 11 = E p (1 -ν pt ν tp ) Γ ;
(1.64)

d 12 = E p (ν p + ν pt ν tp ) Γ ;
(1.65)

d 33 = E t 1 -ν 2 p Γ ;
(1.66)

d 13 = E p (ν tp + ν p ν tp ) Γ ; (1.67) Γ = 1 1 -ν 2 p -2ν pt ν tp -2ν p ν pt ν tp ; (1.68) d 44 = µ t ; (1.69) (d 11 -d 12 ) 2 = E p 2 (1 + ν p ) = µ p .
(1.70)

1.D.2 Material upscaling to reach flexural slip mode

Salamon [1968] proposes a method to upscale a stack of layers, defined by isotropic or tranverse isotropic elastic parameters, into a unique layer defined by transverse isotropic elastic parameters (Appendix 1.D.1). The five elastic parameters (U P upperscript) of the layer resulting of the upscaling are [Salamon, 1968, Crea et al., 1981]:

ν U P p = i f i ν i p E i p 1-(ν i p) 2 i f i E i p 1-(ν i p) 2
(1.71)

ν U P pt = 1 -ν U P p i f i ν i pt 1 -ν i p (1.72) E U P p = 1 -ν U P p 2 i f i E i p 1 -ν i p 2
(1.73)

E U P t = 1 i f i E i p E i p E i t - 2(ν i pt ) 2 1-ν i p + 2(ν U P pt ) 2 (1-ν U P p )E U P p (1.74) µ U P t = 1 i f i µ i t , (1.75)
with i the i th layer of the stack to upscale, i ∈ 1; number of layers . If all the layers are defined by isotropic elastic parameters, each layer i is defined by a single Young's modulus E i , a single Poisson's ratio ν i and a single shear modulus µ i (derives from (1.58)). Thus, the expressions of the five parameters of the upscaled unit are simplified [Chalon et al., 2004, Titeux, 2009, Durand-Riard, 2010, Durand-Riard et al., 2013a]:

ν U P p = i f i ν i E i 1-(ν i ) 2 i f i E i 1-(ν i ) 2
(1.76)

ν U P pt = 1 -ν U P p i f i ν i 1 -ν i (1.
77)

E U P p = 1 -ν U P p 2 i f i E i 1 -(ν i ) 2
(1.78)

E U P t = 1 i f i E i 1 -2(ν i ) 2 1-ν i + 2(ν U P pt ) 2 (1-ν U P p )E U P p (1.79) µ U P t = 1 i f i µ i
.

(1.80)

1.E From current configuration to initial configuration

The following mathematical development is from Belytschko et al. [2013, p. 126].

1.E.1 Volumetric integral conversion

The conversion of the integral of f defined on the volume Ω to the integral defined on the volume Ω 0 is [Belytschko et al., 2013, p. 84]:

Ω f dΩ = Ω 0 f.J.dΩ 0 with J = det (F ), (1.81)
with F the gradient of the displacement from Ω 0 to Ω.

1.E.2 Mass conversion

Let ρ and m respectively denote the density and the mass of a volumetric domain Ω. These physical quantities are related by:

m = Ω ρ dΩ. (1.82)
The mass conservation means that the mass of a domain Ω remains constant during the evolution of this domain. In other terms, the mass in the unrestored state is equal to the mass in the restored state:

m = Ω ρ dΩ = Ω 0 ρ 0 dΩ 0 = m 0 .
(1.83)

Equation ( 1.83) can be used to convert the density from the current state to the initial state in total Lagrangian formulation. Using (1.81) it comes:

Ω 0 ρ 0 dΩ 0 = Ω ρ dΩ = Ω 0 ρJ dΩ 0 .
(1.84)

Section 1.F. Integration by parts and Gauss theorem Hence:

Ω 0 (ρJ -ρ 0 ) dΩ 0 = 0. (1.85)
This relation is valid for any not null domain, so [Belytschko et al., 2013, p. 51]: ρJ = ρ 0 .

(1.86)

1.E.3 Stress conversion

The divergence of the Cauchy stress tensor σ in the domain Ω can be written in the domain Ω 0 in function of the nominal stress tensor P [Belytschko et al., 2013, p 126]:

(∇ • σ) i = ∂σ ji ∂x j = ∂ (J -1 F jk P ki ) ∂x j = J -1 F jk ∂P ki ∂x j + ∂J -1 F jk ∂x j P ki (1.87) = J -1 F jk ∂P ki ∂x j = J -1 ∂x j ∂X k ∂P ki ∂x j = J -1 ∂P ki ∂X k = J -1 (∇ 0 • P ) i .
(1.88)

Hence:

∇ • σ = J -1 ∇ 0 • P . (1.89)
1.E.4 Equation of motion in the initial space

∇ • σ + ρb = ρü + ωρ u J -1 ∇ 0 • P + J -1 ρ 0 b = J -1 ρ 0 ü + J -1 ωρ 0 u ∇ 0 • P + ρ 0 b = ρ 0 ü + ωρ 0 u. (1.90)
1.F Integration by parts and Gauss theorem

Ω 0 δu • ∇ 0 • P dΩ 0 = Ω 0 (δu) i (∇ 0 • P ) i dΩ 0 = Ω 0 δu i ∂P ji ∂X j dΩ 0 = Ω 0 ∂ (δu i P ji ) ∂X j dΩ 0 - Ω 0 ∂δu i ∂X j P ji dΩ 0 = Ω 0 ∂ (δu) i P T ij ∂X j dΩ 0 - Ω 0 ∂δu i ∂X j P ji dΩ 0 = Ω 0 ∂ δu • P T j ∂X j dΩ 0 - Ω 0 (∇ 0 (δu)) T ji P ji dΩ 0 Ω 0 δu • ∇ 0 • P dΩ 0 = Ω 0 ∇ 0 • δu • P T dΩ 0 - Ω 0
(∇ 0 (δu)) T : P dΩ 0 .

(1.91) ∇ 0 • (δu • P ) corresponds to the divergence of a vector. The Gauss theorem relates the integral on a volume of the divergence of a vector to the integral on the closed surface boundary of the volume of this vector such as:

Ω 0 ∇ 0 • δu • P T dΩ 0 = Γ 0 δu • P T dΓ 0 = Γ 0 δu • P T • n 0 dΓ 0 , (1.92)
with dΓ 0 = n 0 dΓ 0 , n 0 being the surface outgoing normal. By definition, the stress vector on a surface expressed in the initial configuration t 0 is related to the nominal stress by: t 0 = P T • n 0 . Thus, the final expression of the integration by parts using Gauss theorem is:

Ω 0 δu • ∇ 0 • P dΩ 0 = Ω 0 δu • t 0 dΩ 0 - Ω 0 (∇ 0 (δu)) T : P dΩ 0 .
(1.93)

1.G Dynamic relaxation formula demonstration

By injecting the expressions of un , ün and C (Equations (1.21) and (1.22)) in (1.20) it comes [e.g., Oakley and Knight, 1995a, Muron, 2005, Belytschko et al., 2013]:

M • un+ 1 2 -un-1 2 ∆t + cM • un+ 1 2 + un-1 2 2 + F int,n = F ext,n 2 + c∆t 2∆t M • un+ 1 2 = 2 -c∆t 2∆t M • un-1 2 + F ext,n -F int,n un+ 1 2 = 2 -c∆t 2 + c∆t M -1 • M • un-1 2 + 2∆t 2 + c∆t M -1 • F ext,n -F int,n un+ 1 2 = 2 -c∆t 2 + c∆t un-1 2 + 2∆t 2 + c∆t M -1 • F ext,n -F int,n . (1.94)
By incoporating this result into the expression of un+ 1 2 (Equation (1.21)) it comes Equation (1.23).

Introduction

Structural restoration is a valuable tool to investigate the geometries of geological structures through time, assess the validity of structural interpretations, and analyze strain fields [e.g., Chamberlin, 1910, Dahlstrom, 1969, Gratier et al., 1991, Léger et al., 1997, Williams et al., 1997, Rouby et al., 2002, Griffiths et al., 2002, Dunbar and Cook, 2003, Muron, 2005, Groshong, 2006, Maerten and Maerten, 2006, Moretti, 2008, Durand-Riard et al., 2010, 2013b, Maerten and Maerten, 2015, Vidal-Royo et al., 2015, Stockmeyer et al., 2017]. Over more than fifteen years, geomechanical restoration approaches that approximate natural rock behavior have been developed to overcome several limitations of traditional geometrical restoration methods [e.g., Fletcher and Pollard, 1999, Maerten and Maerten, 2001, Santi et al., 2003, Muron, 2005, Moretti et al., 2006, Maerten and Maerten, 2006, Moretti, 2008, Guzofski et al., 2009, Durand-Riard et al., 2010, Lovely et al., 2012, Vidal-Royo et al., 2012, Durand-Riard et al., 2013b, Maerten and Maerten, 2015, Vidal-Royo et al., 2015, Tang et al., 2016]. Mechanics-based restoration follows the fundamental physical laws of continuum mechanics, i.e., mass and linear momentum conservations, and invokes a linear or non-linear elastic constitutive relation to govern rock deformation. Boundary conditions are required to unfault and unfold geological structures simultaneously in the simplest manner possible and to obtain a unique solution. Three general types of boundary conditions have been shown to yield geologically reasonable results: (1) an imposed displacement to flatten the uppermost horizon, (2) a set of contacts to ensure fault compliance (neither gap nor penetration between fault blocks, and contact of the uppermost horizon fault cutoff lines), and (3) the definition of pin walls, pin lines and pin nodes to fix degrees of freedom and guarantee that the solution is unique [e.g., Plesch et al., 2007, Guzofski et al., 2009, Vidal-Royo et al., 2012, Durand-Riard et al., 2013a,b, Stockmeyer et al., 2017].

Boundary conditions have an important impact on the restored geometry. While several studies have shown that simple boundary conditions can yield viable restoration results [e.g., Maerten andMaerten, 2006, Guzofski et al., 2009], there are many pitfalls. For example, Lovely et al. [2012] show a simple example in which classical boundary conditions applied to a geomechanical restoration lead to unphysical strain fields and that a different set of boundary conditions significantly changes the resultant strain field (see Figure 1 in Lovely et al. [2012]). This uncertainty in appropriate boundary conditions is particularly problematic if one intends to analyze the corresponding stress or strain for fracture analysis or other purposes [Maerten and Maerten, 2006, Mejía-Herrera et al., 2014, Stockmeyer et al., 2017]. The problems illustrated by Lovely et al. [2012] are: (1) there may be instances when the classical boundary conditions, as defined above, may be unphysical and ( 2) there is no specific guideline to choose appropriate boundary conditions in restoration. Durand-Riard [2010], Lovely et al. [2012] and Durand-Riard et al. [2013b] suggest that these classical boundary conditions may be insufficient to restore geologically consistent and physical strain. They show on synthetic models that a lateral displacement boundary condition along a boundary wall is necessary to recover the expected strain in compressive, extensional, and strike-slip and oblique-slip contexts. In addition, they show that the restoration displacement field is only consistent with the forward displacement field when the amount of the displacement condition on a wall is equal to the amount of forward displacement. These works highlight the need for additional constraints derived from geologic or tectonic insights for mechanics-based restoration. The main challenge for defining these additional constraints is that they require knowledge of the deformation history, which is rarely accessible and, ideally, should be an output of the mechanics-based restoration. Moreover, these studies of boundary conditions were applied to numerical or synthetic models, which are typically idealizations of natural geologic structures, and present additional uncertainties and assumptions (structural interpretation, deformation path, etc.). Models from laboratories are often used to validate restoration methods [e.g., Schultz-Ela, 1992, Yamada and McClay, 2003, Maerten and Maerten, 2006, Groshong et al., 2012, Moretti and Callot, 2012]. As they are laboratory experiments, the forward boundary conditions and the mechanical behaviors are known, and the kinematic evolution of structures may be recorded. Moreover, the interpretation uncertainties of these deformed forward models are generally small, such that the applied boundary conditions can be considered the primary source of uncertainty in restoration attempts. This is a significant benefit compared to restoration attempts of natural structures, where restoration uncertainties result from the interplay of boundary condition uncertainties and of structural interpretation uncertainties [Gratier et al., 1991, Schultz-Ela, 1992].

In this study, we performed a sequential restoration on an analog model deformed in laboratory. Computed tomography (CT) images capture the sequential development of this analog model. These images constrain the forward deformation path of each structure, cap-turing paleo-geometries through time on one edge of the analog model, and, thus, provide a reasonable benchmark for restoration quality and boundary condition testing. In the following, we describe the structural sandbox model and our tests of boundary conditions with the goal of restoring deformed geometries and related fault slip that are consistent with the reference model paleo-geometries. New fault compliance boundary conditions are proposed to handle the complex fault network identified on the CT images. In addition, we propose methods to define lateral displacement boundary conditions without detailed knowledge of the forward deformation path, improving the viability of the 3D geomechanical restoration method for use with natural geologic structures.

Case study and its representativeness 2.1.1 Extensional sandbox model: supra-salt structures

Geologists typically have access only to the present-day state of deformed strata, often times informed by sparse and uncertain data. Such data may be consistent with multiple interpretations that may vary significantly [e.g., [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF], Bond et al., 2007, Wellmann et al., 2010, Bond, 2015, Cherpeau and Caumon, 2015]. Thus, the analysis of rock deformation through time is made difficult by the lack of direct information on paleo-structures and the limitations of the available data. To overcome some of these concerns, laboratory analog models are widely used to model viable deformation styles and paths of natural geologic structures [e.g., McClay, 1990]. They provide a way to follow the evolution of well-known geometries under known physical mechanisms through time. X-ray tomography is a common, non-destructive method to image the 3D structures of a deformed structural sandbox [e.g., Colletta et al., 1991, Callot et al., 2012]. X-ray tomography resolution on an analog model is generally sufficient to study deformed structures with minimal geometric uncertainties. Distinct horizons and faults can be observed due to density contrasts in the model's stratigraphy. Thus, as pointed by Colletta et al. [1991], X-ray tomography is a valuable tool to analyze the temporal evolution of laboratory models. Moreover, an analog model must be properly defined to reproduce the behavior of geological structures. Scale, mechanical materials, physical processes, and timing are examples of parameters to consider for the purpose of assessing the degree to which analog models represent natural structures.

In this paper, we restore a laboratory model analogous to supra-salt extensional structures observed in salt basins around the world, such as the Gulf of Mexico, Angola and Morocco. It is well established that dry sand (no cohesion) is a viable material for modeling brittle and ductile rock deformation in sedimentary systems [e.g., Panien et al., 2006, Victor and Moretti, 2006, Dooley et al., 2007, Callot et al., 2012, Moretti and Callot, 2012, Darnault et al., 2016]. Moreover, Weijermars et al. [1993], Victor and Moretti [2006] and Moretti and Callot [2012], among others, have shown that an analog composed by a stack of sand above silicone can produce structures representative of natural salt basins. Silicone has a very weak rheology relative to sand. Thus, sand layers deform and may penetrate into the silicone. This effectively reproduces the subsurface at the interface between a viscous salt layer and overlying brittle rocks [Weijermars et al., 1993].

Our work is based on a deformed structural sandbox done in laboratory by IFPEN (http: //www.ifpenergiesnouvelles.fr) and C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System (http://www.ccreservoirs.com) (Figure 2.1) to reproduce extensional salt structures. The model box was initially composed of two horizontal layers composing a pregrowth stratigraphy: one of silicone at the bottom with a thickness of 1.8 cm (0.71 in), and one of sand above with a thickness of 4 mm (0.2 in). The initial thickness (along Z direction) of the pregrowth strata is 2.2 cm (0.87 in). Along the Y axis the structural sandbox length is 10 cm (3.9 in), and 18 cm (7.1 in) along the X axis. The model box was inclined by 1.5 • [Weijermars et al., 1993, Victor andMoretti, 2006]. Deformation was induced initially by gravity sliding along this tilt (toward the eastern side on Figure 2.1). On the down-dip end of the model there was no wall to restrain the motion of the materials. As the model deformed, alternating layers of pyrex or sand were deposited (one layer every 16 min in mean), further driving deformation by a combination of gravity spreading and gravity gliding [Victor and Moretti, 2006]. This deposition of successive stratigraphic horizons during deformation represents syn-tectonic strata (i.e., growth strata). As pyrex and sand strata are deposited above the silicone, this experiment describes supra-salt structures. At each depositional time step, the newly deposited sediments filled the available model space. In total, 12 layers were deposited (Figure 2.1) over the course of the forward analog model. The total duration of the experiment is 4h16min. The properties of the silicone, pyrex and sand are provided in Table 2.1. The scaling from the analog model scale to real field scale is approximately 1 cm (0.4 in) for 1 km (0.6 mi), consistently with similar analog models [Ellis and McClay, 1988, Dooley et al., 2007, Wu et al., 2009, Hidayah, 2010, Darnault et al., 2016]. See Hubbert [1937] and Ramberg [1981] for more details about the methods used to define this scaling. It is possible to distinguish silicone, sand, pyrex and the fault offsets by tomographic imaging due to their density contrasts. Indeed, faults are visible in the analog model due to sand and pyrex dilatation (areas of lower density) [Colletta et al., 1991, Cobbold and Castro, 1999, Le Guerroué and Cobbold, 2006, Groshong et al., 2012]. Moreover, sand and pyrex have a sufficient Hounsfield density contrast (Table 2.1) to distinguish them in X-ray tomography [Panien et al., 2006, Darnault et al., 2016], allowing the analysis of fault offset. In addition, we assume that friction is negligible on the edges of the structural sandbox model [Souloumiac et al., 2012].

Interpretation of the structural sandbox model

A 3D X-ray tomography volume of the final state of the deformed box was produced. The volume is defined by an X-ray tomography section every 3 mm (0.1 in) along the Y axis, and an X-ray tomography section every 0.6 mm (0.02 in) along the X and Z axes, producing a tomography volume that can be interpreted with similar methods as a 3D seismic reflection survey (Figure 2.2A). Unfortunately, we had access to only a part of the structural sandbox volume. Indeed CT imaging only recorded through time a specific interval of the structural sandbox. Beyond this interval, down-dip, the analog model continued to deform but was out of the scope of the CT imaging. Although we analyzed the majority of the volume, a part on the eastern side could not be considered in our study. From our interpretation, we built a 3D explicit numerical structural model (boundary representation, Figure 2.2B) using SKUA-GOCAD [Paradigm, 2015], in which horizon and fault surfaces are conformal [Caumon et al., 2004]. The geological model is composed of three primary structures: two grabens in the western and central regions of the model, and a series of west-dipping half-grabens in the eastern region of the model. The layer of silicone representing autochthonous salt is not explicitly represented in our numerical model. Within the analog model, we identified 52 faults. As the purpose of the modeling is the restoration and the strain analysis, faults with small offsets were ignored for simplicity (Figure 2.3B). Additionally, a few faults that define narrow fault blocks were removed and the horizons were made continuous (Figure 2.3C).

As discussed by Vidal-Royo et al. [2012] and Pellerin et al. [2014], this eases the meshing and avoids numerous small volumetric elements in the 3D mesh used in restoration, which reduces the computational time of restorations. The final numerical model is composed of 8 horizons and 22 normal faults. The uppermost layer in the structural model corresponds to the uppermost layer of sand in the experiment. Below, each layer in the model represents a layer of pyrex and the underlying layer of sand. The six pyrex basal horizons were not modeled for simplicity. Figure 2.2A presents the model stratigraphy on the north CT image.

Structural uncertainties

Although the analog structures are well imaged, uncertainties in our structural interpretations exist. This is largely the result of approximating diffuse horizons and faults in the analog model with discrete surfaces in our structural representation. Thus, quantifying our interpretation precision is necessary in order to properly evaluate the quality of subsequent restoration results. Figure 2.4 illustrates four examples of interpretation uncertainties of the deformed analog model. Boundaries between white layers (sand) and black layers (pyrex or silicone) are typically blurred gray (Figure 2.4). The thickness of these gray transition zones provides an estimate of the uncertainty associated with an interpreted horizon between two strata intervals. Although this thickness may vary laterally (Figure 2.4), we estimate an error of 0.1 cm (0.04 in) is a representative uncertainty for all of our interpretations. We will be mindful of this precision as we analyze our restoration results. We also note that other restoration uncertainties, such as finite element approximation or mechanical simplifications, although present, are not considered in our uncertainty analysis.

Restoration settings 2.2.1 Physical volumetric model

We created a 3D mesh (Figure 2.5A) from the structural model (Figure 2.2B) using the Geogram [Lévy, 2015], RINGMesh [Botella et al., 2016, Pellerin et al., 2017], VorteXLib [Botella, 2016a,b] and TetGen [Si, 2015a,b] libraries. It is composed of 647,558 tetrahedra and the average tetrahedron length is 0.25 cm (0.098 in). We made some efforts to reduce the number of tetrahedra, which impacts the restoration computational time, and to avoid imprecision due to a coarse mesh. The VorteXLib library enabled us to develop a 3D mesh maximizing the quality of the tetrahedra (equilaterality) to avoid numerical issues during restorations [Parthasarathy et al., 1994, Shewchuk, 2002, Munson, 2007]. The silicone layer was not represented in the model used for the restoration for two primary reasons. First, its rheology is far weaker than the sand and the pyrex, and thus is not considered to contribute any significant resistance. Second, the viscous behavior of the silicone interval cannot be properly represented by the elastic constitutive law invoked in our restoration method. Thus, we focused on restoration of the sand and pyrex layers that overlie the silicone, with the base of the model being a free surface that represents the top of the silicone [START_REF] Stockmeyer | Interplay Between Extension, Salt and Pre-Existing Structure, Offshore Angola[END_REF]. As sand and pyrex are rheologically similar [Panien et al., 2006], we applied homogeneous elastic properties for the entire model: Young's modulus was set to 70 GPa (10 7 psi) and Poisson's ratio to 0.2 [Holtzman et al., 2009].

Classical boundary conditions

A video of one edge of the analog model was recorded by X-ray computed tomography, allowing us to visualize the deformation and model geometries through time. The deformation front was located on the eastern side of the model. In contrast, the western side was only weakly deformed (Figure 2.2A). Thus, we fixed the western wall in the X and Y directions during the restoration, allowing it to only move vertically (Figure 2.5B). As the experiment is inside a box, no flow occurred in the north-south direction (Y axis) through the northern and southern walls. Therefore, during the restoration, we fixed the northern and the southern walls in Y (Figure 2.5C), as recommended by Durand-Riard [2010] and Durand-Riard et al. [2013b] in other deformation contexts. At each restoration step, we set a datum boundary condition for the uppermost stratigraphic surface because we know the original depositional gradient (Figure 2.5D). As the model had a tilt of 1.5 • toward east, we rotated the entire model before each restoration. This allowed us to set our datuming boundary condition to a constant Z value (Figure 2.5D). We rotated our restored models to their proper geometries for proper comparisons between the numerical models and the CT images. The basal horizon, which defines the interface between sand and silicone, i.e., H1, was defined as a free surface [e.g., [START_REF] Stockmeyer | Interplay Between Extension, Salt and Pre-Existing Structure, Offshore Angola[END_REF].

We also defined fault contact conditions to tie the hanging wall and footwall cutoff lines of the uppermost (flattened) horizon (Figure 2.5E) and to avoid any gap or penetration along fault surfaces (Figure 2.5F). We ensured fault compliance by contact mechanics [Wriggers and Laursen, 2006] which is a master-slave approach adapted to restoration purposes by Muron [2005], and Maerten and Maerten [2006]. This method enables us to tie fault blocks without any friction along fault planes [Muron, 2005, Maerten and Maerten, 2006, Wriggers and Laursen, 2006]. The slave surface cannot penetrate nor have a gap with the master surface but the contrary is possible when faults are curved, owing to limited mesh resolution. No relative displacement constraint between the master and the slave is defined: motion is bilateral and is a consequence of both energy minimization and the constraint for the two sides of the fault to be in contact. In the case of the contact of the fault cutoff lines of the restored (uppermost) horizon (Figure 2.5E), the throw is already defined by the datuming condition applied on the Z-component of the uppermost horizon (Figure 2.5D). The heave is defined by contact mechanics as explained above.

Non-classical boundary condition: imposed shortening condition

As extension clearly occurred during forward deformation, we test an optional lateral shortening condition applied to the down-dip model boundary during restoration. We consider this boundary condition analogous to those suggested by Durand-Riard [2010], Lovely et al. [2012] and Durand-Riard et al. [2013b]. This optional shortening condition is limited to motion along the X axis (i.e., in the west-east direction) and is applied to the eastern wall (Figure 2.5G). In this paper, when we refer to a "no shortening condition", we refer to a restoration scenario without this shortening condition set to the down-dip wall. In these scenarios, the eastern wall is free to move along any direction and the resultant shortening is the output of the restoration, ultimately controlled by the datum and fault slip conditions.

Non-classical boundary conditions: contacts between faults

The complexity of a model increases substantially with the numbers of faults due to the increasing number of the interactions between them and with the horizons (i.e., cutoff relationships, [Pellerin et al., 2015]). The model that we restored presents numerous connections between faults (Figures 2.2B complex fault network is a difficult task during structural modeling, but also during each step of a sequential restoration. To accomplish this task, we present in this paper two additional fault contact boundary conditions that we applied in our restorations. These conditions use contact mechanics, as classical fault contact boundary conditions that we previously defined [Wriggers and Laursen, 2006].

Handling branching faults

In our structural modeling procedure, a fault is represented by two surfaces, one for the hanging wall and one for the footwall. This and the fault contact conditions enable the sliding of the fault blocks along the faults (Figure 2.6). As a result, a branch line between two faults is represented on the main fault by two surface internal borders (black dots in (B1) and (C1) in Figure 2.7 with F1 -the main fault). In the case of branching contacts between faults, discontinuities may occur in the restored state if care is not taken (Figure 2.7B). Therefore, we set contact conditions to tie the internal surface borders and thus avoid internal gaps or overlaps in the restored state (Figure 2.7C).

Handling offset fault surfaces

There were several situations where a fault surface was offset by a different fault. To properly characterize these faulted faults, we split each offset fault into two or more distinct fault surfaces. For example, Figure 2.8 shows that fault F1 is cut and displaced by fault F2. In this case, F1 was represented by two independent faults: F1-hw and F1-fw where the labels -hw and -fw respectively refer to the hanging wall and footwall sides of F2. As a result, F1hw and F1-fw are able to move independently. However, as F1 was originally a single fault surface, and since we know all of the faults are normal faults at all times during the forward model, we know that the slip between F1-hw and F1-fw along F2 should decrease through the restoration until it becomes null. In other words, the distance between F1-hw and F1-fw should decrease along F2 until they merge. Upon removing of all the fault offsets, F1-hw and F1-fw should no longer behave independently, but form a single continuous fault surface. We ensure this condition by a set of contact conditions that aims to tie the different connected components of an offset fault (Figure 2.8). In our case, we were able to quickly determine which restoration steps to apply this contact condition for a particular offset fault (presenting apparent continuity) by investigating the CT images that recorded the forward deformation process.

2.3 Results: restoration of the analog model

Sequential restoration

We performed a partial sequential restoration using RINGMecha [Chauvin and Mazuyer, 2016], a mechanics-based restoration library based on the work of Muron [2005] and Durand-Riard [2010]. We used a time-independent finite element solver to perform the restoration [e.g., Zienkiewicz andTaylor, 2000a,b, Belytschko et al., 2013] with a small deformation assumption. After each restoration step, we removed the uppermost, restored layer before performing the subsequent restoration step. Using the classical and newly defined boundary conditions described above, we performed four steps of sequential restoration for our model (Figures . For all the restoration steps except the fourth one, the majority of the dip slip deltas is within the uncertainty range (-0.2 cm and 0.2 cm).

Validation: quantitative comparison with a reference solution

The visual comparison of the restored models with the references provides valuable insight on the quality of our restorations. However, quantitative analysis is necessary to rigorously and objectively assess the restoration quality [Lingrey andVidal-Royo, 2015, 2016] and uncertainties. To quantify the difference between the restored models and the reference paleogeometries, we measured the magnitude of residual dip slip along faults after a restoration step. Dip slip provides a quantitative measure of the recovered strain along each fault. For each fault that crosses the northern edge, we measured the amount of offset for each horizon on the CT pictures. In addition, we computed the dip slip values on each fault in the numerical models at the considered time steps. The difference of dip slip values between the restored model and the reference, that we call delta, is calculated as

delta = DS res -DS ref , (2.1)
with DS res and DS ref respectively the dip slip in the restored state and the dip slip on the reference CT image. The corresponding distributions for the first four restoration steps (with shortening boundary conditions) are shown in Figure 2.13. To avoid bias in this analysis, these distributions do not include the dip slip measures at the uppermost horizon, as these dip slip values are defined by input boundary conditions (Figure 2.5E). Table 2.2 presents the mean and median values for each dip slip delta distribution shown in Figure 2.13, as well as the percentage of dip slip deltas within the picking uncertainty range estimated to be between -0.2 cm (-0.08 in) and +0.2 cm (+0.08 in). This uncertainty value originates from the picking uncertainty (0.1 cm) applied on the footwall and the hanging wall. It also considers similar uncertainties in our interpretations of paleo-geometries and dip slip magnitudes on the CT images. The distributions show maxima near zero delta. In addition, the majority of the residual dip slip measurements are within the uncertainty range considered (Table 2.2). Therefore, we suggest these restorations, which each included the applied shortening boundary condition to the down-dip model wall, are valid. Nevertheless, some slip measurements from restoration models differ significantly from the reference solution. These are clearly not a common result, except in the fourth restoration step (white intervals in Figure 2.13, Table 2 For each restoration step, the mean and the median of the dip slip delta distribution (Figure 2.13) are indicated, in addition to the percentage of dip slip deltas within the picking uncertainty range.

Estimation of shortening

As shown previously, the amount of extension that is restored without the applied boundary condition to the down-dip model wall consistently underestimates the actual amount of extension that occurred in the forward model. For natural structures, the total amount of extension (or shortening) that occurred to yield the present-day geometry is generally unknown. In these cases, an estimation of the amount of displacement can be attempted using 2D kinematic restoration approaches [e.g., Chamberlin, 1910, Dahlstrom, 1969] or 2D areadepth analysis [START_REF] Epard | Excess area and depth to detachment[END_REF], Groshong et al., 2003, Groshong, 2006, Groshong et al., 2012]. Specific markers, such as channel offsets, can be used if present as proposed by Durand-Riard et al. [2013b].

Methods based on rigid motion and bed length conservation

Table 2.3 presents the incremental shortening evaluated by different methods, in particular fault heave and bed length conservation. The former corresponds to the required horizontal displacement to tie the uppermost horizon parts as a pure rigid motion of the fault blocks. The latter, in addition to joining the uppermost horizon parts, assumes that this horizon conserves its bed length and is restored to horizontal. In this case, the horizontal displacement is equal to the horizontal extension of the analog model in the X direction (i.e., down-dip direction) before restoration minus the sum of the lengths of the uppermost horizon parts. For each of these two methods, we used the unrestored model (with applied shortening boundary condition) geometry at each restoration step. Both methods provide displacement estimates that are significantly less than the expected values (Table 2.3). In other words, rigid motion along faults is not an accurate measure of total tectonic displacement for our model and bed lengths did not remain constant through deformation. This latter conclusion is a known expectation for extensional structures [e.g., Xiao and Suppe, 1992]. There is internal deformation accommodated by structures below image resolution or by deformation of a more continuous nature.

Area-depth method

We applied the area-depth method [START_REF] Epard | Excess area and depth to detachment[END_REF], Groshong et al., 2003, Groshong, 2006, Groshong et al., 2012] to estimate the total forward extension without the need of a reference paleo-geometry. The area-depth method may be used to calculate the magnitude of shortening or extension of a system above a basal detachment. A benefit of this method is that it accounts for the displacement due to faults omitted from the interpretation or tectonic strain that is below our imaging resolution, which has been found to accommodate up to 60% of total extension within a given system [e.g., Kautz and Sclater, 1988 Allmendinger , 1992, Baxter, 1998, Groshong et al., 2003]. The area-depth method is independent of the mechanical processes and is based on assumptions of area conservation and plane-strain, given that a thin detachment level exists. The area-depth method defines for each horizon a regional depth of detachment and a lost area inside the graben (below the regional datum and above the horizon), as shown in Figure 2.14. This lost area is equal to the product of the displacement that produced the graben and the depth to the detachment level. It follows that the total extension is given by the lost area divided by the depth to the detachment. We did not plot an area-depth graph of the entire growth sequence, which would integrate each lost area and each distance from the regional to a reference level, since by definition the layers did not undergo the same magnitude of extension [Groshong et al., 2003]. Indeed, such a plot enables to evaluate the common displacement and the depth to detachment only for pregrowth strata or for no-growth sequences of growth strata [Groshong, 2015]. Thus, we assume that the depth to the detachment is known. In our analog model, the definition of the detachment level is not straightforward, as the silicone layer is thick and may act as a distributed detachment zone. Assuming that no slip occurs along silicone boundaries [Weijermars et al., 1993], we approximated the detachment level to be at the middle of the silicone layer (Figure 2.14). The regional level of each horizon is defined by a straight line dipping 1.5 • (parallel to the detachment) and starting from the intersection between the horizon and the most western fault (Figure 2.14). Our calculations only use the north CT image of the analog final deformation stage. The estimates of the shortening magnitude increments for the first four horizons using the area-depth method and the CT images are given in Table 2.3. The amounts of displacement predicted by the area-depth method are within 15% of the shortening magnitudes provided by the CT images. We consider this a valid estimate given the structural uncertainties.

3D dilatation analysis

We propose a complementary approach to evaluate the model forward extension from calculations of dilatation, where

dilatation = 100 × V r -V u V u , (2.2)
with V r and V u respectively the restored volume and the unrestored volume. As previously discussed, horizontal dilatation is expected during the experiment. Due to the small duration and the scale of the experiment, the strata used in our model are not expected to undergo significant vertical compaction [Schultz-Ela, 1992]. Thus, we expect the volume to increase during forward deformation experiment. A consequence is that the volume should decrease during restoration, resulting in negative dilatation calculations from Equation (2.2). We ran a large number of geomechanical restorations varying the magnitude of shortening imposed as a boundary condition. Figure 2.15 represents the proportion of tetrahedra with a positive dilatation from the unrestored state to the restored state according to magnitudes of the imposed shortening conditions. In this way, we attempt to estimate the magnitude of shortening required to minimize the number of tetrahedra with positive dilatation. Similarly, Durand-Riard [2010], with a contractional model, used a lateral (elongation) displacement to reduce the number of tetrahedra with a negative dilatation. The shortenings in Table 2.3 are displayed in the different graphs of Figure 2.15. As expected, the number of tetrahedra with a positive dilatation decreases when the magnitude of the applied shortening increases. However, in each scenario (Figure 2.15), a plateau of diminishing returns develops with additional applied shortening. The beginning of each plateau, as well as the area-depth estimates, provides a much improved estimate of the shortening magnitude than the other methods investigated above (Table 2.3). While this conclusion is still empirical, we suggest that dilatation may be an effective tool to estimate the magnitude of the lateral displacement boundary condition and to evaluate the validity of the restored state.

Discussions

Reasons for a shortening boundary condition

There are several potential explanations for the requirement of an imposed shortening boundary condition. A first reason is the granular nature of the growth strata. According to Groshong et al. [2003], Yamada and McClay [2003], Le Guerroué and Cobbold [2006], and Moretti and Callot [2012], dilatation is likely to occur in structural sandbox models when granular materials undergo shear. This effect enables faults to develop in unconsolidated materials [e.g., Colletta et al., 1991, Cobbold and Castro, 1999, Le Guerroué and Cobbold, 2006, Groshong et al., 2012]. As deformation progresses, additional shear occurs and more voids develop in the system [Groshong et al., 2003, Le Guerroué andCobbold, 2006]. This disorder is at the origin of an increase of the global volume, and thus must be countered by applied .15: 3D dilatation analysis. Proportion of tetrahedra with a positive dilatation after each restoration step (steps 1 to 4) according to different shortenings (in centimeters). Black dots are data points (restoration simulations). For each restoration, the unrestored state is the restored state at the previous restoration step with the imposed shortening from CT image (and without the uppermost restored layer). Shortenings in Table 2.3 are displayed on each graph. Ref: shortening from the CT image. AD: shortening from the area-depth method. HS: sum of the fault heaves of uppermost horizon fault cutoff lines on the northern wall. BL: shortening from the bed length conservation method. NIS: no imposed shortening on the northern wall (restoration without shortening condition). Section 2.5. Discussions Restoration step 1 2 3 4 Forward dilatation of sand and pyrex strata (%) 3.38 0.73 3.07 3.64 Table 2.4: Incremental forward dilatation of the sand and pyrex strata measured on the CT images. See Equation (2.3) and text for calculation details. The forward dilatation of sand and pyrex strata for each restoration step is positive, testifying an increase of area forward in time.

shortening in the restoration. Table 2.4 provides the forward dilatation v (volumetric strain) of the sand and pyrex strata measured on the CT images for each restoration step using

v = 100 × A d i -A u i A u i , (2.3)
with A u i and A d i respectively the area of the sand and pyrex strata on the CT image at deposition time of the layer i and just before the deposition of the layer just above the layer i. At each step, the forward dilatation is positive, which means the analog model area on the northern edge increased through time. Nevertheless, we have just such an evidence of dilatation on the northern edge thanks to the CT images; this dilatation may or may not be compensated elsewhere within the volume. However, due to the style of deformation of the sandbox model, it is very probable that the forward dilatation observed on the CT images is representative of the entire volume. Rock dilatation may exist in real extensional fields in which sediments contain fluids [e.g., Boerner and Sclater, 1992]. Another important reason which explains the need for the shortening boundary condition is that a part of the fault displacement may not be taken into account. Indeed, all the observed faults are not represented, and there may be faults below tomography resolution. Even if their offsets are small, accumulated fault heaves may represent significant forward extension. This is analogous to nature with faults below seismic resolution, which are not imaged and, thus, unable to be represented at the macro-scale [e.g., Groshong et al., 2003].

Residual amounts of fault dip slip values

Although the distributions in Figure 2.13 are encouraging, the number of inconsistent fault dip slip deltas increases with each successive restoration step. A possible explanation for this is that each residual dip slip on these faults is not corrected between restoration steps to fit the dip slip observed on CT images, leading to the accumulation of errors. Another possible explanation is that some faults are kept within the volumetric model whereas they were not present in the analog model at the time corresponding to the restoration step. As faults behave as sliding surfaces in our volumetric mesh, small artificial slip may be present on these surfaces, leading to local inconsistent shear strain. We kept these faults to avoid rebuilding a new structural model, which can be quite time-intensive for such complex fault networks [Zehner et al., 2015]. Another observation of these distributions suggests that the restorations seemed to have recovered too much dip slip (numerous negative deltas). Since the forward deformation path involves friction on faults, this result may be due to the frictionless contacts in our mechanics-based restoration method [Wriggers and Laursen, 2006].

Mismatches with the area-depth method

As mentioned previously, the area-depth method provides a reasonable estimate of the incremental extension that occurred during the forward model (Table 2.3). However, the estimates are not perfect, in particular for the restoration step 4 (Figures 2.12 and 2.15). Several factors may explain these errors. First, material dilatation observed on the CT images and attested by several authors [Yamada andMcClay, 2003, Le Guerroué andCobbold, 2006] is inconsistent with the area conservation hypothesis underlying the area-depth method. As the upper horizons have accumulated less dilatation than the bottom horizons, the constant area hypothesis deteriorates with each successive restoration step. Second, the silicone layer could migrate laterally and blend with the sand and pyrex, leading to area changes. In the analog model, from the CT image of the restoration step 4 to the CT image representing the first unrestored state (northern edge), we calculated a forward dilatation of the silicone to ∼6.8%. Third, as a part of the analog model on the eastern side was not available for analysis in our CT tomography images, we could not integrate this data in our area-depth computations. Fourth, the definition of the detachment level, even based on several reasonable assumptions, is uncertain. Fifth, as the units of sand and pyrex can penetrate into the silicone, the resulting subsidence modifies the definition of the regional levels. This effect is equivalent to the "floating regional" mentioned by Groshong [2015] for a buckle-style fold above a thick salt unit.

Boundary conditions

This study suggests that for extensional systems, the combination of classical boundary conditions and a new lateral displacement boundary condition along the dominant transport direction (Figure 2.5) may yield consistent restored geometries. In this paper, we also propose the use of novel contact conditions to ensure consistent restoration of complex branching and crossing fault geometries. These new constraints enable effective sequential restoration of four steps of the analog model. Without them, only two steps could have been performed, and quality of these restorations would have been reduced. We believe that the boundary conditions presented in Figure 2.5 can be applied in compressive contexts with an elongation displacement condition instead of the shortening condition. Indeed, Durand-Riard [2010] shows that an elongation condition is necessary to properly restore a fault-bend fold model. An estimation of the elongation may be done using the area-depth method [Groshong et al., 2012]. In case of strike-slip faults, displacement conditions parallel to the strike direction should also be considered, as shown by Durand-Riard et al. [2013b].

In the literature and in this paper, all the boundary conditions correspond to displacement conditions except for the mechanical contact conditions which are a mix between displacement and traction conditions [Muron, 2005, Wriggers and Laursen, 2006, Maerten and Maerten, 2006]. Such displacement conditions may lead to unphysical strain fields [Lovely et al., 2012]. In reality, rock deformation is a consequence of force constraints. Maerten and Maerten [2006] suggest the possibility of employing mechanical boundary conditions that incorporate the far field stress as an additional boundary condition. The main difficulty of this technique would be to know the intensity of the forces to apply [Muron, 2005]. A first start could be to use the determined displacement condition for a model (e.g., Figure 2.15) and convert it to a force: dilatation multiplied by Young's modulus in linear elasticity. In addition, the overburden force is not incorporated in our geomechanical restoration method. As our experiment was gravity-driven, it would be interesting to add an overburden body force to the finite element procedure to analyze its impact on the restored geometries.

Conclusions

The restoration of an analog model, in which the structural uncertainties are limited and paleogeometry is well known, enabled us to define effective boundary conditions that yield optimal restored models using mechanics-based restoration. For extensional structures, a shortening boundary condition was applied to obtain a good fit with reference paleo-geometries. Such a condition may be estimated by the area-depth method. Our experiments suggest that an analysis of the volumetric dilatation can complement the estimate of the shortening boundary condition magnitude. Moreover, to handle complex fault networks, we propose the application of contact conditions on internal fault borders and between fault connected components. Ultimately, the methods developed in this paper, in particular the lateral displacement boundary condition, should lead to improved results if applied to geomechanical restorations of natural structures.

Chapter 3

Comparison between mechanics-based and GeoChron-based restorations. Application to a structural sandbox model

Introduction

Since the beginning of the last century, several restoration methods have been developed to analyze the paleo-geometry of rock units through time [e.g., Chamberlin, 1910, Dahlstrom, 1969]. Classical restoration methods are geometric or kinematic, and mainly used on cross sections or on surfaces representing the horizons [e.g., Gratier, 1988, Williams et al., 1997, Rouby et al., 2002, Griffiths et al., 2002, Gjerde, 2002, Dunbar and Cook, 2003, Groshong, 2006]. More recently, such techniques have been extended to meshed volumes [Massot, 2002, Muron, 2005, Medwedeff et al., 2016]. Recently, Medwedeff et al. [2016] propose a 3D geometric restoration method based on the Geo-Chronological model [Mallet, 2004, Moyen et al., 2004, Moyen, 2005, Mallet, 2014]. The Geo-Chronological (GeoChron) model is a 3D chronostratigraphic theoretical framework [e.g., de Groot et al., 2006, Monsen et al., 2007, Wu and Hale, 2015, Labrunye and Carn, 2015, Karimi and Fomel, 2015] which corresponds to the mathematical formulation of the Wheeler space [Wheeler, 1958].

These restoration methods approach the mechanical behavior of rocks by geometrical approximations. For instance, the classical area or volume conservation is an approximation of mass conservation. Several authors pointed out that real rock deformation cannot be assessed without geomechanical rules [Fletcher andPollard, 1999, Gjerde et al., 2002] and that the classical restoration methods need to incorporate real mechanics [Muron, 2005, Maerten and Maerten, 2006, Guzofski et al., 2009]. 3D mechanics-based restoration methods have been developed since the 2000's to propose a restoration technique integrating rock mechanical behavior [e.g., Maerten and Maerten, 2001, Santi et al., 2002, Muron, 2005, Moretti et al., 2006, Maerten and Maerten, 2006, Plesch et al., 2007, Moretti, 2008, Guzofski et al., 2009, Durand-Riard, 2010, Durand-Riard et al., 2010, 2011, 2013a,b, Chauvin et al., 2017]. However, mechanics-based restoration presents several drawbacks. First, the mechanical rock behavior relies on elasticity which is a limited representation of rock undergoing large deformation over geologic time. The actual rock behaviors (e.g., plasticity or visco-elasticity) are not taken into account [Durand-Riard et al., 2010, Lovely et al., 2012]. Thus, the real mechanics is still simplified [Gerbault et al., 1998]. Second, it is not straightforward to choose the rock mechanical properties, i.e., Young's modulus and Poisson's ratio (or equivalent), within a geological model due to the uncertainties. Third, the boundary conditions are controversial. In mechanics-based restoration, the boundary conditions, such as unfolding and unfaulting, are defined to ensure a consistent restored state. It is well established now that classical boundary conditions are unphysical and may lead to an invalid strain field [Durand-Riard, 2010, Lovely et al., 2012, Durand-Riard et al., 2013b, Chauvin et al., 2017]. Unphysical deformation and so stress field is problematic to determine potential fracture areas [Macé, 2006, Maerten and Maerten, 2006, Mejía-Herrera et al., 2014, Stockmeyer et al., 2017].

It is difficult to objectively determine which restoration method is the most appropriate for a specific case study. Both geometrical and geomechanical restorations have proven their applicability in different case studies. Moreover, in real geological cases, as the paleo-geometries are unknown, assessing the validity of restoration can be difficult. The geologist can only estimate the consistency of a restored model, which is a subjective choice depending on the amount of available data and on their uncertainties. The preferred restoration method is the one which produces the most "consistent" restored model. The consistency here means that no available data invalidates the restored state. New external data can afterward question the initial restoration.

On account of the difficulty to select a restoration method in particular, we propose a comparison of two methods. The purpose of this paper is not to provide a complete comparison between geometric/kinematic methods and geomechanical methods in general. We compare the restoration results provided by a geomechanical approach [Chauvin et al., 2017] and a GeoChron-based method [Medwedeff et al., 2016] on the same structural sandbox model. The paleo-geometries of the sandbox are known through time on a cross section, defining a reference solution for restoration purposes. In addition, as the sandbox deformation was performed in laboratory, the structural uncertainties are relatively low. Thus, it is possible to objectively determine which restoration method provides the restored model which is geometrically the closest to the reference solution. After a theoretical description of both methods, we present a geometrical comparison of the sandbox restored states. Finally we discuss about the impact of mechanical properties in restoration and about the applicability of both methods.

Restoration methods

Mechanics-based restoration

Mechanics-based restoration corresponds to solution of a geomechanical boundary value problem. It relies on continuum mechanics. Both mass and linear momentum are conserved. Generally, these laws are used on a static (time-independent) form [e.g., Maerten and Maerten, 2001, Muron, 2005, Moretti et al., 2006, Maerten and Maerten, 2006, Guzofski et al., 2009]:

∇ • σ + ρb = 0, (3.1)
with σ the Cauchy stress tensor, ∇• the divergence operator, ρ the rock density and b the gravity. Elastic mechanical properties are defined within the geological model to mimic rock behavior. A linear or non-linear law is used to relate the strain to the stress. Elasticity is chosen for reversibility and simplicity reasons [Muron, 2005, Moretti et al., 2006, Maerten and Maerten, 2006]. The simplest elastic rock materials are isotropic and defined by Young's modulus and Poisson's ratio or equivalent [e.g., Sokolnikoff, 1956, Gercek, 2007]. These isotropic parameters can vary within the geological model. Boundary conditions, i.e., displacements and tractions, are set to ensure the geometrical consistency of the restored state. They are generally displacement conditions, typically to flatten and to remove the offsets of the uppermost horizon. These boundary conditions and Equation (3.1) constitute the system of equations corresponding to a mechanics-based restoration problem.

The typical elastic restoration workflow consists in six steps: (1) construction of a structural model (boundary representation), ( 2) construction of a volumetric mesh of the structural model, (3) definition of the mechanical properties inside the volume (constant or varying Young's modulus and Poisson's ratio), ( 4) definition of the boundary conditions, (5) resolution by finite elements of the mechanical problem to get the 3D restoration displacement field, ( 6) application of the computed displacement on the volume to get the restored state. The previous steps aim to reverse the deformations undergone by the uppermost layer. In the case of deformation events prior to the deposition of the uppermost layer, or in case of synsedimentary deformation (growth strata), the paleo-geometries of underlying layers are obtained by removing the uppermost restored layer and restarting restoration process from steps 4 to 6. This process is sequentially applied on all strata from top to bottom.

As the comparison in this paper is based on the restorations of Chauvin et al. [2017], we focus on the specific geomechanical restoration method developed by Muron [2005] and Chauvin et al. [2017]. It is based on a time-independent finite element solver [Zienkiewicz andTaylor, 2000a,b, Belytschko et al., 2013]. The resolution consists in solving Equation (3.1) in each elementary volume which composes the 3D mesh discretizing the geological model. The resolution is global, i.e., the contribution of each tetrahedron is assembled into a system of equations representing the entire model. Chauvin et al. [2017] chose the small deformation assumption and use the linear elasticity relation known as Hooke's law to relate stress and strain [e.g., Ramsay and Huber, 2000]. The global system of equations is expressed by the following matrix equation:

K • u = F ext , (3.2)
where K is the stiffness matrix representing the mechanical behavior of the entire model, u corresponds to the restoration displacement field and F ext denotes the external forces applied on the model. In (3.2) only the displacement field u is unknown. This equation is solved by numerical methods (here the conjugate gradient [e.g., M.R and Stiefel, 1952, Bathe, 2014]).

Sliding between fault blocks is allowed by splitting each fault into two surfaces corresponding to the fault mirrors. To handle fault compliance and to remove the offsets of the uppermost horizon throughout restoration, contact mechanics is used [Wriggers and Laursen, 2006]. The contact method is a master/slave approach which constrains the nodes of the slave fault surface to be on the master fault surface [Muron, 2005, Wriggers and Laursen, 2006, Maerten and Maerten, 2006]. That prevents any penetration or gap of the slave into the master, but the contrary is possible. This contact method is bilateral, i.e., the master and the slave move, and their relative motion depends on the global energy minimization (the master is not fixed relatively to the slave).

GeoChron-based restoration

In this section, we present a recent restoration method developed by Medwedeff et al. [2016] and based on the Geo-Chronological model [Mallet, 2004, Moyen et al., 2004, Moyen, 2005, Mallet, 2014].

The GeoChron model

The Geo-Chronological model is a method to calculate the chronostratigraphic coordinates of a geological model [Mallet, 2004, Moyen et al., 2004, Moyen, 2005, Mallet, 2014]. These coordinates correspond to the 3D mathematical formulation of the Wheeler space. In this space, the vertical axis corresponds to the time, and the horizontal plane defines the "paleogeographic coordinates" [Mallet, 2014, p. 10]. As the vertical axis is time, all the horizons are horizontal in this space. Moreover, as Wheeler space is the deposition space, horizons are undeformed and continuous, i.e., unfaulted.

The GeoChron model defines three chronostratigraphic coordinates, u, v, and t, defining the uvt (deposition) space. u and v are the "paleo-geographic coordinates", and t corresponds to the time coordinate. The t field is first built from the horizons [Mallet, 2014, p. 54]. Indeed, each horizon defines an iso-value constraining the interpolation of the t field. There are two mathematical definitions to generate the u and v fields from the t field, each corresponding to a tectonic style: minimal deformation and flexural slip [Mallet, 2014, p. 53-54 and chap. 2.3]. Minimal deformation style minimizes the strain in all directions and is mathematically defined by [Mallet, 2014, p. 53, 71-72]:

∇ u 1 ; ∇ v 1 ; ∇ u • ∇ v 0 ; ∇ t • ∇ u 0 ; ∇ t • ∇ v 0 , (3.3)
with ∇ the gradient operator. Flexural slip minimizes the strain along bedding interfaces and is mathematically represented by [Mallet, 2014, p. 54, 72-74]:

∇ H u 1 ; ∇ H v 1 ; ∇ H u • ∇ H v 0, (3.4)
with ∇ H the projection of the gradient ∇ orthogonally to the horizon H. These equations are solved by Discrete Smooth Interpolation [Mallet, 1989[Mallet, , 1992[Mallet, , 1997[Mallet, , 2014]]. For more details about the uvt field generation, in particular on other potential constraints, see Mallet [2014, chap. 1.9].

Several applications of the GeoChron model are proposed by Mallet [2014]. The primary application is the generation of a structural model from sparse and uncertain data [Mallet, 2014, Section 1.10]. This technology is developed by Paradigm [2015] in the SKUA R structural modeling workflow. See Mallet [2014, Part 2] for other applications.

A restoration method based on the GeoChron model

The Wheeler space cannot be assimilated to a restored space for several reasons [Mallet, 2014, Medwedeff et al., 2016]. First, the temporal axis in the Wheeler space is intrinsically not equivalent to a depth metric axis. Second, at deposition time of a geological layer, the remaining layers are not necessary unfolded and unfaulted (e.g., growth structures). Finally, as every single chronostratigraphic horizon is horizontal in the Wheeler space, thickness variations along each layer are not preserved. Medwedeff et al. [2016] develop a geometric restoration method which is based on the GeoChron model, and avoids the limitations described previously for restoration purposes. First, Medwedeff et al. [2016] build a structural model with the SKUA R structural modeling workflow. Thus, the input unrestored model of the GeoChron-based restoration method is a meshed volume, here composed by tetrahedra, in which the horizons are represented by isovalues of a scalar field [Frank et al., 2007], here the t field. This volume represents the present-day structural model. Second, within this volume, a uvt field is generated with only the horizon to restore as constraint for the t field. Thus, this new uvt field is different from the one generated to build the complete unrestored stratigraphy. The other horizons are not used to constrain the t field. Additional constraints can be defined on the uvt field to, for example, integrate strike-slip constraints, or to enforce a specific kinematic behavior (e.g., flexural slip or inclined shear). Third, the volume is transformed into the uvt space. To enable the sliding along faults, nodes along fault surfaces are duplicated [Medwedeff et al., 2016]. The compliance between fault blocks is ensured by additional constraints called veclinks [Ait Ettajer, 1995, Moyen, 2004]. The key of the GeoChron-based restoration method is the use of only a single horizon as constraint of the t field. Indeed, the subsequent t field gradient is almost constant [Medwedeff et al., 2016]. Thus, the horizons, defined by isovalues in the unrestored volume, define layers in the uvt space in which the thickness variations are preserved. That permits to properly handle growth structures. The last step is a scaling of the time field to match the Z field, by preserving volume or vertical thickness during restoration. Such constraints approximate the minimization of the deformation during restoration [Medwedeff et al., 2016]. At the end, the deformed volume contains the restored geometry of the geological model at a specific deposition time, from which restored faults and horizons may be extracted.

Restoration of a structural sandbox model

Analog model: a structural extensional sandbox model

The comparison of the two restoration methods is based on the restoration of a laboratory sandbox model. The forward deformation of the sandbox was performed by IFPEN (http: //www.ifpenergiesnouvelles.fr) and C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System (http://www.ccreservoirs.com). This sandbox was initially composed by a basement of silicone (analogous to salt) and a layer of sand just above it (Figure 3.1). Then the sandbox was deformed by gravity (box inclination of 1.5 • ). During the deformation, alternating layers of pyrex or sand were deposited to mimic growth structures. The initial dimensions of the sandbox were 18 cm along the X axis, 10 cm along the Y axis and 2.2 cm along the Z axis. There is a scaling from the sandbox scale to real geological cases: 1 cm is equivalent to 1 km. This model is analogous to supra-salt structures which can be found in real geological basins such as in Gulf of Mexico or offshore Angola. For more details about the experiment and its applicability to real case studies, see Chauvin et al. [2017] and references therein.

Geological model

At the end of the experiment, a 3D computed tomography (CT) volume of the structural sandbox, equivalent to a seismic cube, is generated. From this volume, Chauvin et al. [2017] built a structural model with SKUA-GOCAD software [Paradigm, 2015]. This structural model (Figure 3.2) is the base model used for both restorations. However, each restoration method (GeoChron-based and mechanics-based) calls for a specific input. Indeed, for the mechanics-based restoration method a 3D tetrahedral mesh conformal to both horizons and faults is required. Chauvin et al. [2017] used the software VorteXLib [Botella, 2016a,b] which depends on the libraries RINGMesh [Pellerin et al., 2015], Geogram [Lévy, 2015] and TetGen [Si, 2015b,a]. VorteXLib software, before meshing the structural model with tetrahedra, remeshes the structural model surfaces. Thus, the structural model used in the mechanicsbased restoration is slightly different from the one produced by SKUA-GOCAD. In addition, Medwedeff et al. [2016] used a SKUA model [Paradigm, 2015], as input for the GeoChronbased restoration, generated from the initial structural model. Thus, the geological surfaces of the SKUA model (horizons and faults) may slightly vary from the initial structural model and the remeshed model used in the geomechanical restoration. Moreover, as horizons in the SKUA model are implicit (isovalues of a scalar field) and as the silicone is not represented in the numerical model, the horizon H1 was slightly and locally changed around faults to truncate them. This avoids faults terminating down-dip internal to the volume, and consequently fault block connection at the bottom of the model which would prevent correct sliding during restorations.

Figure 3.3 is a visual comparison between both input 3D models along three cross sections. The two input models are very similar. Figure 3.4 quantifies the mismatch between both input models and localizes the highest gaps. The principal dissimilarities are near the bottom horizon since it was reworked as said previously. The top of the model has also high dissimilarities. This is due to missing points on the top horizon, i.e., H8, to compute the distance with the explicit model. This is a numerical artifact since the top horizon was extracted from an isovalue which is numerically at the same location as the top of the mesh.

Restoration settings

Mechanics-based restoration settings

The boundary conditions used in the mechanics-based restoration of the sandbox model are in Chauvin et al. [2017] (see Figure 5 in Chauvin et al. [2017]). In short, the western wall is fixed in X and Y, the north and the south walls are fixed in Y, the uppermost horizon is flattened and there is an optional shortening along X on the eastern wall. Moreover, fault contact boundary conditions are defined to ensure the consistency between fault blocks in the restored state. In particular, Chauvin et al. [2017] used specific contact conditions to handle branching faults and faults which are cut and displaced by later faults. Mechanically the entire model is homogeneous and elastic with a Young's modulus of 70 GPa and a Poisson's ratio of 0.2 [Holtzman et al., 2009].

GeoChron-based restoration settings

We restored the sandbox model with the GeoChron-based restoration method [Medwedeff et al., 2016]. We chose the minimal deformation style. This tectonic style is more appropriate than flexural slip in this case study since the forward deformation is extensional and the analog model has no mechanical layering. For the time to depth scaling, we chose to preserve vertical thicknesses.

Restoration comparison

In this section we geometrically compare the restored states provided by both restoration methods. For simplicity, and to minimize the effect of error accumulation throughout restoration, we only compare the restored states for the first three restoration steps.

Geometrical comparison

To be able to compare both restoration methods, we scaled the GeoChron-based restored models to fit the W/E extension of the mechanics-based restored models. Concretely, we multiplied the X-component of each vertex of the GeoChron solid by a factor. This factor is equal to the extension of the geomechanical restored model, obtained with a proper shortening boundary condition (CT images), divided by the extension of the GeoChron-based restored model. respectively correspond to a visual comparison for each restoration step between both restored models on three cross sections. Visually, after scaling of the restored state, both restoration methods provide a very similar solution. To quantify the mismatches between both solutions, we computed the restoration displacement difference within the overall volume. Figures 3.8 -3.9 -3.10 respectively represent for each restoration step the norm of the restoration displacement difference within the elastic solid and its distribution. The majority of the norms is below 0.1 cm which is the picking uncertainty distance evaluated by Chauvin et al. [2017]. The extremal norms, above 0.2 cm, are mainly located around faults. Thus, both restoration methods provide very similar solutions in the continuous parts of the geological model.

The main mismatches are located near the faults. That is explained by a problem of transfer of the restoration vectors from the GeoChron solid to the geomechanical solid (in the unrestored state). Indeed, some values are assigned to the wrong fault block as illustrated in Figure 3.11. This effect is shown for a fault of the sandbox model in Figure 3.12b and Figure 3.12d. To deal with that, for analyzing the restoration displacement mismatches on faults, the transfer is performed using the fault surfaces with the information of the block side. As said previously, a fault in a model to restore is decomposed into two surfaces, one for each fault block. Let + and -denote both blocks. Restoration vectors are transferred from the fault surface on the GeoChron solid corresponding to the block + to the fault surface on the geomechanical solid corresponding to the block + (same principle for the side -). Thus, the previous issue of transfer disappears (Figure 3.12c and Figure 3.12e). Figures 3.13 -3.14 -3.15 respectively represent for each restoration step the norm of the restoration displacement difference on the faults and its distribution. These figures show that even near the faults the restoration displacement vectors are very similar for both restoration methods.

Fault compliance

Fault compliance is a validity criterion for a restored model. Both restoration methods define constraints to ensure neither gap nor penetration between fault mirrors. Figure 3.16 quantifies the fault compliance for each restoration step for both restoration methods. The histograms on this figure correspond to the distributions of the distance between the two mirror surfaces of a fault. Considering the limit of 0.1 cm defined by Chauvin et al. [2017], both restoration methods provide an acceptable compliance for faults. It can be noted that the mechanicsbased restoration, with a fault contact method based on contact mechanics [Wriggers and Laursen, 2006, Muron, 2005, Maerten and Maerten, 2006], provides a better fault compliance than the GeoChron-based method. However, hard conclusions cannot be made since fault compliance is in part dependent on the implementation choice. In addition, considering the geological model scale, several distances may be equivalent. One millimeter and one meter are mathematically very different, but both are negligible in comparison to a kilometer scale. Thus, even if the fault compliance is not perfect, the compliance is still admitted below a certain distance. Furthermore, fault compliance is impacted by the mesh size which is not the same between both models.

There is an issue related to faults cut and displaced by later faults. The continuity between fault parts is not ensured in the GeoChron-based restoration method, which provokes wrong fault configurations. In Figure 3.17a the orange fault is initially normal. After restoration (Figure 3.17c) it becomes reverse and it is not continuous. This is due to the sliding along the blue fault (Figure 3.17b) which is completely independent of the orange fault. This issue is observed during the restoration of the sandbox model by the GeoChron approach (Figure 3.18). Chauvin et al. [2017] avoid this issue by applying line contact boundary conditions. Thus, the mechanics-based restoration has the advantage to define custom boundary conditions to handle complex fault network. Similar conditions could be implemented in the GeoChron method in the future.

Some caution must be taken with contact conditions in general for both restoration methods. As they force the fault compliance, they may hide errors in the structural interpretations. In the case of a forced fault compliance, fault geometry and the restored strain should be analyzed to check for potential issues.

Extension recovery

In this section, we compare the restored models without integrating the knowledge of the sandbox forward extension, i.e., no scaling for the GeoChron-based restorations and no shortening boundary condition for the geomechanical restorations. We recomputed the geomechanical restorations without any shortening condition, i.e., the unrestored state is from the previous restoration step without any shortening boundary condition. There is no correction of the unrestored models throughout the sequential restoration as done by Chauvin et al. [2017] to a proper comparison with the GeoChron-based restoration method without any scaling. Thanks to the CT images it is possible to know the amount of extension of the structural sandbox model after restorations (along the X axis). For each restoration step, Table 3.1 displays the total forward extension of horizons H6, H7, and H8. This table presents the measurements on CT images, i.e., the reference extensions. The differences of west/east extension between the restored states and the last stage of the analog model deformation for both restoration methods are also indicated. It is clear that both restoration methods fail to recover the expected extension of the analog model. The results are better in the GeoChronbased method than in the geomechanical restoration. In the case of the mechanics-based restoration, a shortening boundary condition is necessary, for each restoration step, to properly recover the forward extension of the analog model [Chauvin et al., 2017]. Concerning the GeoChron-based restoration, the forward extension can be handled throughout restorations by an appropriate scaling [Medwedeff et al., 2016](Section 3.3.1). Both restoration methods need a specific shortening condition to be able to capture continuous deformation, extension of missing faults, and potential dilatation within the analog model [Chauvin et al., 2017]. CT images 1.44 (100%) 2.15 (100%) 4 (100%) GeoChron no W/E scaling 0.9 (63%) 1.715 (80%) 2.757 (69%) Mechanics no imposed shortening 0.718 (50%) 1.168 (54%) 2.238 (56%) Table 3.1: Total W/E forward analog model extensions for each restoration steps. Values from the CT images are the reference extensions. For each extension in centimeter, the percentage is relative to the corresponding reference extension. For the geomechanical restorations, for each restoration step there is no imposed shortening, i.e., the unrestored model is the restored model at the previous restoration step without shortening boundary condition.

Impacts of the mechanical properties in the geomechanical restoration

One of the major conceptual differences between GeoChron-based restoration and geomechanical restoration concerns how the deformation imposed on the restored horizon is propagated within the layers. In this section, we propose a sensitivity analysis of the restoration results by varying the elastic material properties. More precisely, we recompute the first three restoration steps with different sets of Young's modulus and Poisson's ratio, considered either homogeneous or heterogeneous within the sandbox model. We apply the same boundary conditions, including the shortening boundary condition evaluated from CT images. As the actual sand and pyrex properties in the analog model are the same [Panien et al., 2006], these tests are not expected to recover the reference solution. This can be interesting in real geological cases where mechanical layering may exist and should affect the deformation patterns.

Impact of Young's modulus

We geometrically compare different restored models with the reference restored model (Ref in Table 3.2). Only Young's modulus differs (models 1-6 in Table 3.2). Table 3.2 also presents the total restoration dilatation after each restoration step. The dilatation is always based on the initial unrestored volume (last deformation stage of the structural sandbox) using:

dilatation = 100 × V restored -V initial V initial . (3.5)
In (3.5) V initial only considers the layers which are present at the restoration step, i.e., the removed layers are not included in the volume calculation to have a consistent dilatation. Figure 3.19 presents two distance distributions between the reference restored model (the one from Chauvin et al. [2017]) and each of the other restored models. One distance is the norm of the difference of restoration displacements. The other distance is the volume dilatation relative to the reference restored model:

relative dilatation = 100 × V -V ref V ref . (3.6)
For each restoration test, the deviation to the reference is always significantly below the picking uncertainty distance limit of 0.1 cm defined by Chauvin et al. [2017]. Therefore, Young's modulus has almost no influence on the restored geometry and on the global dilatation. This is quite logical because Young's modulus relates the strain to the stress and that the strain is highly constrained by displacement boundary conditions. Thus in that case the geometry is poorly affected by Young's modulus, but the stress is. .32 -1.32 -1.32 -1.32 -1.32 -1.32 -1.32 Total dilatation step 2 (%) -1.85 -1.86 -1.85 -1.85 -1.85 -1.85 -1.85 Total dilatation step 3 (%) -3. 51 -3.51 -3.50 -3.47 -3.53 -3.50 -3.50 Table 3.2: Mechanical parameters of different geomechanical restorations. Only Young's modulus differs from the reference sequential restoration (E = 70 GPa, ν = 0.2) to test the influence of this parameter. Poisson's ratio is always equal to 0.2. The total restoration dilatation is displayed for each restoration step. The dilatation is relative to the initial unrestored model (last deformation stage of the structural sandbox). Young's modulus has no impact on the dilatation. Norm of the difference of restoration displacement between a restored model and the reference restored model (log10 scale, in cm)

Step 1

Step 2

Step 3

Volume dilatation relative to the reference restored model (log10 scale, in %) .84 -1.62 -1.32 -0.95 -0.47 -0.15 Total dilatation step 2 (%) -2.55 -2.27 -1.85 -1.34 -0.67 -0.24 Total dilatation step 3 (%) -4.81 -4.26 -3.53 -2.56 -1.33 -0.54 Table 3.3: Mechanical parameters of different geomechanical restorations. Only Poisson's ratio differs from the reference sequential restoration (E = 70 GPa, ν = 0.2) to test the influence of this parameter. Young's modulus is always equal to 70 GPa. The total restoration dilatation is displayed for each restoration step. The dilatation is relative to the initial unrestored model (last deformation stage of the structural sandbox). Poisson's ratio significantly impacts on the global dilatation.

ν = 0.2 E (Pa) ν = 0.2 ν = 0.2 10

Impact of Poisson's ratio

We compare different restorations for which only Poisson's ratio was varied (Table 3.3). Figure 3.20 presents the norm of the difference of restoration displacements of each restored model with the reference restored models. The volume dilatation relative to the reference restored models is also presented. For Poisson's ratio, the deviations are more important and are sometimes above the picking uncertainty limit. Poisson's ratio has a more significant impact on the displacement field than Young's modulus, but differences are relatively small for each node. However, the absolute dilatation (Table 3.3) decreases when Poisson's ratio increases. This is also logical since the lower Poisson's ratio is, the more the global volume decreases (for a non-null strain the volume does not change only if Poisson's ratio is equal to 0.5). Thus, although the impact of Poisson's ratio is locally quite small here on the displacement field, the consecutive change of the entire volume after restoration may be significant and must be checked. If the boundary conditions are considered as valid, Poisson's ratios which provide a non-geologically realistic dilatation must not be chosen.

Heterogeneous mechanical properties

We recomputed the first three restoration steps on the same structural sandbox but with heterogeneous mechanical properties. We chose the elastic properties of a sandstone (E = 10 GPa, ν = 0.2) and a shale (E = 70 GPa, ν = 0.4) as reported by Fossen [2016, p. 104]. We voluntarily chose extremal elastic behaviors to analyze the impact of a heterogeneous medium on the restoration. We tested two models (Table 3.4). The first one is an alternation of sandstone and shale layers. The second model has four upper layers of shale and three lower layers of sandstone. Figure 3.21 shows that for both tested models, after each restoration step, the majority of norms are below the picking uncertainty limit.

Step 3 of the model presenting thick layers is an exception with a large spike above 0.1 cm. The mechanical heterogeneity here implies a more significant difference with the reference homogeneous restored models.

Dilatation in the GeoChron-based restorations

Table 3.5 presents the total dilatation of the restored model after each restoration step of the GeoChron method. If no scaling is performed to fit CT images, the volume does not change after each restoration. Even if the conservation of volume is a common assumption in structural restoration, the GeoChron-based method cannot alone recover rock dilatation. However, the dilatation of the restored models can be handled with a proper scaling. The scaling we used permited to have a visual fit with CT images, but the resulting dilatation is very high. The dilatation values are above the ones provided by a geomechanical restoration with a Poisson's ratio of 0.01 (Table 3.3), which is logical since the scaling was only applied along the X axis without any compensation in the other directions (as null Poisson's ratio). Obviously, this is unrealistic. In the general case, a proper scaling in all directions must be performed to get a proper dilatation. We did not scale along the Z axis since the resulting displacement to preserve volume is at the picking uncertainty limit defined by Chauvin et al. [2017], and therefore this displacement would not change the geometric comparison between both restoration methods.

Discussions

Two restoration methods: "equivalent" restored states

The fact that both mechanics-based and GeoChron-based restoration methods furnish an "equivalent" restored state is reassuring. Neither of the methods is far from the reference Norm of the difference between both restoration displacement fields (log10 scale, in cm)

Step 3

Step 2

Step 3.4). The results are for the first three restoration steps. The distributions correspond to the norm of the difference of restoration displacements with the reference restored models.

Step 1

Step 2

Step 3 Not scaled Scaled Not scaled Scaled Not scaled Scaled Total dilatation (%) 0.00 -1.82 -0.00 -2.31 0.00 -5.30

Table 3.5: Total restoration dilatation for each restoration step of the GeoChron-based method. "Not scaled" refers to the restored model obtained after a restoration with the GeoChron method alone. "Scaled" refers to the restored model obtained after a restoration with the GeoChron method followed by a scaling to fit the CT images.

solution and both methods can be used. The similarity between both solutions may be simply explained. First, both restoration methods rely on minimization of the deformation. Second, the sandbox model is mechanically homogeneous, which makes the use of a global deformation style appropriate. Moreover, as the used mechanics is elasticity and that all the boundary conditions are imposed displacements, Young's modulus does not have a significant effect on the restored displacement. Only the stress is affected. If the boundary conditions were traction conditions, Young's modulus would have surely significantly affected the final geometry. These considerations may explain the similarities between both restoration solutions of the structural sandbox model. However, even if Young's modulus does not seem to have an important impact on the displacement field, for stress analysis purpose it still remains fundamental.

An incorrect Young's modulus would lead to an inconsistent stress field which has a major impact for instance on the prediction of fracture areas [Macé, 2006, Maerten and Maerten, 2006, Mejía-Herrera et al., 2014, Stockmeyer et al., 2017]. Nevertheless, the ability of the mechanics-based restoration to capture a physical stress field is controversial because of the elastic nature of the used mechanics and of the unphysical boundary conditions [Lovely et al., 2012]. Finally, as the strain mathematically only depends on the restoration displacement field, even if this field is obtained by a geometrical approach, such as the GeoChron method, it is possible to compute a consistent stress field if Young's modulus is properly chosen [Mallet, 2014, chap. 8].

Impacts of the mechanical properties

In the studied sandbox model, the strain provided by a geomechanical restoration on a homogeneous case is poorly affected by Young's modulus. Poisson's ratio has more impact, in particular on the volume dilatation. This is also a conclusion of Durand-Riard [2010].

In the case of knowledge about rock dilatation, the mechanics-based restoration enables to control, thanks to the boundary conditions and Poisson's ratio, the restored model volume, whereas the GeoChron approach preserves the original volume (assumption generally chosen in restoration), in the case of no postprocessing (here to horizontally scale the GeoChron restored model with the geomechanical restored model). A dilatation analysis study must be performed to assess consistent dilatation globally within a restored model and also locally.

Even in presence of mechanical heterogeneity, we did not get major changes in the restored models. This may be due to the displacement boundary conditions, which are numerous and highly constrain the strain. Moreover Poisson's ratio impact is function of the strain. In the case of more strain along the X axis, Poisson's ratio would have a more significant impact. It must be the case for the remaining restoration steps not performed in this study. Some caution must be taken on these results since they are based on the restoration of an extensional sandbox model. In actual cases, deformation and mechanics are far more complex. Several authors show the advantages of geomechanical modeling over geometric/kinematic methods in several contexts [e.g., Fletcher and Pollard, 1999, Maerten and Maerten, 2006, Guzofski et al., 2009].

Flexibility versus practicality

The GeoChron-based restoration has the main advantage to be easy to use and geologically accurate. In the current development made by Medwedeff et al. [2016] and Paradigm [2015], the number of settings which must be defined by a geologist is small: choice between two deformation styles, and choice between the preservation of the volume or the vertical thickness in the scaling of the t field to depth. This is an advantage over the mechanics-based restoration which needs a consistent set of boundary conditions (qualified by some authors as non-physical [Lovely et al., 2012]) and rock mechanical properties. However, this can be a disadvantage in a complex area which needs more flexibility in the restoration settings. For instance, the GeoChron-based restoration only proposes two extremes of deformation style (flexural slip and minimal deformation). It would be necessary to use an intermediate between these two end-members to cover all possible deformation styles [Mallet, 2014, p. 74]. In addition, in the GeoChron-based restoration the deformation style must be known whereas the mechanics-based restoration aims to be a restoration method in which the deformation style is a result of a mechanical simulation.

The GeoChron-based restoration also has the advantage to be fast on several aspects. First it uses a SKUA model [Paradigm, 2015] of which the generation and edition have numerous advantages over the traditional explicit methods in particular on the time of construction and the sensitivity on the structural data [e.g., Frank et al., 2007, Caumon et al., 2013]. This is a gain of time over the mechanics-based restoration which is generally based on an explicit structural model, with the exception of Durand-Riard et al. [2010]. Moreover, as the mechanics-based restoration is a sequential restoration and is sensitive to the mesh quality, a remeshing step may be necessary after a few sequential restoration steps [Chauvin et al., 2017]. Rebuilding an explicit structural model and remeshing it may be a very time-consuming task. The GeoChron-based restoration avoids this issue by restoring a horizon in a single step. Nevertheless, in this geometric restoration, in the workflow defined by Medwedeff et al. [2016], it is not real sequential restoration, but a "pseudo-sequential restoration", as the restored state is not used as unrestored state for a next restoration step. Mechanics-based restoration, in its classical use, has the advantage to perform a real sequential restoration. Finally, the GeoChron-based restoration, based on numerical optimization, was much faster in this case study (several minutes to restore several horizons) than geomechanical restoration (a couple of hours to restore each horizon). Even if software performance depends on many parameters and on the implementation, the capability to quickly validate or invalidate a structural model and go back to the initial interpretation is fundamental for testing multiple realizations and handling the structural uncertainties.

Conclusions

The restoration of a structural sandbox model in which the structural uncertainties are low and the paleo-geometries are known on a cross section enabled us to compare two restoration methods (mechanics-based and GeoChron-based). In this case study, both restoration methods provide a very similar restored state in terms of geometry and displacement field. There are several reasons for this similarity. First, although the sandbox model has complex structures, its deformation is quite simple. Second, both restoration methods are mainly based on similar geometric considerations (e.g., datuming) and on the minimization of the deformation. Finally, the sandbox is mechanically homogeneous and the mechanical parameters have a small impact on the displacement field. As there is neither pure traction boundary condition nor body force condition, and as the strain is much constrained by imposed displacement and fault contact conditions, Young's modulus has a very small impact on the displacement field. Poisson's ratio has a larger (yet not very significant) impact, and affects the total dilatation. The GeoChron-based method has significant advantages in terms of meshing ease, practicality and computational time. Nevertheless, the current implementation of this method has almost no flexibility, whereas custom boundary conditions can be set in the mechanics-based restoration method. Finally, these restoration comparison conclusions cannot be simply transposed on models presenting high mechanical contrasts and complex deformations, which can be the case on real geological case studies.

• Geomechanical second order tensors field visualization through the segmentation and animation of geological 3D models, Côme Lebreton (Master student, 2015).

• Integration of paleomagnetic data in 3D mechanics-based restoration, Anne Raingeard (ENSG 2 nd year student, 2015).

• Accounting for pore pressure during decompaction, Laure Pizzella (Master student, 2016).

• Implementation of RINGMecha Graphical User Interface, Anne Raingeard (ENSG 3 rd year student, 2016).

• Paleo-topography boundary condition in 3D mechanics-based restoration, Camille Philippe (ENSG 2 nd year student, 2016).

• Validation of a finite element program, Charles Vouaux (ENSG 2 nd year student, 2016).

Developed software

During this thesis, we developed the following computer programs:

• Restorationlab: 3D mechanics-based restoration SKUA-GOCAD plugin1 [Muron, 2005, Durand-Riard, 2010].

• RINGMecha: finite element library for mechanical simulations2 . This work was in collaboration with Antoine Mazuyer.

• RINGMesh: open-source mesh structure library3 . This work was in collaboration with Jeanne Pellerin, Arnaud Botella, Antoine Mazuyer, François Bonneau, Pierre Anquez, Margaux Raguenel and Gautier Laurent. See Pellerin et al. [2017] for more details.

These programs, and the related trainings, are available for the members of the RING-Gocad consortium in http://www.ring-team.org/software. RINGMesh is available on the RING team website or on bibucket: https://bitbucket.org/ring_team/ringmesh.

Perspectives of this thesis

In this thesis we tried to provide some answers to the current issues of the 3D mechanics-based restoration. In this section, we propose several perspectives in mechanics-based restoration. They are related to the meshing constraints, the choice of the boundary conditions and the physical applicability of the restoration.

Meshing constraints

Problematic

Geomechanical restoration methods require a mesh of the geological model to restore. The mesh is built prior the first step of restoration. It may be necessary to rebuild a mesh between two restoration steps of a sequential restoration to geometrically correct the restored mesh or to improve the numerical quality of the deformed mesh elements [e.g., Parthasarathy et al., 1994, Shewchuk, 2002, Botella, 2016]. However, building or rebuilding a mesh is not straightforward. In the restorations presented in Chapter 2, we could not perform an entire sequential restoration because of the issues on the mesh element quality and of the difficulty to rebuild a 3D numerical model. Moreover, meshing constraints may also force some simplifications on the numerical model, which may have an impact on the geological consistency of the numerical model [Vidal-Royo et al., 2012].

Constraints on the structural model

The different methods of geomechanical restoration need a 3D numerical structural model. This model can be directly used for restoration purpose or as input to produce a 3D mesh. There are several constraints on the construction of a valid geological boundary representation [Caumon et al., 2004, Zehner et al., 2015, Anquez et al., 2016]. Among them, two have an important impact on the usability of the geomechanical restoration. The first constraint is to ensure the conformity between the surfaces of the structural model. Two surfaces are conformal if their intersection is a line which is present in both surface meshes [Caumon et al., 2009]. The second constraint is that the boundary representation must be water-tight, i.e., there is no communication between the interior of the numerical model and the exterior. These constraints make the construction of a structural model a complicated task. In addition, for restoration purposes, along each fault surface the mesh is duplicated to enable the sliding. Thus there are two surfaces for each fault and the conformity is lost after restoration [Anquez et al., 2016]. Moreover, fault compliance, although good, cannot be numerically perfect. Therefore, the generation of a boundary representation from a restored fault network is not trivial.

Constraints on the volumetric mesh

Concerning the volumetric mesh, its generation, respecting the topology of a structural model, is the topic of several works in the literature (see Botella [2016] and references therein). A difficult task is to respect all the geological interfaces, i.e., faults and horizons. Mesh algorithms may fail to reproduce some geometrical configurations, such as low angles or very thin layers. Geometrical/topological simplifications are often required before meshing, leading to a simplification of the geological interpretation [Pellerin et al., 2014, Raguenel et al., 2016]. Smaller volumetric elements, such as tetrahedra, help to respect the structural model. But such a method implies a high number of elements to honor thin features, which significantly increases computer memory and computational time requirements. Another constraint, specific to the restoration, is the potential need of mesh nodes on particular locations to define local boundary conditions such as the piercing points proposed by Durand-Riard et al. [2013].

Some ideas to address meshing issues

Several methods are used or could be used to avoid, or at least reduce, meshing constraints. Implicit modeling may ease the generation of a structural model [e.g., Chilès et al., 2004, Cowan et al., 2004, Frank et al., 2007, Calcagno et al., 2008, Paradigm, 2015]. In such method, horizons are defined by isovalues of a scalar property [Frank et al., 2007]. However, the conversion of an implicit model to an explicit model is not simple. This problem is more related to structural geomodeling than to the restoration but it highlights the issues of current restoration requirements. Implicit restoration reduces the constraints on the generation of a structural model since the horizons are not represented by surfaces [Durand-Riard, 2010, Durand-Riard et al., 2010]. However conformal fault surfaces are still needed. The use of a fully implicit model, i.e., with implicit horizons and faults, as input for restoration may solve the meshing issues. The extended finite element method (X-FEM) is an extension of the finite element method to implicit models. Extended finite element method [e.g., Moës et al., 1999, Moës and Belytschko, 2002, Belytschko et al., 2013] was used by Siavelis et al. [2010Siavelis et al. [ , 2011]], Siavelis [2011], Annavarapu et al. [2013] and Siavelis et al. [2013] to model the forward evolution of sedimentary basins and complex fault contacts with large fault sliding. Even if Siavelis [2011] did not perform a mechanics-based restoration using X-FEM, the same principle as forward simulation can be applied, easing mesh generation and contacts between faults. A more advanced method to avoid meshing problems is the use of meshless methods [Belytschko et al., 1996, Müller et al., 2004, Maerten, 2014, Renaudeau et al., 2016]. During this thesis, we tried to develop implicit boundary conditions on points located anywhere within a 3D mesh [Chauvin andCaumon, 2014, 2015b] by generalizing the work of Durand-Riard [2010] and Durand-Riard et al. [2010]. These conditions are transferred to the surrounding nodes using tetrahedral barycentric coordinates. This approach also aimed to bring flexibility of the boundary condition definition, opening the path to study subsequent uncertainties [Chauvin and Caumon, 2014]. We did not develop this work in this manuscript because it did not succeed, surely because the definition of "boundary conditions" within a mesh is not considered in the classical finite element method. In addition, our approach to flatten a horizon applies boundary conditions on all the nodes of the tetrahedra containing the implicit points, leading surely to an overconstrained system. An alternative method would consider only the nearest nodes as done by Durand-Riard et al. [2010] to flatten a horizon. Finally, external works on meshing enable now the use of multi-element mesh, such as hex-dominant mesh [Botella et al., 2014b[Botella et al., , 2015] ] and local adaptive mesh [Botella et al., 2014a]. Concerning the contruction/reconstruction of a structural model, a thesis is currently trying to generate a valid geological boundary representation from a set of non-conformal surfaces [Anquez et al., 2016]. That reduces the difficulties on meshing, avoids some simplifications of geological contacts and may harness the number of mesh elements.

Physical considerations

Choice of the boundary conditions

Several studies, including this thesis, expose the problem of defining consistent boundary conditions [Durand-Riard, 2010, Lovely et al., 2012, Durand-Riard et al., 2013]. We proposed a set of boundary conditions, validated on an extensional analog system with known deformations. A future work would be to further test them in other deformation contexts [Durand-Riard, 2010, Durand-Riard et al., 2013] and in natural case studies. In addition, in our work we used the area-depth method to estimate the shortening boundary condition. We compute the lost area divided by the depth to the detachment. Unfortunately, as the analog model has a growth stratigraphy, we could not use an area-depth graph which uses the lost area and the depth to detachment of several horizons, leading to more precise result [Groshong et al., 2003]. The use of such a graph in case of pre-growth strata or in case of no-growth intervals [Groshong, 2015], can be used to estimate the magnitude of the lateral displacement boundary condition. Moreover, the area-depth method is perfectly adapted to compressive contexts [Groshong and Epard, 1994, Groshong, 2006, Groshong et al., 2012, Groshong, 2015]. Concerning the novel fault contact conditions, we directly applied them on a complex fault network. To properly evaluate their importance, tests on small synthetic models should be performed. For instance, analyzing the geometrical evolution of two or three connected faults could provide insight on the validity of these novel contact conditions. Another work would be to define global fault contact condition. Currently, contact conditions are set to tie the two mirrors of a fault [e.g., Muron, 2005]. Thus the contacts are defined for each fault, and a fault mirror is constrained to stay on the other fault mirror, no sliding is allowed over several faults. The fault network should be seen as a unique sliding surface and not as an ensemble of discrete sliding surfaces. Such a global contact condition should improve the contact between connected faults and handle zip junctions [Platt andPasschier, 2016, Passchier andPlatt, 2016]. Numerically, instead of applying a fault contact fault by fault, it should be applied on all the faults at the same time considering all the faults as a unique surface.

Uncertainties of the mechanical parameters

Several challenges on restoration come from the used physics. For instance, elastic models are used whereas tectonic deformations are plastic. This assumption is from the need of reversibility of the method and for simplicity reasons [e.g., Muron, 2005, Moretti et al., 2006, Maerten and Maerten, 2006, Guzofski et al., 2009]. Restoration could be coupled with forward mechanical simulation integrating plasticity [Siavelis, 2011]. Indeed, from a restored model using elasticity, a coherent forward simulation could be realized to validate the restored model. However, such an approach would require several iterations before getting a forward model matching with the present day geometry, which is time-consuming. Moreover, that implies a good knowledge of the mechanisms which generated the deformations through time.

Concerning the elastic parameters, i.e., Young's modulus and Poisson's ratio, more sensitivity studies should be performed to assess their importance on the restored geometries. For instance, our tests on an extensional analog model could be extended to other deformation contexts (compressive, strike-slip, etc.) and on models presenting high mechanical contrasts. Moreover, new comparisons with the GeoChron-based restoration method could provide more insight on the limits of both methods. Finally, the choice of the elastic parameters is impacted by the rock uncertainties within a geological domain. Indeed, petrophysical data are only observed along wells and outcrops. Uncertainties arise from the interpretation of these data [e.g., Prasad et al., 2002, Chang et al., 2006, Ameen et al., 2009] and the petrophysical interpolations/simulations within the entire geological domain [e.g., Bertoncello et al., 2008, Lallier et al., 2009, Abdideh and Ghasemi, 2014, Bennis et al., 2014]. Finally, there are also uncertainties due to a necessary upscaling/homogenization of the mechanical parameters in order to reduce the number of volumetric mesh elements [e.g., Chalon et al., 2004, Bayuk et al., 2008, Durand-Riard, 2010, Durand-Riard et al., 2010, Capdeville et al., 2015, Cupillard and Botella, 2015].

"Odi panem quid meliora. Ça veut rien dire, mais je trouve que ça boucle bien."

[Le roi Loth, Kaamelott, Livre V, tome 1].

Applicabilité de la restauration géomécanique : conditions aux limites, réseau de failles et comparaison avec une méthode géométrique Résumé : La restauration structurale a pour objectifs de déterminer la géométrie passée des roches et de valider les interprétations structurales. Les méthodes classiques sont basées sur des hypothèses géométriques et/ou cinématiques, et imposent un style de déformation. Les méthodes géomécaniques, en intégrant le comportement élastique des roches et les lois fondamentales de conservation mécanique, visent à résoudre les problèmes des méthodes classiques. Toutefois, il y a des incertitudes sur le choix des paramètres élastiques, et les contraintes de maillage rendent difficile l'utilisation de cette méthode comme un outil de validation des interprétations structurales. Le choix d'une méthode de restauration en particulier est rendu difficile par le fait qu'il y ait plusieurs approches de restauration géomécanique, en plus des nombreuses méthodes géométriques et cinématiques. Cette thèse présente en premier lieu une revue des différentes méthodes géomécaniques 3D visant à déplisser et annuler l'action des failles dans un modèle géologique. L'objectif de cette revue est de présenter les forces ainsi que les limites, théoriques et pratiques, de chaque méthode. Dans un second temps, à travers la restauration d'un modèle analogique (sandbox ), nous présentons nos travaux sur le choix de conditions aux limites appropriées pour obtenir un modèle restauré cohérent. Ce modèle structural expérimental a été déformé en laboratoire et présente plusieurs analogies avec des structures extensives postérieures à une base salifère. Grâce à l'observation de l'évolution temporelle de la géométrie du modèle analogique sur une coupe, nous montrons qu'une condition aux limites correspondant à un raccourcissement latéral est nécessaire. Ce raccourcissement peut être estimé par la méthode de la surface transférée. De plus, nous définissons de nouvelles conditions aux limites de contacts de failles pour restaurer correctement le réseau de failles complexe du modèle analogique. Ces nouvelles conditions lient les bords internes des surfaces de failles et connectent les composantes connexes des failles coupées et déplacées par des failles plus récentes. Troisièmement, le test de différents paramètres élastiques indique que le module de Young, défini homogène au sein d'un modèle géologique, n'a quasiment pas d'effet sur le champ de déplacement. Toutefois, le coefficient de Poisson a un impact significatif sur la dilatation volumique. Dans un dernier temps, nous comparons la restauration géomécanique avec une méthode géométrique qui repose sur un modèle chronostratigraphique (GeoChron) qui fait une bijection de chaque point du sous-sol avec son équivalent dans l'espace de dépôt (Wheeler). Nous montrons que les deux approches de restauration fournissent des modèles restaurés du modèle analogique qui sont similaires géométriquement. La méthode géométrique a de nombreux avantages pour obtenir rapidement et avec précision le modèle restauré, mais elle manque de flexibilité sur le choix des contraintes de la déformation. La force de la méthode géomécanique est de pouvoir définir des conditions aux limites personnalisées et des comportements mécaniques spécifiques pour gérer les contextes mécaniquement complexes.

Mots-clés : restauration géomécanique, conditions aux limites, modèle analogique, GeoChron, réseau de failles Applicability of mechanics-based restoration: boundary conditions, fault network and comparison with a geometrical method Summary: Structural restoration aims to recover rock paleo-geometries and to validate structural interpretations. The classical methods are based on geometric/kinematic assumptions and impose a style of deformation. Geomechanical methods, by integrating rock elastic behavior and fundamental mechanical conservation laws, aim to solve issues of classical methods. However several studies show that the geomechanical restoration lacks physical consistency in particular because of the boundary conditions. There are uncertainties on the choice of the elastic properties, and the meshing constraints limit this method to be used as a validation tool of structural interpretations. The choice of a specific restoration method is difficult because there are many geomechanical restoration approaches, in addition to the numerous geometric/kinematic methods. Firstly, this thesis presents a review of the various 3D geomechanical methods to unfold and unfault a 3D geological model. The objective is to present their, theoretical and practical, strengths and limits. Secondly, through the restoration of a structural sandbox model, we worked on the choice of adequate boundary conditions to get a proper restored model. This structural sandbox model was deformed in laboratory and presents several analogies with supra-salt extensional structures. Thanks to the observation of the analog model geometry through time on a cross section, we show that a lateral shortening boundary condition is necessary. We show that this shortening can be estimated by the area-depth method. Moreover we define new fault contact conditions to handle complex fault networks. These novel conditions tie internal fault borders and join parts of offset faults. Thirdly, the test of several elastic parameters shows that Young's modulus, homogeneous within a geological model, has almost no effect on the restoration displacement field. However, Poisson's ratio has a significant impact on the volume dilatation. Finally, we compare the mechanics-based restoration method with a geometric-based method relying on a chronostratigraphic model (GeoChron) mapping any point of the subsurface to its image in depositional (Wheeler) space. We show that both methods provide a geometrically similar restored state for the analog model. The geometric method has numerous advantages to quickly and accurately get a restored model, but it lacks flexibility on the choice of the deformation constraints. The geomechanical restoration method force is to define custom boundary conditions and specific mechanical behaviors to handle complex contexts.

Keywords: mechanics-based restoration, boundary conditions, analog model, GeoChron, fault network

  model (dark blue) and the mechanics-based input model (orange) . . . . . . . 3.4 Distribution of the distance between the two initial models of each restoration method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Visual comparison of the restored state between the GeoChron-based restoration (dark blue) and the mechanics-based restoration (orange) for the first restoration step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v LIST OF FIGURES 3.6 Visual comparison of the restored state between the GeoChron-based restoration (dark blue) and the mechanics-based restoration (orange) for the second restoration step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Visual comparison of the restored state between the GeoChron-based restoration (dark blue) and the mechanics-based restoration (orange) for the third restoration step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1 : Modèle structural du modèle analogique. Le modèle est composé de 8 horizons et 22 failles normales. Il y a deux grabens (à l'ouest et au centre) et un demi-graben (à l'est).

2. 4 . 1

 41 La méthode de la surface transférée (area-depth method )

3

  http://www.ring-team.org/software/ring-libraries/44-ringmecha xxii Section 3. Comparaison entre la restauration géomécanique et la restauration basée GeoChron :

1. 1 . 1 . 1

 111 Equations of motion: mass and linear momentum conservations

Figure 1 . 1 :

 11 Figure 1.1: Restoration traditional boundary conditions. Displacement and/or contact conditions are set to unfold and/or unfault a geological model. Data courtesy of Total.

Figure 1 . 2 :

 12 Figure 1.2: Equivalent mechanical properties. Faulted model with an implicit horizon (black thick line). This horizon is the boundary of two layers (red and blue). Average Young's modulus and Poisson's ratio is defined for each tetrahedron crossed by the implicit horizon.

Figure 1 .

 1 Figure 1.3: Sequential restoration workflow. a) Boundary representation of a geological model (structural model). b) 3D mesh, here filled by tetrahedra. c) Mechanical parameters (Young's modulus E and Poisson's ratio ν) assigned in the mesh. These parameters may vary within the volume. d) Definition of the restoration boundary conditions to unfault and unfold the model. e) Computation of the restoration displacement field by a numerical method. f) Restored state. g) Removal of the uppermost layer. The sequential restoration continues at step d). Data courtesy of Total.

Figure 1 . 4 :

 14 Figure 1.4: Mesh discretization. a) Theoretical shape of the medium Ω 0 . b) Approximation of the shape of the medium Ω 0 . This new domain Ω0 depends on the mesh element choice Ωe 0 . The equation of motion is solved on each mesh element.

Section 1 . 2 .Figure 1

 121 Figure 1.5: Mass-spring approach. a) Mesh of a 2D surface, composed by triangles. b) Example of spring and mass point disposition.Here the mass points are on mesh nodes, and the springs correspond to the mesh edges. c) Displacement of a node of the mesh (green node) due to boundary conditions. The springs connected to this node (green springs) are in elongation. Blue nodes and blue springs do not move yet. d) The elongation of the springs in c) generates forces, controlled by the spring stiffnesses, acting on the nodes connected to them (green nodes). The motion of these nodes implies a length change of the springs connected to them (green springs). Thus the motion propagates, spring by spring, with a damping which progressively dissipates the energy, until the system reaches an equilibrium.

Chapter 1 .

 1 Mechanics-based restoration review 0 Newton-Raphson: Find such as : function to solve such as : displacement at step

Figure 1 . 7 :

 17 Figure 1.7: Newton-Raphson algorithm principle. From an initial solution u 0 the real solution is progressively computed. Each step of the algorithm computes the tangent of the function f at the abscissa u i (previous solution) to define a new approximation of f such as f (u i+1 ) = 0. The process stops when u i+1 -u i is small enough.

FaultFigure 1 . 8 :

 18 Figure 1.8: Node to edge condition. A s is a node on the slave fault cutoff line. It is projected onto the master fault cutoff line. The image of A s onto the master edge A m 0 A m 1 is A m (defined by its barycentric coordinates in the edge).

Figure 1 . 9 :

 19 Figure 1.9: Node to triangle condition. A s is a node on the slave fault mirror. It is projected onto the master fault mirror. The image of A s onto the master triangle A m 0 A m 1 A m 2 is A m which is defined by its barycentric coordinates. Data courtesy of Total.

Figure 1 .

 1 Figure1.10: Non-linear geometrical contact. First step of an iterative contact process. The contact is between the red line (slave) and the blue line (master). Red dots are the slave line nodes. a) Slave and master lines to connect (as an unrestored state). b) Definition of the projection vector P for each slave nodes (no specific rule). c) Fictitious step in which the nodes are moved to cancel the projection vector gap without motion of the master. d) Master line moves due to the mechanical parameters and the other boundary conditions. Slave nodes are still not on the master line. This highlights that the contact process is not linear. From d) new projection vectors are defined and a new mechanical simulation is performed until the slave nodes are on the master line. Step d) is as an iteration in the Newton-Raphson loop.

Figure 1 .

 1 Figure 1.11: Stress vector decomposition. The stress vector of the master t m is decomposed into a normal component t m n and a tangential component t m t .In frictionless contact, t m t is null. nm and τ m are respectively the normal and the tangential vectors at the point xm on the master boundary. nm is oriented toward the exterior of the master domain Ω m . These relations are the same for the slave.

Figure 1 .

 1 Figure1.13: Non-linear mechanical contact. First step of an iterative contact process. The contact is between the red line (slave) and the blue line (master). Red dots are the slave line nodes. For a better understanding, the master is fixed in the figure even if it is deformed and moved in real mechanical contact processes. a) Slave and master lines to connect (as an unrestored state). b) Definition of the normal gap g n for each slave nodes. c) Fictitious step in which the nodes are moved to cancel the normal gap without tangential motion. d) Frictionless contact does not define the tangential gap. Slave nodes tangentially move as a result of the mechanical simulation (consequence of the chosen mechanics and the other boundary conditions). Slave nodes are still not on the master line. This highlights that the contact process is not linear. From d) new normal gaps are defined and a new mechanical simulation is performed until the slave nodes are on the master line. Step d) is as an iteration in the Newton-Raphson loop.

  Assuming small deformations, there is a linear relation between the Cauchy stress σ and the linear Green-Lagrange strain . Let D denote the fourth order elasticity tensor. Hooke's law is defined by [e.g.,Ramsay and Huber, 2000]: σ = D : .(1.42)

Figure 1 .

 1 Figure 1.14: Transverse isotropic axes. x and y axes define isotropic (stratigraphic) planes. z is the axis of anisotropy (∼stratigraphy normal).

Figure 2 . 1 :Figure 2 . 2 :

 2122 Figure 2.1: Scheme of the analog model experiment. An initial pregrowth stratigraphy composed by a layer of silicone and a layer of sand is deformed by gravity. The initial dimensions of the structural sandbox are: 18 cm (7 in) along X axis, 10 cm (4 in) along Y axis, and 2.2 cm (0.9 in) along Z axis. At each time step, a layer of sand or pyrex is deposed to generate syn-sedimentary deformations. Twelve layers are deposed forming a growth stratigraphy (3.1 cm, 1.2 in). Figure created from C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System and IFPEN documentation.

Figure 2 Figure 2 . 4 :

 224 Figure 2.3: CT image interpretation and model simplifications. Two examples of interpretation and simplification on the northern edge of the 3D tomography (interpreted as a seismic cube). (A) North CT image at the final stage of the structural sandbox experiment. Each example is in black square. (B)-(C) Interpretation examples composed from the left to the right by: the uninterpreted CT image part, its high order interpretation, and its simplification. The simplified interpretation is the one used for the restorations. Dashed lines are the interpreted elements not kept in the final model. Continuous lines are the interpreted elements kept in the final model. Faults are in black. Horizons are in white or gray. (B) Interpretation and simplifications on a part of the western graben. (C) Several faults delimiting relatively small fault blocks were neglected to facilitate 3D mesh generation. Horizons were made continuous, introducing local mismatches with the data. CT data courtesy of IFPEN and C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System.

Figure 2 . 5 :

 25 Figure 2.5: Restoration boundary conditions. (A) Volumetric model at the final stage of the structural sandbox experiment. The presented boundary conditions are applied on this model. Equivalent conditions are assigned to restore the other horizons (sequential restoration). (B) Western wall: no motion is allowed along X and Y directions. (C) Northern and southern walls: no motion is allowed along Y direction. (D) Uppermost horizon: flattening. This condition explicitly constrains Z component of mesh nodes. (E) Uppermost horizon parts are tied. It is a fault contact condition. For each fault cutting the uppermost horizon, the footwall cutoff line and the hanging wall cutoff line are tied together. This condition is applied on X and Y directions. Z component is ensured by the flattening, see (D). (F) For each fault, fault mirrors are tied (neither gap nor penetration). (G) Eastern wall: shortening along X direction. This condition is optional.

Figure 2 . 6 :

 26 Figure 2.6: Fault definition in two sides. (A) Fault and horizon topology as in a boundary representation of a geological model. Fault F (single surface) splits a horizon H in two parts. (B) Fault and horizon topology as in a volumetric model to restore. Fault F corresponds to two surfaces, one for the hanging wall (F -) and one for the footwall (F + ). Arrows represent fault contact boundary conditions used in restoration. These conditions tie both surfaces which compose F, allowing sliding without friction.

Figure 2 . 7 :

 27 Figure 2.7: Contact condition to connect fault internal borders. The subfigures have different views as indicated by the axes. (A) Restored model (first step) with an eastern shortening of 1.44 cm.The studied contact is between F1 and F2: F2 branches onto F1. (B) Contact between F1 and F2 after restoration when no contact condition (B1) is set to ensure a proper contact between these faults: there is a hole between the connected components that compose F1 (B2 and B3). (C) Contact between F1 and F2 after restoration when a contact condition (C1) is set to ensure a proper contact between these faults (black arrows). Continuity (within the contact precision scale) exists between the connected components that compose F1 (C2 and C3). In (B1) and (C1) gray arrows represent classical fault contact conditions, and in (B2), (B3), (C2) and (C3) X-coordinate contours are displayed on fault surfaces. + and -signs are an arbitrary convention to make the distinction between fault sides (seeFigure 2.6)

Figure 2 . 8 :Figure 2 . 9 :Figure 2 .Figure 2 .Figure 2 .

 2829222 Figure 2.8: Contact condition to fit footwall and hanging wall fault surfaces. (A) Initial unrestored model. This example focuses on a fault F1 cut into two parts by another fault: a hanging wall part F1-hw and a footwall part F1-fw. The hanging wall and footwall definition of F1 is relative to the fault that cuts F1: F2. In all the remaining subfigures, F2 is visible by its border, and X-coordinate contours are displayed on F1 surface. (B) Initial unrestored shape of F1. It is composed of two disconnected connected components. (C) Shape of F1 after the first restoration step. F1 is still into two parts. (D) Shape of F1 after the second restoration step. Both connected components of F1 are partially connected on the southern side. (E) Shape of F1 after the third restoration step. It seems visually that F1 should be continuous. (F) Shape of F1 after the third restoration step as (E) with additional contact constraints to ensure continuity between the hanging wall and the footwall.

Figure 2 .

 2 Figure2.13: Dip slip delta distributions. Difference of dip slip measures between the restored model and the reference for the first four restoration steps (with shortenings). For all the restoration steps except the fourth one, the majority of the dip slip deltas is within the uncertainty range (-0.2 cm and 0.2 cm).

Figure 2 .

 2 Figure2.14: Area-depth method principle. Area-depth method applied on the north CT image at the final stage of the structural sandbox experiment. Examples of the calculus for two horizons: H8 and H3. For each horizon a lost area (white area with dashed border line) is computed. The lost area for a horizon is the area above this horizon and below the regional level of this horizon. The height to the detachment is computed for each horizon. Total shortening (from the beginning of the deformation) undergone by a horizon is its lost area divided by its height to the detachment. CT data courtesy of IFPEN and C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System.

Figure 2

 2 Figure2.15: 3D dilatation analysis. Proportion of tetrahedra with a positive dilatation after each restoration step (steps 1 to 4) according to different shortenings (in centimeters). Black dots are data points (restoration simulations). For each restoration, the unrestored state is the restored state at the previous restoration step with the imposed shortening from CT image (and without the uppermost restored layer). Shortenings in Table2.3 are displayed on each graph. Ref: shortening from the CT image. AD: shortening from the area-depth method. HS: sum of the fault heaves of uppermost horizon fault cutoff lines on the northern wall. BL: shortening from the bed length conservation method. NIS: no imposed shortening on the northern wall (restoration without shortening condition).

1

  Article to be submitted to AAPG Bulletin.Authors: Benjamin P. Chauvin 1 , Peter J. Lovely 2 , Stanislas N. Jayr 3 , Guillaume Caumon 1 GeoRessources, Université de Lorraine / CNRS / CREGU, ENSG, Vandoeuvre-lès-Nancy, France 2 Chevron ETC-Integrated Exploration Research Team, Houston, TX 77002, USA 3 Chevron ETC-Earth Science Computational Technology, Houston, TX 77002, USA Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 3.1 Restoration methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.1.1 Mechanics-based restoration . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.1.2 GeoChron-based restoration . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.1.2.1 The GeoChron model . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.1.2.2 A restoration method based on the GeoChron model . . . . . . . 105 3.2 Restoration of a structural sandbox model . . . . . . . . . . . . . . . . . 106 3.2.1 Analog model: a structural extensional sandbox model . . . . . . . . . . . . 106 3.2.2 Geological model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 3.2.3 Restoration settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 3.2.3.1 Mechanics-based restoration settings . . . . . . . . . . . . . . . . . 109 3.2.3.2 GeoChron-based restoration settings . . . . . . . . . . . . . . . . . 109 3.3 Restoration comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 3.3.1 Geometrical comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 3.3.2 Fault compliance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 3.3.3 Extension recovery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 3.4 Impacts of the mechanical properties in the geomechanical restoration 123 3.4.1 Impact of Young's modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 3.4.2 Impact of Poisson's ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 3.4.3 Heterogeneous mechanical properties . . . . . . . . . . . . . . . . . . . . . . 128

Section 3. 2 .

 2 Restoration of a structural sandbox model

Figure 3 . 1 :Figure 3 . 2 :

 3132 Figure 3.1: Scheme of the structural sandbox experiment. Based on the documentation furnished by IFPEN and C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System.

Figure 3 . 4 :

 34 Figure 3.3: Visual comparison of the unrestored state between the GeoChron-based input model (dark blue) and the mechanics-based input model (orange). In each sub-figure both models are displayed. a) Overlaid 3D models. The CT image and the two rectangles with black borders are the cross sections presented in b), c) and d) respectively on the north, the center and the south of the model. Both models are visually very similar along the model. CT data courtesy of IFPEN and C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System.

Figure 3 . 6 :Figure 3 . 7 :Figure 3

 36373 Figure 3.5: Visual comparison of the restored state between the GeoChron-based restoration (dark blue) and the mechanics-based restoration (orange) for the first restoration step. In each sub-figure both models are displayed. a) Overlaid 3D restored models. The CT image and the two rectangles with black borders are the cross sections presented in b), c) and d) respectively on the north, the center and the south of the model. Both restored models are visually very similar along the model. CT data courtesy of IFPEN and C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System.

Chapter 3 .Figure 3 . 9 :

 339 Figure 3.9: Volumetric mismatch between both restorations at restoration step 2. a) Distribution of the norm of the difference of restoration displacement between the restored models obtained by the mechanics-based restoration and the GeoChron restoration (step 2). b)-f) Location of the highest distances within the unrestored model (near faults).

Figure 3 Figure 3

 33 Figure 3.10: Volumetric mismatch between both restorations at restoration step 3. a) Distribution of the norm of the difference of restoration displacement between the restored models obtained by the mechanics-based restoration and the GeoChron restoration (step 3). b)-f) Location of the highest distances within the unrestored model (near faults).

  Figure 3.12: Example of bad property transfers in the analog model. a) Analog model at the last deformation stage. b) and d) Both sides, respectively F -and F + , of a fault F in the geomechanical structural model. The displayed property is the Z component of the restoration displacement field obtained with the GeoChron-based method. The transfer of the restoration displacement, from the GeoChron volumetric mesh to the geomechanical volumetric mesh, failed: there is no continuity in the restoration displacement field. Restoration vectors are in white. c) and e) Same fault with a correct transfer: the restoration displacement field is continuous on each fault side.

Figure 3 Figure 3

 33 Figure 3.13: Norm of the difference between restoration displacements on faults (step 1). a) Distribution of the norm of the difference of restoration displacement between the restored faults (two surfaces for each fault) obtained by the mechanics-based restoration and the GeoChron restoration (step 1). b)-f) Location of the highest distances within the unrestored faults.

Figure 3 <Figure 3

 33 Figure 3.15: Norm of the difference between restoration displacements on faults (step 3). a) Distribution of the norm of the difference of restoration displacement between the restored faults (two surfaces for each fault) obtained by the mechanics-based restoration and the GeoChron restoration (step 3). b)-f) Location of the highest distances within the unrestored faults.

Figure 3

 3 Figure 3.18: Problem of no relation between fault parts (offset fault) in the sandbox model. a) Analog model at the last deformation stage. b) Example of a fault cut and displaced by a later fault in the unrestored state. c) and d) Restored state (step 3) after respectively a GeoChron-based restoration or a geomechanical restoration. Gaps and penetrations are present in the GeoChronbased restoration, and are absent in the geomechanical restoration. e) and f) are respectively the same as c)and d), but F 4 -and the part of F 3 located on the footwall side of F 4 were removed to ease the observation of the inconsistent reverse shape.

Figure 3 .

 3 Figure 3.19: Impact of Young's modulus on the restoration displacement field (geomechanical restoration). Each curve corresponds to a restoration test with a different Young's modulus. The norm of the difference of the restoration displacement for each restoration test with the reference model (E = 70 GPa, ν = 0.2) are in the first row. The dilatation in percent for each restoration test relatively to the reference model (E = 70 GPa, ν = 0.2) are in the second row. Young's modulus has almost no influence on the restoration displacement, and no influence on the change of volume after restoration.

  difference of restoration displacement between a restored model and the reference restored model (log10 scale, in cm)

Figure 3 .

 3 Figure 3.20: Impact of Poisson's ratio on the restoration displacement field (geomechanical restoration). Each curve corresponds to a restoration test with a different Poisson's ratio. The norm of the difference of the restoration displacement for each restoration test with the reference model (E = 70 GPa, ν = 0.2) are in the first row. The dilatation in percent for each restoration test relatively to the reference model (E = 70 GPa, ν = 0.2) are in the second row. Poisson's ratio has a low but significant impact on the restoration displacement, and is responsible on the change of volume after restoration.

Figure 3 . 21 :

 321 Figure 3.21: Impact of heterogeneous mechanical properties in the geomechanical restorations of the structural sandbox. Two mechanical models are tested (Table3.4). The results are for the first three restoration steps. The distributions correspond to the norm of the difference of restoration displacements with the reference restored models.

  3.16 Fault compliance distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Crossing faults restoration issue . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Problem of no relation between fault parts (offset fault) in the sandbox model 3.19 Impact of Young's modulus on the restoration displacement field (geomechanical restoration) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Impact of Poisson's ratio on the restoration displacement field (geomechanical restoration) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Impact of heterogeneous mechanical properties in the geomechanical restora-Restoration numerical method comparison . . . . . . . . . . . . . . . . . . . . 1.2 Comparison of the finite element methods used in geomechanical restoration . 1.3 Summary of the fault contact methods . . . . . . . . . . . . . . . . . . . . . . 1.4 Mathematical symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Mechanical characteristics of the materials used in the structural sandbox experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Characteristics of the dip slip delta distributions . . . . . . . . . . . . . . . . . 2.3 Incremental shortenings obtained by different methods . . . . . . . . . . . . . 2.4 Incremental forward dilatation of the sand and pyrex strata measured on the CT images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

tions of the structural sandbox . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Tables 1.1 3.1 Total W/E forward analog model extensions for each restoration steps . . . . . 3.2 Mechanical parameters of different geomechanical restorations . . . . . . . . . 3.3 Mechanical parameters of different geomechanical restorations . . . . . . . . . 3.4 Elastic parameters for the study of the impact of heterogeneous mechanical properties in restorations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Total restoration dilatation for each restoration step of the GeoChron-based method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Comparaison entre la restauration géomécanique et la restauration basée GeoChron : application sur le modèle analogique structural aux approches géomécaniques. Nous avons restauré géomécaniquement le modèle analogique en mettant des propriétés mécaniques hétérogènes. Les résultats obtenus montrent que le déplacement nodal diffère significativement même si l'écart n'est pas très grand. Par contre, le changement de volume est important. Il semblerait que le module de Young ait peu d'impact sur le modèle restauré, sûrement à cause des nombreuses conditions aux limites en déplacement. En ce qui concerne le coefficient de Poisson, la dilatation est fortement influencée localement et globalement par ce paramètre élastique.

xxiv Section 3.

  Geometric fault contact in the static method . . . . . . . . . . . . Geomechanical fault contact methods . . . . . . . . . . . . . . . . . . . . . 1.4.2.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2.2 Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Einstein notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A.3 Tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A.4 Transpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A.5 Matrix product . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A.6 Scalar product . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A.7 Product between a vector and a matrix . . . . . . . . . . . . . . . . . . . . 1.A.8 Double dot product . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A.9 Gradient of a vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	). O. C. Zienkiewicz and R. L. Taylor. The finite element method, volume 1, the basis. Butterworth-Heinemann, Oxford, United Kingdom, 5th edition, 2000a. O. C. Zienkiewicz and R. L. Taylor. The finite element method, volume 2, solid mechanics. Butterworth-Heinemann, Oxford, United Kingdom, 5th edition, 2000b. Chapter 1 Mechanics-based restoration: review of the different methods Contents Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Mechanics-based restoration: toward a restoration process including rock mechanical behavior . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 A restoration based on continuum mechanics . . . . . . . . . . . . . . . . . 1.1.1.1 Equations of motion: mass and linear momentum conservations . 1.1.1.2 A restoration based on elasticity assumption . . . . . . . . . . . . 1.1.2 Boundary conditions to ensure a geologically consistent restored state . . . 1.1.2.1 Displacement boundary conditions: unfolding and numerical stability 1.1.2.2 Unfaulting based on geometrical and/or physical boundary conditions 1.1.3 Numerical representation of a geological model in restoration . . . . . . . . 1.2 Numerical methods for solving a geomechanical problem . . . . . . . . 1.2.1 Geomechanical restoration solved by the finite element method . . . . . . . 1.2.1.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1.2 Application of the finite element method to the mechanics-based restoration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Geomechanical restoration solved by the boundary element method . . . . 1.2.3 Geomechanical restoration solved by the mass-spring method . . . . . . . . 1.2.4 Comparison between the different numerical methods . . . . . . . . . . . . 1.3 Geomechanical restoration solved by different finite element methods . 1.3.1 Static finite element method . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1.1 Global approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1.2 Local approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1.3 Application of boundary conditions: displacement, traction and body forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Dynamic relaxation method . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2.2 Definition of the Dirichlet boundary conditions . . . . . . . . . . . 1.3.3 Conclusions on the various finite element approaches . . . . . . . . . . . . . 1.4 Handling faults in geomechanical restoration solved by the finite ele-ment method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Geometric fault contact methods . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1.3 Geometric fault contact in the dynamic relaxation method . . . . 1.4.2 1.4.3 Some considerations on the fault contacts . . . . . . . . . . . . . . . . . . . 1.4.3.1 Master/slave choice . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.3.2 Physical consistency of the fault contacts . . . . . . . . . . . . . . 1.4.4 Comparison between the different methods . . . . . . . . . . . . . . . . . . Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A.1 Mathematical symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A.2 1.A.10 Gradient of a scalar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A.11 Divergence of a matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.B Elastic stress/strain laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.B.1 Hooke's law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.B.2 Stress tensor conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.B.3 Neo-Hookean law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.C Isotropic materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.C.1 Elastic parameter simplification in the isotropic case . . . . . . . . . . . . . 1.C.2 Conversion of the isotropic elastic parameters . . . . . . . . . . . . . . . . . 1.C.3 Isotropic elasticity tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.D Transverse isotropic materials used to reach flexural slip mode . . . . . 1.D.1 Transverse isotropic materials . . . . . . . . . . . . . . . . . . . . . . . . . . 1.D.2 Material upscaling to reach flexural slip mode . . . . . . . . . . . . . . . . . 1.E From current configuration to initial configuration . . . . . . . . . . . . 1.E.1 Volumetric integral conversion . . . . . . . . . . . . . . . . . . . . . . . . . 1.E.2 Mass conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.E.3 Stress conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.E.4 Equation of motion in the initial space . . . . . . . . . . . . . . . . . . . . . 1.F Integration by parts and Gauss theorem . . . . . . . . . . . . . . . . . . 1.4.1.2 Chapter 1. Mechanics-based restoration review 1.G Dynamic relaxation formula demonstration . . . . . . . . . . . . . . . . .

Table 1 .

 1

	1. Mechanics-based restoration review
	Criteria	methods Numerical	method element Finite	method element Boundary	method Mass-spring
	mechanics Use continuum	Yes		set of springs Simplified by a
		Mesh	mesh Volumetric	representation Boundary	volumetric mesh representation or Boundary
	constraints Meshing	High	Low	the mesh type Depends on
	time Computational	High		Low
	Formulation	3D	2.5D	Multiple 1D
	Mechanical
	anisotropy heterogeneity /	Easy	Difficult
	integration
	Memory cost	High	Low	the mesh type Depends on
			Muron,
	Authors	Moretti et al., Maerten and Maerten,	Gjerde	Midland Valley
			Tang et al.

1: Restoration numerical method comparison. The presented comparison is based on the published methods used for geomechanical restoration purposes. More details on the finite element method are in

Section 1.3

Table 1 . 2

 12

		Finite element approach
		Static		Dynamic
		Local	Global	relaxation
	Resolution type	equations (matrix) System of	Recurrence relation
	Temporal resolution		No
	integration Non-linearity		Yes
	for uniqueness Necessity of fixed dofs	No		Yes
	Boundary condition	Direct		Incremental
	Computational time	High		Very high
	Memory requirements	Low	High	Low
	Solver stability	High		Low
	Authors	Maerten and Maerten	Muron Moretti et al.,	Muron Santi et al.,

: Comparison of the finite element methods used in geomechanical restoration. This comparison is based on the pusblished methods used for geomechanical restoration purposes: static methods and the dynamic relaxation method.

Table 1 .

 1 .3. 3: Summary of the fault contact methods. This table only relates to the published methods of fault contact conditions in geomechanical restorations based on a finite element (FE) solver. See the text for more details.

	Approach
	Geometrical	Mechanical

Table 1 . 4

 14

	Chapter 1. Mechanics-based restoration review
	Symbol	Definition
	σ	Cauchy stress tensor/matrix
		Linear Green-Lagrange strain tensor/matrix
	P	Nominal stress tensor/matrix
	P T	First Piola-Kirchhoff tensor/matrix
	S	Second Piola-Kirchhoff tensor/matrix
	D	4th-order elasticity tensor
	λ	First Lamé parameter
	µ	Second Lamé parameter
	E	Young's modulus
	ν	Poisson's ratio
	u	Displacement vector
	u	Velocity vector
	ü	Acceleration vector
	t	Time
	∇	Gradient operator
	∇•	Divergence operator
	M	Mass matrix
	K	Stiffness matrix
	C	Damping matrix
	F int	Internal force vector
	F ext	External force vector
	b	Body forces
	ρ	Density
	K	Bulk modulus

: Mathematical symbols. This table lists and defines the mathematical symbols used in this chapter. See text within Chapter 1 for more details.

Table 2 .

 2 1: Mechanical characteristics of the materials used in the structural sandbox experiment. The relative density (unitless) of a material corresponds to the density of this material divided by the density of water. The Hounsfield density (in HU, no SI equivalence) is a measure of the X-ray attenuation in a medium. Data from C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System and IFPEN documentation.

		behavior Rheological	density Relative	size (µm) Grain	angle ( • ) friction Internal	(µPa) Cohesion	or cps) (mPa.s Viscosity	(HU) density Hounsfield
	Sand	Brittle	1.3-1.5 (0.004 in) 100	40	psi) 1-2 (10 -10 -2 × 10 -10		500
	Pyrex	Brittle	1.2	(0.004 in) 100	32-36	psi) > 5 (5 × 10 -10		150
	SGM36 Silicon	Ductile	0.97				5.10 7	95

  2.9-2.12), yielding a restoration of more than half of the growth strata interval. Restorations with a shortening boundary condition are in Figures 2.9E, 2.10E, 2.11E, and 2.12E. As we had the CT images of the paleo-states of the northern wall, we evaluated the shortenings by following a marker on the eastern wall. The measured, incremental shortenings for each restoration step are: 1.44 cm (0.567 in), 0.71 cm (0.28 in), 1.85 cm (0.728 in) and 1.85 cm (0.728 in). Qualitatively, the general consistency on the northern edge between the restorations with prescribed shortening and the reference CT images is quite good, indicating a robust and accurate restoration. Restorations without a shortening boundary condition are shown in Figures 2.9F, 2.10F, 2.11F, and 2.12F. In these models, we only avoided the shortening boundary condition for the last restoration step. For example, the result shown in Figure2.9F was not used as the starting model for the restoration in Figure2.10F. The starting model for Figure2.10F was generated by removing the restored, uppermost layer from the model shown in Figure2.9E. In this way, we attempt to avoid propagating errors. For each restoration step that does not include the shortening boundary condition (Figures 2.9F, 2.10F, 2.11F, and 2.12F), it is clear that there was not enough extension restored to be considered an acceptable restoration result. In each case, the restored faults are too far down-dip relative to the reference position obtained from the CT tomography video. In contrast, the restorations that included the shortening boundary condition provide a better qualitative match between the restored models and the reference CT images (Figures 2.9E, 2.10E, 2.11E, and 2.12E). For our model with the prescribed material properties, the classical boundary conditions alone are not sufficient to produce a reasonable restoration result.

Table 2 .

 2 .2, see discussions). 2: Characteristics of the dip slip delta distributions.

	Section 2.4. Estimation of shortening

Table 2 .

 2 , Marrett and 3: Incremental shortenings obtained by different methods. For each restoration step, several geometric methods are used, in addition to the measure on CT images, to assess the shortening magnitude (see text for details). Each percentage is relative to the corresponding reference shortening measured on the CT image.

					Horizon	H8	H7	H6	H5
	Incremental	shortenings in cm (in, %)	estimated by	different methods	CT image Area-depth (0.5717, 101%) (0.567, 100%) 1.44 1.452 heave sum Fault (0.266, 47%) 0.676 conservation Bed length (0.235, 41%) 0.597 condition (0.283, 50%) shortening No imposed 0.718	(0.28, 100%) 0.71 (0.317, 114%) 0.806 (0.172, 61%) 0.436 (0.101, 36%) 0.257 (0.136, 49%) 0.345	(0.728, 100%) 1.85 (0.7114, 98%) 1.807 (0.4157, 57%) 1.056 (0.363, 50%) 0.922 (0.328, 45%) 0.834	(0.728, 100%) 1.85 (0.6354, 87%) 1.614 (0.237, 33%) 0.602 (0.195, 27%) 0.495 (0.181, 25%) 0.459

Table 3 .

 3 Comparison between mechanics-based and GeoChron-based restorations 4: Elastic parameters for the study of the impact of heterogeneous mechanical properties in restorations. Two rock configurations are used for the structural sandbox restorations.

		Models	Alternating layers	Big layers
	Layers	Mechanical parameters E (GPa)	ν	E (GPa) ν
		L7	10	0.2	70	0.4
		L6	70	0.4	70	0.4
		L5	10	0.2	70	0.4
		L4	70	0.4	70	0.4
		L3	10	0.2	10	0.2
		L2	70	0.4	10	0.2
		L1	10	0.2	10	0.2

http://www.ifpenergiesnouvelles.fr

http://www.ccreservoirs.com

http://www.ring-team.org/software/ring-libraries/44-ringmecha

http://www.ring-team.org/software/skua-gocad-plugins/38-restorationlab

the successor of RestorationLab

also called effective tangent stiffness matrix[Belytschko et al., 2013] or effective matrix[Jacob and Ebecken, 1994]

http://www.ring-team.org/software/skua-gocad-plugins/38-restorationlab

http://www.ring-team.org/software/ring-libraries/44-ringmecha

http://www.ring-team.org/software/ringmesh

Acknowledgements

For non-French people Acknowledgements Acknowledgements Acknowledgements

Ackowledgements

This work was performed as part of the RING project at Université de Lorraine. We would like to thank the industrial and academic sponsors of the RING-GOCAD Consortium managed by ASGA for their support, and Chevron for funding the Ph.D. of Benjamin Chauvin. We also acknowledge Paradigm for the SKUA-GOCAD software and API, and Inria for the Geogram library used in RINGMesh. We thank IFPEN and C&C Reservoirs, DAKS TM -Digital Analog Knowledge System, for the analog model data set. We also acknowledge Justin Herbert, Donald Medwedeff and Richard Groshong for discussions, feedback and involvement in this work. Frantz Maerten, Robert Worthington and Frank Zwaan provided helpful comments that improved the quality of this manuscript.

Acknowledgements

This work was performed in the frame of the RING project at Université de Lorraine. We would like to thank for their support the industrial and academic sponsors of the RING-GOCAD Consortium managed by ASGA, and Chevron for funding the PhD of Benjamin Chauvin. We also acknowledge Paradigm for the SKUA-GOCAD software and API. We thank IFPEN and C&C Reservoirs, DAKS TM -Digital Analogs Knowledge System, for the analog model data set and Chevron for the results of the GeoChron restorations on the structural sandbox model.

Contributions

The different contributions of this thesis are:

• A review on the different geomechanical restoration methods.

• The requirement of a lateral boundary condition (shortening) to properly restore an extensional structural model.

• The estimation of the shortening boundary condition magnitude using the area-depth method or the restoration dilatation.

Authors: Benjamin P. Chauvin 1 , Peter J. Lovely 2 , Joseph M. Stockmeyer 3 *, Andreas Plesch 3 , Guillaume Caumon 1 , John H. Shaw 3

General conclusions

Contributions of this thesis

This thesis provided several keys to the current issues of the mechanics-based restoration. Our main contributions are related to the choice of the boundary conditions and on the comparison of restoration methods. We also provided insight on the choice of elastic parameters.

Definition of proper boundary conditions

We tested different sets of boundary conditions in the restoration of an extensional structural sandbox model. Thanks to computed tomography (CT), the deformation path was known on a cross section. We could objectively determine which boundary conditions geometrically provided the best restored state in comparison to the CT images.

Lateral displacement boundary condition

We showed that a restoration with the classical boundary conditions, i.e., a datuming and fault contact boundary conditions, did not permit to properly capture the forward extension of the analog model. The use of a shortening boundary condition was essential to fit the restored geometries with the CT images. The magnitude of this shortening condition corresponded to the forward displacement observed on the CT images. These conclusions are consistent with others studies made on synthetic models by Durand-Riard [2010], Lovely et al. [2012] and Durand-Riard et al. [2013], respectively on a compressive, an extensive and a strike-/obliqueslip contexts. In our case, the magnitude of the displacement condition to apply was deduced from the CT images. In natural case studies, it is not straightforward to get such information.

We showed that the area-depth method [e.g., Epard andGroshong, 1993, Groshong et al., 2003] provides a good estimate of the shortening magnitude. We also proposed a tetrahedral dilatation analysis to assess this magnitude and to check the validity of the restored state [Durand-Riard, 2010].

Contact conditions to handle complex fault networks

We defined two novel contact boundary conditions to handle the complex fault network present in the structural analog model. The first condition connects fault surface internal borders to handle faults that branch onto other faults. The second condition ensures the continuity of an offset fault. These conditions use the same contact method to tie fault cutoff lines by contacts mechanics [Muron, 2005, Maerten and Maerten, 2006, Wriggers and Laursen, 2006]. Thanks to these conditions, several steps of the analog model sequential restoration were possible.

Impact of the elastic parameters

The mechanics-based restoration uses elastic parameters, i.e., Young's modulus and Poisson's ratio, to take into account rock mechanical behavior. We studied the impact of these parameters on the restored state of the analog model.

Impact of Young's modulus

We showed that on a homogeneous case, Young's modulus has a small impact on the displacement field. We justify that by the fact that all the boundary conditions are displacement conditions. As the strain is highly constrained by these conditions, Young's modulus only changes the stress in consequence. Althought the displacement field is not significantly affected, the Young's modulus still has an important impact on the stress field. An erroneous stress may lead to erroneous interpretations such as a bad prediction of fracture areas [Macé, 2006, Maerten and Maerten, 2006, Mejía-Herrera et al., 2014, Stockmeyer et al., 2017].

Impact of Poisson's ratio

We showed on a homogeneous case that Poisson's ratio has more significant impact on the displacement field than Young's modulus. Even if the nodal displacement does not seem to be highly modified when Poisson's ratio is changed, the subsequent change of volume, locally and globally, is very important. A dilatation analysis must be performed to reject Poisson's ratios which provide an irrealistic dilatation.

Comparison of restoration methods

A multitude of restoration methods exists, defined by geometric, kinematic and/or geomechanical rules. The choice of a specific method of restoration is not always straightforward. So we reviewed the existing geomechanical restoration methods and we did an experimental comparison between our geomechanical method [Muron, 2005] and the geometrical approach relying on the Geo-Chronological model [Mallet, 2014, Medwedeff et al., 2016].

Comparison between the geomechanical approaches

We theoretically compared the different geomechanical approaches. All these methods are based on fundamental equations of continuum mechanics. We investigated the differences on the numerical methods used to solve a restoration problem: the mass-spring method [Terzopoulos et al., 1987, Macaulay et al., 2015], the boundary element method [Gjerde, 2002], and the finite element method [e.g., Zienkiewicz and Taylor, 2000a,b, Belytschko et al., 2013, Bathe, 2014]. For the latter method, several approaches exist, all neglecting the temporal part of the resolution to focus on the steady-state: the static method solved locally [Maerten and Maerten, 2006], the static method solved globally [Muron, 2005, Moretti et al., 2006], and the dynamic relaxation method [Santi et al., 2003, Muron, 2005]. We also explained the different approaches to handle fault contacts: by a geometrical method [Muron, 2005, Moretti et al., 2006] or by a mechanical method [Muron, 2005, Maerten and Maerten, 2006, Wriggers and Laursen, 2006, Tang et al., 2016]. All these approaches have the same purpose but each has geometrical and physical implications on the restored state. The choice of a specific method is dependent on the needs and constraints of the structural geologist and of the case study.

Comparison with a geometric method based on GeoChron

We geometrically compared the analog model restored states provided by a geomechanical restoration approach [Muron, 2005, Chauvin et al., 2017] and the GeoChron-based restoration [Mallet, 2014, Medwedeff et al., 2016]. We showed that both restoration methods provided very similar paleo-geometries. We attribute this simulitude to the mechanical simplicity of the structural sandbox model, i.e., no mechanical contrast and a deformation process simpler than in nature, and to the deformation minimization in the GeoChron-based method. We also showed that, as the geomechanical restoration, the GeoChron-based restoration does not properly capture the forward extension of the analog model. Dilatation, and extension due to missing faults which are below the seismic data resolution, are not taken into account [Kautz and Sclater, 1988, Marrett and Allmendinger, 1992, Baxter, 1998, Groshong et al., 2003].

A scaling post-process is needed on the GeoChron restored states, similarly to the shortening boundary condition in the mechanics-based restoration. The geomechanical restoration approach has the main advantage to define custom boundary conditions and rock elastic behaviors. Conversely, the GeoChron-based approach has almost no flexibity and assumes the knownledge of a tectonic style applied globally. However, this restoration method is very interesting since it quickly provides the first order paleo-geometries and highly reduces several issues met in the geomechanical approach such as the meshing constraints.

Publications associated to this thesis During this thesis, I supervised several Master students in different topics related to the restoration in general. These projects were not directly related to this thesis. In total, we supervised 9 projects:

• A new method to reconstruct eroded paleotopographies using mass balance principle, Gabriel Godefroy (Master student, 2014).

• Geometrical errors after restoration, Claire Launoy (Master student, 2015).

• Integration of physical parameters into decompaction during the restoration, Simon Fuet (Master student, 2015).