Alison Malcolm

Thomas Bohlen

Alexandrine Gesret

Isabelle Lecomte

Guillaume Caumon

Paul Cupillard

Paul Sava

for accepting the invitation to evaluate my work. Isabelle Lecomte also reviewed chapter III.2 when it

My list of remarkable teachers would not be complete without mentioning Pauline Collon, who taught me C++ programming. I think the most important skill I learned from the ENSG is probably C++ programming. Thank you Pauline! Pauline, thank you also for my final year internship at CEREMA. There is also François Bonneau, who I have known since I was a student at the ENSG. François has always been available to discuss both scientific and non-scientific matters; I feel like he has also supervised me in some sense during my stay in Nancy. When I was still a student at the ENSG, I did a summer internship at the Laboratoire de Géologie de Lyon (thanks to Paul Cupillard) where I was supervised by Benoit Tauzin, who taught me a lot about inverse problems. Over the course of my PhD, I had the opportunity to supervise three master students: Jean Langanay, Thomas Beraud, and Morgan Thierry-Coudon; I learned as much from them as they learned from me. Thomas Beraud made significant contributions to chapter I.3 and Morgan Thierry-Coudon made significant contributions to chapter II.2. Other colleagues who made significant contributions to this manuscript include Gautier Laurent (chapter II.1), Julien Renaudeau (chapter II.1. I also got the chance to collaborate with Julien on one of his papers), Mustapha Zakari (chapter II.2), and Jonathan Edwards (chapter III.2). Melchior Schuh-Senlis also made important contributions by proofreading part of this manuscript and proofchecking my algebra in chapters I.2 and I.3. There are many more colleagues who have contributed to this manuscript, some of whom are mentioned in the ac-knowledgements section of relevant chapters throughout the manuscript. I'm very grateful to my PhD advisory committee composed of Ludovic Métivier and Pierre Thore. I got stuck a number of times during my PhD and discussing with Ludovic Métivier and Pierre Thore was often all that I needed to make some progress. Equally important, Ludovic Métivier and Pierre Thore pointed out many flaws in my reasoning, thus enabling me to propose less flawed arguments.

I would like to thank the following members of the Association Scientifique pour la Géologie et ses Applications -ASGA -: Armelle Faynot, Christophe Antoine, Fatima Chtioui, Maryse Villeneuve, Nicole Galloy and Sophie Romain, for making my life in Nancy much easier (and consequently happier). I also acknowledge the RING-GOCAD Consortium, managed by ASGA, for funding my PhD. I enjoyed my stay in Nancy greatly due to my thoughtful colleagues in the RING-team (in the following list I omit those already mentioned above) : Antoine Mazuyer, Arnaud Botella, Benjamin Chauvin, Corentin Gouache, Gabriel Godefroy, Margaux Raguenel, Marion Parquer, Nicolas Clausolles, Nicolas Mastio, Paul Baville, Pierre Anquez (whith whom I got an opportunity to collaborate on a paper), and Yves Frantz. Thanks guys: I'm really lucky, thankful, and happy to have been part of your team! Last but not least, I would like to thank my family: my parents, two sisters, one brother, and my newly born niece, for their support; they have ALWAYS been supportive and I would not have finished my master's degree at the ENSG had it not been for their support. Speaking of everlasting support, I would also like to thank my lifelong friends, who are now practically family, Ralph P. and Michel N.

"Restons factuels les mecs!" -Guillaume Caumon "Rien n'est simple!" -Alfred B. (a.k.a my dad) I would like to dedicate this thesis to my mother who was convinced I had what it took from the beginning (before I started university). Thanks mum!

Résumé

Les modèles géologiques sont couramment utilisés pour estimer les ressources souterraines, pour faire des simulations numériques, et pour évaluer les risques naturels; il est donc important que les modèles géologiques représentent la géométrie des objets géologiques de façon précise. La première étape pour construire un modèle géologique consiste souvent à interpréter des surfaces structurales, telles que les failles et horizons, à partir d'une image sismique; les objets géologiques identifiés sont ensuite utilisés pour construire le modèle géologique par des méthodes d'interpolation. Les modèles géologiques construits de cette façon héritent donc les incertitudes d'interprétation car une image sismique peut souvent supporter plusieurs interprétations structurales. Dans ce manuscrit, j'étudie le problème de réduire les incertitudes d'interprtation à l'aide des données sismiques. Particulièrement, j'étudie le problème de déterminer, à l'aide des données sismiques, quels modèles sont plus probables que d'autres dans un ensemble des modèles géologiques cohérents. Ce problème sera connu par la suite comme le problème d'évaluation des modèles géologiques par données sismiques.

La première partie du manuscrit est consacrée à l'imagerie sismique. J'y propose de se servir de la "reverse time migration" (RTM) pour préconditionner l'inversion des formes d'ondes. Les expériences numériques montrent que le préconditionneur proposé accélère l'inversion des formes d'ondes linéarisée ("least squares reverse time migration") et l'inversion des formes d'ondes non-linéaire ("full waveform invesion") d'au moins un ordre de grandeur. Je justifie la performance numérique positive du préconditionneur proposé en montrant algébriquement qu'un filtre passe-bas d'une image RTM peut approximer les éléments de la diagonale de la matrice hessienne de la fonction-objectif sous des hypothèses appropriées. Cependant, je ne suis toujours pas en mesure de proposer une interprétation physique du filtre passe-bas et de la manière dont il relie l'image RTM aux éléments de la diagonale de la matrice hessienne. Ensuite, je propose un opérateur généralisé du "Kirchhoff extended imaging" pour la construction de modèles de vitesses; l'opérateur est généralisé dans le sens où il décrit simultanément plusieurs extensions dans le "data-domain" (ex. extension en source, extension en distance source-récepteur, extension en angle de réflexion) et des extensions dans l'"image-domain" (ex. extension en temps, extension en décalage spatial). Il existe au moins deux avantages de l'opérateur généralisé proposé: premièrement, il permet une implémentation unifiée pour plusieurs extensions (c'est-à-dire une seule implémentation valable pour plusieurs extensions); deuxièmement, lopérateur aboutit à un algorithme du gradient unifié de la "migration velocity analysis" (MVA). J'illustre la sensibilité de l'opérateur généralisé proposé aux vitesses imprécises en utilisant cet opérateur dans une expérience MVA basée sur des rayons sismiques.

La deuxième partie du manuscrit est consacrée à la modélisation structurale, en particulier à l'interpolation structurale implicite. J'y propose une nouvelle méthode sous le nom de "Finite Difference Structural Implicit Modeling" (FDSIM). Il existe au moins deux avantages de la FD-SIM: premièrement, il est relativement facile d'implémenter cette méthode sur ordinateur et de l'optimiser car elle est basée sur les différences finies sur des grilles régulières; deuxièmement, comme la FDSIM gère les discontinuités par rastérisation, cette méthode s'est montrée capable de gérer facilement des réseaux de failles très complexes. Le principal inconvénient de la méthode est qu'elle peut nécessiter une résolution très fine en fonction de la complexité du réseau de failles, conduisant parfois à des limites de mémoire. Je propose également de nouveaux opérateurs de régularisation; la particularité de ces opérateurs est qu'ils n'ont pas besoin dêtre implémentés sur les noeuds aux bords, une propriété très convenable pour la modélisation structurale implicite où les conditions aux limites sont généralement inconnues. Je propose ensuite une autre méthode sous le nom de "Finite Element Structural Implicit Modeling" (FESIM). La FESIM est basée sur une implémentation éléments finis des opérateurs de régularisation dernièrement proposés. Je montre que la méthode des éléments finis classique connue pour résoudre les problèmes aux limites doit être légèrement modifiée pour la modélisation implicite où les conditions aux limites sont généralement inconnues.

La troisième partie du manuscrit est consacrée à l'évaluation des modèles structuraux à l'aide des données sismiques. J'y introduis et formalise le problème d'évaluation des interprétations structurales à l'aide des données sismiques. Je propose de résoudre ce problème par génération des données sismiques synthétiques pour chaque interprétation structurale dans un premier temps, ensuite d'utiliser ces données synthétiques pour calculer la fonction-objectif pour chaque interprétation; cela permet de classer les différentes interprétations structurales. La difficulté majeure d'évaluer les modèles structuraux à l'aide des données sismiques consiste à proposer des fonctions-objectifs adéquates. Je propose un ensemble de conditions qui doivent être satisfaites par la fonction-objectif pour une évaluation réussie des modèles structuraux à l'aide des données sismiques. Je montre que, puisqu'il n'est pas possible de satisfaire ces conditions en utilisant les données de type "vertical seismic profile" (VSP), il n'est pas possible d'évaluer les interprétations structurales à l'aide des données VSP dans le cas le plus général. Ces conditions imposées à la fonction-objectif peuvent en principe être satisfaites en utilisant les données sismiques de surface ("surface seismic data"). Cependant, en pratique il reste tout de même difficile de proposer et de calculer des fonctions-objectifs qui satisfassent ces conditions. Je termine le manuscrit en illustrant les difficultés rencontrées en pratique lorsque nous cherchons à évaluer les interprétations structurales à l'aide des donnés sismiques de surface. Je propose une fonction-objectif générale faite de deux composants principaux: (1) un opérateur de résidus qui calcule les résidus des données (c'est-à-dire la différence entre données synthétiques et données observées), et (2) un opérateur de projection qui projette les résidus de données depuis l'espace de données vers l'espace physique (c'est-à-dire le sous-sol). Cette fonction-objectif est donc localisée dans l'espace car elle génère des valeurs en fonction de l'espace. Cependant, je ne suis toujours pas en mesure de proposer une implémentation pratique de cette fonction-objectif qui satisfasse les conditions imposées pour une évaluation réussie des interprétations structurales; cela reste un sujet de recherche.

Abstract

Subsurface structural models are routinely used for resource estimation, numerical simulations, and risk management; it is therefore important that subsurface models represent the geometry of geological objects accurately. The first step in building a subsurface model is usually to interpret structural features, such as faults and horizons, from a seismic image; the identified structural features are then used to build a subsurface model using interpolation methods. Subsurface models built this way therefore inherit interpretation uncertainties since a single seismic image often supports multiple structural interpretations. In this manuscript, I study the problem of reducing interpretation uncertainties using seismic data. In particular, I study the problem of using seismic data to determine which structural models are more likely than others in an ensemble of geologically plausible structural models. I refer to this problem as appraising structural models using seismic data.

The first Part of the manuscript is devoted to seismic imaging. I first propose to use reversetime migration (RTM) as a preconditioner for waveform inversion. Numerical experiments show that the proposed preconditioner accelerates both linearized waveform inversion (least squares reverse time migration) and nonlinear waveform inversion (full waveform inversion) by at least an order of magnitude. I justify the positive numerical performance of the proposed preconditioner by showing algebraically that a low pass filter of the RTM image can approximate the diagonal elements of the Hessian matrix of the objective function under appropriate assumptions. However, I am still unable to propose a physical meaning of the low pass filtering and how it relates the RTM image to the elements of the diagonal of the Hessian matrix. Then, I propose a generalized extended Kirchhoff imaging operator for velocity modeling; the operator is generalized in the sense that it describes multiple data-domain extensions (e.g. shot, offset, and angle extensions) and image-domain extensions (e.g. time-lag and space-lag extensions) simultaneously. The advantages of the proposed generalized extended operator are twofold: firstly, it allows a unified implementation for multiple extensions (i.e. a single implementation that is valid for multiple extensions); secondly, the operator leads to a unified gradient-based migration velocity analysis (MVA) scheme. I confirm the ability of the proposed generalized extended operator to capture image distortion caused by inaccurate velocity by applying it to a ray-based MVA experiment.

The second Part of the manuscript is devoted to structural modeling, particularly implicit structural interpolation. I introduce Finite Difference Structural Implicit Modeling (FDSIM). The advantages of FDSIM are twofold: firstly, it is relatively easy to implement and to optimize since it is based on finite differences on regular grids; secondly, because it handles discontinuities by rasterization, FDSIM has shown to easily handle very complex fault networks. The main disadvantage of the method is that it may require a very fine resolution depending on the complexity of the fault network, sometimes leading to memory limits. I also propose new regularization operators; the particularity of these operators is that they do not need to be implemented on boundary nodes, a property which is very appealing in implicit modeling where boundary conditions are usually unknown. I then introduce Finite Element Structural Implicit Modeling (FESIM). FESIM is based on a finite element implementation of the newly proposed regularization operators. I show that the conventional finite element familiar for solving boundary value problems has to be slightly modified for implicit modeling where boundary conditions are usually unknown.

The third Part of the manuscript is devoted to appraising structural models/interpretations using seismic data. I introduce and formalize the problem of appraising structural interpretations using seismic data. I propose to solve the problem by generating synthetic data for each structural interpretation and then to compute misfit values for each interpretation; this allows us to rank the different structural interpretations. The main challenge of appraising structural models using seismic data is to propose appropriate data misfit functions. I derive a set of conditions that have to be satisfied by the data misfit function for a successful appraisal of structural models. I argue that since it is not possible to satisfy these conditions using vertical seismic profile (VSP) data, it is not possible to appraise structural interpretations using VSP data in the most general case. The conditions imposed on the data misfit function can in principle be satisfied for surface seismic data. In practice, however, it remains a challenge to propose and compute data misfit functions that satisfy those conditions. I conclude the manuscript by highlighting practical issues of appraising structural interpretations using surface seismic data. I propose a general data misfit function that is made of two main components: (1) a residual operator that computes data residuals, and (2) a projection operator that projects the data residuals from the data-space into the image-domain. This misfit function is therefore localized in space, as it outputs data misfit values in the image-domain. However, I am still unable to propose a practical implementation of this misfit function that satisfies the conditions imposed for a successful appraisal of structural interpretations; this is a subject for further research. La modélisation structurale est la pratique de construire des modèles géologiques du soussol. Un exemple d'un modèle géologique du sous-sol en 3D est illustré dans la Figure 1. Les modèles géologiques servent dans plusieurs applications stratégiques comme dans l'évaluation des risques naturels, ou encore dans l'estimation des ressources naturelles souterraines. Il est donc important que les modéles géologiques représentent la géométrie des structures géologiques souterraines de façon précise. Dans ce manuscrit, nous allons nous intéresser surtout à des domaines stratigraphiques dans lesquels les modèles structuraux sont souvent construits à partir des données interprétées sur des images sismiques; en particulier, nous allons nous intéresser à la modélisation structurale implicite (voir par exemple Caumon et al., 2013) telle qu'illustrée dans la Figure 2. Malheureusement, les modèles structuraux construits de la manière décrite précédemment peuvent être sujets à des incertitudes d'interprétation structurale (voir par exemple Bond et al., 2007) telles qu'illustrées dans la Figure 3. Dans ce manuscrit, nous allons étudier le problème de reduire ce type d'incertitudes à l'aide des données sismiques. En particulier, nous allons étudier le problème de se servir des données sismiques pour determiner quels modèles sont plus probables que d'autres dans un ensemble des modèles structuraux géologiquement cohérents; ce problème sera connu sous le nom du problème d'évaluation des modèles structuraux par données sismiques par la suite. Dans la première partie du manuscrit je propose deux contributions dans le domaine de l'imagerie sismique. Premièrement, dans le chapitre I.2, je propose d'utiliser la "reverse-time migration" (RTM) (voir par exemple Baysal et al., 1983, pour L'idée de base repose sur l'observation que la formule du "Kirchhoff imaging" (voir par exemple Etgen et al., 2009, pour

Résumé long

où β est le paramètre d'extension. Dans l'équation 1, I(x) désigne l'image sismique calculée, s la position de la source, r la position du récepteur, W (s, r, x) une fonction de pondération, T (s, r, x) le temps de trajet d'un rayon sismique issu de la source s puis réfléchi sur le point x et enfin enregistré au récepteur r, et D(s, r, t) les données sismiques. La fonction ζ(s, r, x) dans l'équation 2 est un paramètre libre permettant à l'utilsateur de définir la nature du paramètre d'extension β; par exemple, si on définit ζ(s, r, x) = |s -r|, l'image construite par l'équation 2 est par définition une "surface-offset extended image". L'équation 2 a le mérite d'aboutir à une formulation unifiée de la "migration velocity analysis" (MVA) (voir par exemple Woodward et al., 2008, pour plus d'informations sur la MVA), c.-à-d. l'équation aboutit potentiellement à une implémentation unique valable pour plusieurs types de MVA selon la définition de ζ(s, r, x) (voir la section I.3.6 pour plus de détails). Dans la deuxième partie du manuscrit je propose deux contributions dans le domaine de la modélisation structurale implicite. Premièrement, dans le chapitre II.1, je propose une nouvelle méthode pour la modélisation implicite par différences-finies. Cette méthode a le mérite de gérer facilement les réseaux complèxes de failles, comme l'illustre la Figure 5 , car la méthode gère les discontinuités par rastérisation; en plus, la méthode est simple à implémenter et optimiser. Dans ce travail, je propose aussi de nouveaux opérateurs de régularisation, par exemple l'opérateur

R(φ) =            ∂ 2 x φ ∂ 2 z φ 1 2 (∂ 2 x + ∂ 2 z -2∂ x ∂ z )φ 1 2 (∂ 2 x + ∂ 2 z + 2∂ x ∂ z )φ
en 2D (voir chapitre II.1 pour la 3D), où φ(x) est la fonction implicite qui nous intéresse. La particularité de ces opérateurs est qu'ils n'ont pas besoin d'être implémentés sur les noeuds aux bords, une propriété très convenable pour la modélisation structurale implicite où les conditions aux limites sont généralement inconnues. Deuxièmement, dans le chapitre II.2, je propose une autre méthode pour la modélisation implicite par éléments-finis. Cette méthode est basée sur l'implémentation éléments-finis de nouveaux opérateurs proposés dans le chapitre II.1. Je montre que la méthode des éléments finis classique connue pour résoudre les problèmes aux limites doit être légèrement modifiée pour la modélisation implicite où les conditions aux limites sont généralement inconnues.

Φ i (Ω) = x∈Ω E i (x) = x∈Ω P R(d i , d o ) , (3)
                     a i (s, r, t) = R d i (s, r, t), d o (s, r, t) = d i (s, r, t) -d o (s, r, t) E i (x) = P R(d i , d o ) = argmin ∆m(x) s,r,t R 2 a(s, r, t, ∆m), a i (s, r, t) = argmin ∆m(x) s,r,x,t δ(x -r)R 2 a(s, x, t, ∆m), a i (s, x, t) , subject to L(m + ∆m)a(x, t, ∆m; s) = δ(x -s)w(t) Φ i (Ω) = x∈Ω E i (x)
.

Dans cette équation, s désigne la position de l'ondelette source w(t), r la position des récepteurs, m le modèle initial de vitesses lisse, L l'équation d'onde, et R 2 un opérateur de résidus de type "phase-shift" (voir Ma and Hale, 2013;[START_REF] Warner | Adaptive waveform inversion: Theory[END_REF][START_REF] Yang | Recent advances in imaging crustal fault zones: a review[END_REF], pour quelques exemples des opérateurs de résidus de type "phase-shift").

Enfin, dans la quatrième partie du manuscrit, je présente brièvement SIGMA, un code développé pendant ma thèse. Tous les résultats présentés dans ce manuscrit sont reproductibles et ont été produits par SIGMA. Structural modeling is the practice of building subsurface geological models. An example of a 3D subsurface geological model is illustrated in Figure 7. Geological models are useful in a number of strategic applications such as the evaluation of natural risks and the estimation of subsurface natural resources. It is therefore important that geological models represent the geometry of subsurface geological structures accurately. In this manuscript, we will focus mainly on stratigraphic domains where structural models are often built from data interpreted from seismic images; in particular, we are going to focus on structural implicit modeling (see for example Caumon et al., 2013) as illustrated in Figure 8. Unfortunately, structural models built as just described can be subjected to structural interpretation uncertainties (see for example Bond et al., 2007) as illustrated in Figure 9. In this manuscript, we are going to study the problem of reducing these types of uncertainties using seismic data. In particular, we are going to study the problem of using seismic data to determine which models are more likely than others from a set of geologically consistent structural models; this problem will be referred to as the problem of appraising structural models using seismic data in what follows. In the first Part of the manuscript I propose two contributions in the field of seismic imaging. Firstly, in chapter I.2, I propose to use reverse time migration (RTM) (see for example Baysal et al., 1983, for more information on RTM) as a preconditioner to speed up waveform inversion (see for example Virieux and Operto, 2009, for more information on waveform inversion). The underlying idea is based on the observation that the RTM image gives a good estimation of the seismic illumination in the subsurface for a given acquisition geometry. This method has the advantage of being simple to implement while still remaining numerically efficient as illustrated in Figure 10. Secondly, in chapter I.3, I propose a generalized formula for Kirchhoff extending imaging (see for example Symes, 2008, for more information on the concept of extended imaging). The underlying idea is based on the observation that the standard Kirchhoff imaging formula (see for example Etgen et al., 2009, for more information on Kirchhoff imaging)

Extended Abstract

I(x) = s,r W (s, r, x) ∂ ∂t D s, r, T (s, r, x) = s,r F (s, r, x) (4)
can be extended by writing

I(x, β) = s,r δ β -ζ(s, r, x) F (s, r, x), ( 5
)
where β is the extension parameter. In Equation 4, I(x) designates the output seismic image, s the source location , r the receiver location, W (s, r, x) a weighting function , T (s, r, x) the travel time of a seismic ray from the source s then reflected off the point x and finally recorded at the receiver r, and D(s, r, t) the input seismic data. The function ζ(s, r, x) in Equation 5is a free parameter that allows the user to define the nature of the extension parameter β; for example, if one defines ζ(s, r, x) = |s -r|, the image volume output by Equation 5is by definition a surface-offset extended image. Equation 5 has the advantage of leading to a unified formulation of gradient-based migration velocity analysis (MVA) (see for example Woodward et al., 2008, for more information on MVA), i.e. the equation potentially leads to a unique implementation that is valid for multiple types of MVA depending on the definition of ζ(s, r, x) (see section I.3.6 for more details). In the second Part of the manuscript I propose two contributions in the field of structural implicit modeling. Firstly, in chapter II.1, I propose a new method for implicit modeling based on finite differences. This method has the advantage of readily handling complex fault networks, as illustrated in Figure 11, since the method handles discontinuities by rasterization; furthermore, the method is easy to implement and optimize. In this work, I also propose new regularization operators, for example the operator

R(φ) =            ∂ 2 x φ ∂ 2 z φ 1 2 (∂ 2 x + ∂ 2 z -2∂ x ∂ z )φ 1 2 (∂ 2 x + ∂ 2 z + 2∂ x ∂ z )φ
in 2D (see chapter II.1 for 3D), where φ(x) is the implicit function of interest. The particularity of these operators is that they do not need to be implemented on boundary nodes, a property which is very appealing in implicit modeling where boundary conditions are usually unknown. Secondly, in chapter II.2, I propose another method for implicit modeling based on finite elements. This method is merely a finite element implementation of the new regularization operators proposed in chapter II.1. I show that the conventional finite element familiar for solving boundary value problems has to be slightly modified for implicit modeling where boundary conditions are usually unknown. The third Part of the manuscript is devoted to appraising structural models/interpretations using seismic data. I propose two contributions in this field. Firstly, in chapter III.2, I develop a quantitative theory that formalizes the problem of appraising structural models using seismic data. The global workflow proposed to tackle this problem is illustrated in Figure 12: the main idea is to generate synthetic data for each structural interpretation, and then to compare these data to observed data by computing an appropriate data misfit function. I derive a set of conditions that have to be satisfied by the data misfit function for a successful appraisal of structural models, then I argue that a misfit function that is likely to satisfy those conditions has the form

Φ i (Ω) = x∈Ω E i (x) = x∈Ω P R(d i , d o ) , (6)
where E(x) designates the interpretation error map, R a residual operator that computes the difference between observed d o and the synthetic data d i computed in the i th structural model, P a projection operator that projects data residuals from the data-space into the imagedomain; Ω is the region we would like to evaluate in the subsurface. Secondly, in chapter III.3, I study the practical aspects of the problem and show that it remains a challenge to propose an implementation of Equation 6 that satisfies the conditions imposed for a successful appraisal of structural models. Although I fail to demonstrate a practical implementation that satisfies the required conditions, I conclude that the most promising implementation has the form

                     a i (s, r, t) = R d i (s, r, t), d o (s, r, t) = d i (s, r, t) -d o (s, r, t) E i (x) = P R(d i , d o ) = argmin ∆m(x) s,r,t R 2 a(s, r, t, ∆m), a i (s, r, t) = argmin ∆m(x) s,r,x,t δ(x -r)R 2 a(s, x, t, ∆m), a i (s, x, t) , subject to L(m + ∆m)a(x, t, ∆m; s) = δ(x -s)w(t) Φ i (Ω) = x∈Ω E i (x)
.

In this equation, s designates the location of the source wavelet w(t), r the location of receivers, m the initial smooth imaging velocity model, L the wave equation, and R 2 a phase-shit residual operator (see Ma and Hale, 2013;[START_REF] Warner | Adaptive waveform inversion: Theory[END_REF][START_REF] Yang | Recent advances in imaging crustal fault zones: a review[END_REF] for some examples of phase-shift residual operators).

Finally, in the fourth Part of the manuscript, I briefly present SIGMA, a software developed during my PhD. All the results presented in this manuscript are reproducible and were obtained using SIGMA.

Introduction

Structural modeling is the practice of building subsurface structural models. A structural model is a set geologically consistent structural surfaces such as faults and horizons (e.g. Caumon et al., 2009); an example of a structural model is shown in Figure 13 Subsurface structural models are routinely used for resource estimation, numerical simulations, and risk management; it is therefore important that structural models represent the geometry of subsurface geological objects accurately. In this manuscript, we will only discuss structural models that are built from input data interpreted from seismic images, as shown in Figure 14. Seismic images are obtained by seismic imaging. Seismic imaging is the practice of using of elastic waves propagated in the subsurface to image geological structures of the subsurface. A typical seismic imaging workflow is illustrated in Figure 15. Seismic waves are first propagated in the subsurface (e.g. Figure 15a) and then recorded at the surface of the Earth as seismic data (e.g. Figure 15b) after scattering off heterogeneities in the subsurface. Recorded data are then used to estimate a migration velocity model (e.g. Figure 15c) which is used to migrate the observed data to obtain a seismic image of the subsurface (e.g. Figure 15d). The reader interested in learning more on seismic imaging is referred to Baysal et al. (1983); Claerbout (1985); Miller et al. (1987); Bleistein et al. (2001); Xu et al. (2001); [START_REF] Sava | Angle-domain common-image gathers by wavefield continuation methods[END_REF]; Biondi (2006); Plessix (2006); Etgen et al. (2009). Structural models built from seismic interpretations may be subjected to (structural) uncertainties inherited from seismic imaging. We speak of structural uncertainties when the position and/or the presence of structural surfaces is uncertain in the subsurface (e.g. Thore et al., 2002;Wellmann et al., 2010). We can distinguish two types of structural uncertainties inherited from seismic imaging 1. Structural uncertainties due to the non-uniqueness of the imaging velocity. The problem of estimating the imaging velocity from seismic data is ill-posed. Its solution usually involves some user-dependent parameters which can be tuned to give different equally acceptable imaging velocity models (e.g. Woodward et al., 2008;Messud et al., 2017;Beraud et al., 2018). These velocity models can lead to different structural models. For example, consider the velocity modeling experiment of Beraud et al. (2018) in Figure 16 carried out using the data set in Figure 15b. Figure 16a shows the initial velocity model used for velocity modeling. Figures 16b,16c, and 16d show the final velocity model obtained by ray-based migration velocity analysis (e.g. Stork, 1992;Woodward et al., 2008) with a "low" regularization weight, a "moderate" regularization weight, and a "high" regularization weight, respectively. These 3 velocity models are considered "accurate" enough for imaging because they all flatten common-image-gathers as shown in Figures 17b,17c, and 17d. However, the 3 velocity models lead to different structural models as illustrated in Figure 18.

2. Interpretation uncertainties due to a limited resolution of seismic images. In complex geological settings where the seismic image has a limited resolution, it is usually possible to propose multiple structural interpretations of the same seismic image (e.g. Bond et al., 2007;Bond, 2015). For example, Figure 19 shows an alternative structural interpretation of the seismic image in Figure 14a by a different interpreter. In this manuscript, we will assume that we have a unique imaging velocity model so that we can only focus on interpretation uncertainties. In particular, given a set of geologically consistent structural models from a single seismic image, we will study the problem of determining which structural models are more likely than others using seismic data. For example, can we use seismic data to determine which structural model is more accurate between Figure 14b and Figure 19b? This problem will be referred to as appraising structural models/interpretations using seismic data hereafter. Our strategy for appraising structural interpretations using seismic data is illustrated in Figure 20. First, we assume that the seismic data used for seismic imaging are available. The seismic image is given to different interpreters who will propose different structural interpretations. The structural interpretations are then used to create subsurface structural models. The structural models are in turn used to create velocity models. Finally, the velocity models are used to generate synthetic data for each structural interpretation. The synthetic data are then compared with observed data in order to determine which interpretations are more likely. The main challenge of the proposed workflow is to find an appropriate way to compare observed data with synthetic data; that is, to compute appropriate data misfit values. This problem will be investigated in detail in chapters III.2 and III.3, where we show that appropriate data misfit functions should be localized in space (i.e. data misfit function as a function of space). The localization is typically achieved by projecting data residuals from the data-space into the image-domain; the projection operators are borrowed from the field of seismic imaging.

Appraising structural models using the workflow proposed above is therefore an integrative problem which involves seismic imaging and structural modeling. It is therefore reasonable to give a short review of seismic imaging and structural modeling before diving into the problem of appraising structural models using seismic data. As we were reviewing the seismic imaging and structural modeling literature, it occurred to us that we could propose alternative solutions, and sometimes simpler ones, to existing problems. Consequently, instead of giving a traditional literature review, we will present alternative solutions to some existing problems in seismic imaging and in structural modeling. As a result, this manuscript has 3 main parts:

1. Part I: Notes on Seismic Imaging.

Three optimization problems

Parts I, II, and III of the manuscript begin with an introductory text which presents seismic imaging, implicit structural modeling, and appraising structural interpretations using seismic data as optimization problems:

• The optimization approach to seismic imaging is a more modern and general (i.e. abstract) way of presenting seismic imaging.

• Presenting the implicit structural modeling problem as an optimization problem has the advantage of highlighting the challenges of interpolation for structural modeling compared to the standard scattered data interpolation problem (e.g. Anjyo et al., 2014); in particular, interpolation for structural modeling has additional optimization constraints.

• The problem of appraising structural interpretations using seismic data can also be seen from an optimization point of view; the advantage of doing so is to highlight the challenges of the problem, which appear as optimization constraints.

Contributions Seismic imaging

In chapter I.2, I propose to use reverse time migration to precondition waveform inversion. I start the chapter with a short review of similar existing preconditioning techniques in order to highlight the difference with the proposed technique.

In chapter I.3, we propose a generalized Kirchhoff extended imaging formula. While Kirchhoff extended imaging is not new, our generalized extension formula allowed us to derive a generalized differential semblance ray-based migration velocity analysis formula. The generalized extension formula also has a practical advantage: it allows a unified implementation (i.e. the same code) of multiple extended imaging operators (shot extension, offset extension, angle extension, ...) depending on input parameters.

Structural modeling

In chapter II.1, we propose an implicit structural modeling (i.e. structural interpolation) method based on finite differences. This technique has the advantage of being simple to implement and being very robust in handling discontinuities. We also derive new regularization operators in this chapter. This chapter resulted in the following publications:

-[Journal paper]: Irakarama M., Laurent G., [START_REF] Renaudeau | Implicit Structural Modeling by Minimization of the Bending Energy with Moving Least Squares Functions[END_REF] Finite difference implicit structural modeling of geological structures. In prep.

-[Conference paper]: Irakarama M., Laurent G., Renaudeau J., Caumon G. (2018) Finite difference implicit modeling of geological structures. In:. 80th EAGE Conference & Exhibition, Copenhagen, Denmark.

In chapter II.2, we propose an implicit structural modeling (i.e. structural interpolation) method based on a finite element implementation of the new regularization operators proposed in chapter II.1.

Appraising structural models using seismic data

In chapter III.2, we introduce and formalize the problem of appraising structural interpretations/models using seismic data. This chapter mainly focuses on theoretical aspects of the problem and resulted in the following publications:

-[Journal paper]: Irakarama M., Cupillard P., Caumon G., Sava P., Edwards J. (2019) Appraising structural interpretations using seismic data -theoretical elements. Geophysics, 84:N29-N40. Nominated for Geophysics Bright Spots.

-[Conference paper]: Irakarama M., Cupillard P., Caumon G, Sava P. (2017) Appraising structural models using seismic data: problem and challenges. In: 87th SEG Annual Meeting, Houston, Texas.

-[Conference paper]: Irakarama M., Cupillard P., Caumon G, Sava P. (2017) Appraising structural interpretations using seismic data misfit functionals. In: 79th EAGE Conference & Exhibition, Paris, France.

In chapter III.3, we highlight the practical challenges of appraising structural models using seismic data. We also consider appraising structural models by waveform inversion.

Other related papers

The following publications also benefited from this PhD project:

-

Part I

Notes on Seismic Imaging

Introductory Comments

Structural models are usually built from interpretations of seismic images. This first Part of the thesis is therefore devoted to seismic imaging. The problem of seismic imaging is that of transforming seismic traces, such as those illustrated in Figure 21a, into an image of the subsurface, such as that shown in Figure 21b, where geological structures can be identified. Traditionally, seismic imaging has been developed from concepts based on geophysical intuition [START_REF] Claerbout | Toward a unified theory of reflector mapping[END_REF]Schneider, 1978;Claerbout, 1985). Arguably, a more modern way to present the seismic imaging problem is to formulate it as an optimization problem (Tarantola, 1984a,b;Plessix, 2006;Fichtner et al., 2006;[START_REF] Sava | A comparative review of wavefield tomography methods[END_REF]. I follow the optimization approach to seismic imaging in this thesis as I find it to more accessible to readers, like myself, with a limited grasp of the enormous seismic imaging literature. Plessix (2006) provides a good summary for the optimization approach to seismic imaging. In this thesis, the generic (L2-norm based) seismic imaging problem is formulated as

         m(x) = argmin m(x) s,r,t D obs (r, t, s) -D syn (r, t, s, m) 2 = argmin m(x) s,r,t,x δ(x -r) D obs (x, t, s) -D syn (x, t, s, m) 2 subject to D syn (x, t, s, m) = F m(x), s , (7)
where D obs (r, t, s) denotes the observed data from a source fired at s and recorded at receiver r, D syn (r, t, s, m) denotes the corresponding synthetic data, F is the forward modeling operator, and m(x) is the subsurface Earth model. When the operator F is linear in m(x), problem 7 leads to what is known as migration and least-squares migration (Tarantola, 1984a;[START_REF] Zhang | A stable and practical implementation of least-squares reverse time migration[END_REF]; when the operator F is nonlinear in m(x), problem 7 leads to what is known as (full) waveform inversion (Tarantola, 1984b;[START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF].

Chapter I.1

Seismic Imaging

I.1.1 Introduction

Seismic imaging is a vast subject; in this chapter, I only give a summary of what I consider to be relevant to this thesis (in particular Part III). I present seismic imaging as an optimization problem (Plessix, 2006;Fichtner et al., 2006;Virieux and Operto, 2009). The reader who wishes to get more insights and a more intuition-based description of seismic imaging is referred to Claerbout (1985); Bleistein et al. (2001); Biondi (2006). Seismic imaging methods can be grouped into two main categories (Etgen et al., 2009):

1. In ray-based methods (Miller et al., 1987;[START_REF] Nemeth | Least-squares migration of incomplete reflection data[END_REF][START_REF] Thierry | Fast 2-D ray+Born migration/inversion in complex media[END_REF][START_REF] Xu | Fast migration/inversion with multivalued rayfields: Part I-Method, validation test, and application in 2D to Marmousi[END_REF], the forward modeling operator F in equation 7 is derived from ray-theory and/or the eikonal equation.

2. In wave-based methods [START_REF] Kosloff | Migration with the full acoustic wave equation[END_REF]Baysal et al., 1983;[START_REF] Mcmechan | Migration by extrapolation of time-dependent boundary values[END_REF][START_REF] Mulder | A comparison between one-way and two-way wave-equation migration[END_REF][START_REF] Zhang | Practical issues in reverse time migration: true amplitude gathers, noise removal and harmonic source encoding[END_REF][START_REF] Zhang | A stable and practical implementation of least-squares reverse time migration[END_REF], the forward modeling operator F in equation 7 is derived from the wave equation.

In this chapter, we consider seismic imaging under the acoustic approximation; that is, the forward modeling operator F in equation 7 is derived from the acoustic wave equation

L(m)P (x, t; s) = 1 m 2 (x) ∂ 2 ∂t 2 -∆ P (x, t; s) = w(t)δ(x -s), (8)
where P (x, t; s) denotes the synthetic wavefield for a source fired at s, and w(t) is the source wavelet. Equation 8 is nonlinear in the subsurface model m(x); it can be linearized by decomposing the subsurface model into a smooth varying term m 0 (x) and an oscillatory term m 1 (x) (Tarantola, 1984a): m(x) = m 0 (x) + m 1 (x). Assuming that the synthetic wavefield can also be decomposed in a similar manner: P (x, t; s) = P 0 (x, t; s) + P 1 (x, t; s), it follows from the Taylor expansion 1

m 2 (x) = 1 [m 0 (x) + m 1 (x)] 2 ≈ 1 m 2 0 (x) -2 m 1 (x) m 3 0 (x) that L(m 0 )P 0 (x, t; s) + L(m 0 )P 1 (x, t; s) = w(t)δ(x -s) + 2 m 1 (x) m 3 0 (x) ∂ 2 ∂t 2 P (x, t; s) ≈ w(t)δ(x -s) + 2 m 1 (x) m 3 0 (x)
∂ 2 ∂t 2 P 0 (x, t; s), leading to the linear (linear in m 1 ) Born modeling wave equation

   L(m 0 )P 0 (x, t; s) = w(t)δ(x -s) L(m 0 )P 1 (x, t; s) = 2 m 1 (x) m 3 0 (x) ∂ 2 ∂t 2 P 0 (x, t; s) . (9)
In equation 9, the smooth varying part of the model m 0 (x) is assumed to be known and the oscillatory part of the model m 1 (x) is the unknown subsurface model of interest; P 1 (x, t; s) represents the synthetic wavefield.

In the next two sections, I will consider the case where the forward modeling operator F in equation 7 is the solution of equation 9 and equation 8 to estimate the subsurface model from the data set in Figure 21a. This data set was generated using finite differences (Virieux, 1984) in the portion of the Marmousi model [START_REF] Versteeg | The marmousi experience: Velocity model determination on a synthetic complex data set[END_REF] shown in Figure 22a.

I.1.2 Linearized waveform inversion

Linearized seismic imaging techniques are obtained from equation 7 by equating the forward modeling operator F to the solution of the Born equation 9. In particular, the seismic imaging problem can be stated as

   m 1 (x) = argmin m 1 (x) Φ(m 1 ) = argmin m 1 (x) s,r,t,x δ(x -r) D o (x, t, s) -P 1 (x, t, s, m 1 ) 2 subject to L(m 0 )P 1 (x, t, s, m 1 ) = 2 m 1 (x) m 3 0 (x) ∂ 2 ∂t 2 P 0 (x, t; s) , (10)
where m 0 (x) and P 0 (x, t; s) are as defined in section I.1.1. In what follows, I omit parameters to simplify the notation; I also drop the subscript 1 in m 1 and P 1 . The problem 10 can be solved by iterative methods like, for example, the conjugate gradient method [START_REF] Shewchuk | An Introduction to the Conjugate Gradient Method Without the Agonizing Pain[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF]. Such iterative methods require an explicit computation of the gradient of the objective function Φ(m) with respect to the model m, which can be obtained using the adjoint-state method (Plessix, 2006;Fichtner et al., 2006). Equation 10 leads to the Lagrangian functional

L = 1, δ(x -r) D o -P 2 s,r,t,x + λ, L(m 0 )P -2 m m 3 0 ∂ 2 ∂t 2 P 0 s,t,x , (11)
where λ is the Lagrange multiplier and we have used the notation f, g x,t,.. = x,t,.. f g. The gradient of the Lagrangian with respect to m is given by

∂ ∂m L = -2δ(x -r) D o -P , ∂ ∂m P s,r,t,x + λ, L(m 0 ) ∂ ∂m P s,t,x + λ, - 2 m 3 0 ∂ 2 ∂t 2 P 0 s,t,x . (12)
Using the adjoint of the acoustic modeling operator L * , that is

λ, L(m 0 ) ∂ ∂m P s,t,x = L * (m 0 )λ, ∂ ∂m P s,t,x , (13)
equation 12 leads to

∂ ∂m L = ∂ ∂m P, L * (m 0 )λ - r 2δ(x -r) D o -P s,t,x - λ, 2 m 3 0 ∂ 2 ∂t 2 P 0 s,t,x . (14)
In practice, applying the adjoint operator L * translates to solving the acoustic wave equation 8 backwards in time (e.g. Fichtner et al., 2006). It follows from equation 14 that

   L * (m 0 )λ = 2 r δ(x -r) D o -P ∂ ∂m L = -λ, 2 m 3 0 ∂ 2 ∂t 2 P 0 s,t,x = ∂ ∂t λ, 2 m 3 0 ∂ ∂t P 0 s,t,x . ( 15
)
The derivative of the objective function Φ(m) with respect to m is given by

∂ ∂m Φ(m) = ∂ ∂m L , (16)
which essentially follows by noting from equation 11 that Φ(m) = L when P satisfies the Born equation 9. The reader is referred to Plessix (2006) and the appendix of Virieux and Operto (2009) for more details on this argument. Now that an expression of the derivative of the objective function with respect to the model is available (eq. 16), the problem in equation 10 can be solved iteratively using an algorithm of the form

             L(m 0 )P n = 2 m n m 3 0 ∂ 2 ∂t 2 P 0 L * (m 0 )λ n = 2 r δ(x -r) D o -P n ∆m n = ∂ ∂t λ n , 2 m 3 0 ∂ ∂t P 0 s,t m n+1 = m n + ∆m n , ( 17
)
Observation 1. The term ∂ ∂m L in equation 15 was computed for the entire domain; if it were computed at a fixed point x the integral over x would be absent. This is the reason why there is no integral over x in the discrete (pointwise) gradient ∆m n in equation 17.

Figure 23a shows the subsurface image obtained from the data in Figure 21a after 10 conjugate gradient iterations of the algorithm in equation 17, using the model in Figure 22b as m 0 .

Migration

Traditionally, seismic imaging seems to have been developed from concepts based on geophysical intuition [START_REF] Claerbout | Toward a unified theory of reflector mapping[END_REF]Schneider, 1978;[START_REF] Gazdag | Wave equation migration with the accurate space derivative method[END_REF]Baysal et al., 1983;Claerbout, 1985) under the name of migration. It was later shown that migration is in fact the first iteration of the algorithm in equation 17 (Tarantola, 1984a). In particular, at the first iteration of the algorithm in equation 17, we have m 0 = P 0 = 0, leading to the familiar reverse time migration formula (e.g. [START_REF] Zhang | Practical issues in reverse time migration: true amplitude gathers, noise removal and harmonic source encoding[END_REF])

   L * (m 0 )λ = 2 r δ(x -r)D o m = ∂ ∂t λ n , 2 m 3 0 ∂ ∂t P 0 s,t = -λ n , 2 m 3 0 ∂ 2 ∂t 2 P 0 s,t . (18)
An alternative approach to present migration is as the adjoint of any linearized forward seismic modeling operator (e.g. [START_REF] Claerbout | EARTH SOUNDINGS ANALYSIS: Processing versus Inversion[END_REF]. Consider for example the Born modeling formula appearing in equation 10, it is shown in the appendix (section I.1.5) that its solution can be written as

P (r, t; s) = -2 x,t ′ G(x, t -t ′ ; r) P0 (x, t ′ ; s) m(x) m 3 0 (x) , (19)
where the Green's function G(x, t; r) is the solution to

L(m 0 )G(x, t; r) = δ(x -r)δ(t).
The adjoint of the Born modeling equation 19 leads to another familiar reverse time migration formula (e.g. Symes, 2008)

m(x) = -2 r,s,t,t ′ G(x, t -t ′ ; r) P0 (x, t ′ ; s) 1 m 3 0 (x) P (r, t; s), (20)
where, in this case, P (r, t; s) is replaced by D o (r, t; s) in practice. The equivalence between equation 18 and equation 20 is shown in the appendix (section I.1.5). Figure 23b shows the subsurface image obtained by equation 20 (or alternatively equation 18).

I.1.3 Nonlinear waveform inversion

Nonlinear seismic imaging is obtained from equation 7 by equating the forward modeling operator F to the solution of the acoustic equation 8. In particular, the seismic imaging problem can be stated as

m(x) = argmin m(x) s,r,t,x δ(x -r) D o (x, t, s) -P (x, t, s, m) 2 subject to L(m)P (x, t, s, m) = w(t)δ(x -s) , (21)
In what follows, I omit parameters to simplify the notation. The procedure to solve problem 21 is very similar to that described in the previous section for solving the linearized problem 10. We first define the Lagrangian

L = 1, δ(x -r) D o -P 2 s,r,t,x + λ, L(m)P -w(t)δ(x -s)] s,t,x
; then we compute its derivative with respect the model m(x)

∂ ∂m L = -2δ(x -r) D o -P , ∂ ∂m P s,r,t,x + λ, ∂ ∂m L(m) P s,t,x + λ, L(m) ∂ ∂m P s,t,x . (22)
Noting that

∂ ∂m L(m) = - 2 m 3
∂ 2 ∂t 2 , and using the adjoint as defined in equation 13, equation 22 leads to

∂ ∂m L = ∂ ∂m P, L * (m)λ -2 r δ(x -r) D o -P s,t,x - λ, 2 m 3 ∂ 2 ∂t 2 P s,t,x . (23)
It follows from equation 23 that

   L * (m)λ = 2 r δ(x -r) D o -P ∂ ∂m L = -λ, 2 2m 3 ∂ 2 ∂t 2 P s,t,x = ∂ ∂t λ, 2 m 3 ∂ ∂t P s,t,x
, which leads to an iterative solution to problem 21 of the form

             L(m)P n = w(t)δ(x -s) L * (m)λ n = 2 r δ(x -r) D o -P n ∆m n = ∂ ∂t λ n , 2 m 3 ∂ ∂t P n s,t m n+1 = m n + ∆m n . ( 24
)
Figure 24a shows the subsurface image obtained after a multi-scale inversion [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF] using equation 24: we start with 10 conjugate gradient iterations up to 6 Hz, then 10 iterations up to 8 Hz, 10 Hz, and finally up to 12 Hz. In practice, both linearized and nonlinear inversion require an efficient preconditioner. For example, Figure 24b shows the subsurface image that would have been obtained if no preconditioner was used in Figure 24a. Preconditioning is discussed in chapter I.2.

I.1.4 Discussion

In section I.1.2, we presented a nonlinear conjugate gradient solution to linearized inversion; this was mainly to highlight the similarity with nonlinear inversion. Linearized inversion can also be performed using linear conjugate gradients. The requirements for solving a problem with linear conjugate gradients are:

1. The forward modeling operator should be linear with respect to the model m(x), which is the case for the Born modeling equation 19.

2. The adjoint of the forward modeling operator should be available, which is equation 20 in our example.

In practice, the linearized inversion scheme in equation 17 where all shot gathers are inverted simultaneously did not work very well. The image in Figure 23a was obtained by inverting each shot gather independently and then stacking the resulting partial images. Figure 25 compares images obtained by shotwise inversion against images obtained by simultaneous inversion of all the shot gathers.

The discussion in section I.1.2 suggests that a seismic image has the same units as velocity (m 1 in equation 9). It is customary however to plot seismic images without units. This is because seismic images usually undergo a series of post-processing steps for better visual effects. Figure 26 shows an application of some common post-processing steps to the reverse-time migration image in Figure 23b. In chapter I.3, it is shown that replacing P 0 (x, t ′ ; s) with G(x, t ′ ; s) in equations 19 and 20 (which corresponds to replacing the source wavelet w(t) in equation 9 by δ(t)), and then replacing the Green's functions G(x, tt ′ ; r) and G(x, t ′ ; s) by their ray-based counterpart, equations 19 and 20 become Kirchhoff modeling and Kirchhoff migration operators, respectively. Figure 27 shows the subsurface images obtained by ray-based least-squares migration and ray-based migration, after post-processing as described in Figure 26. In this paragraph, we demonstrate the Born modeling equation 19. We start from the Born modeling formula in equation 10, that is

L(m 0 )P (x, t; s) = 2 m(x) m 3 0 (x) ∂ 2 ∂t 2 P 0 (x, t; s). (25)
For simplicity, we define

g(x, t; s) = 2 m(x) m 3 0 (x) P0 (x, t; s),
allowing us to express equation 25 as

L(m 0 )P (x, t; s) = g(x, t; s) (26) = x ′ ,t ′ δ(x -x ′ )δ(t -t ′ )g(x ′ , t ′ ; s) = x ′ ,t ′ L(m 0 )G(x, t; x ′ , t ′ )g(x ′ , t ′ ; s). (27)
In equation 27, we have defined

L(m 0 )G(x, t; x ′ , t ′ ) = δ(x -x ′ )δ(t -t ′ ). (28)
It follows from equation 27 that L(m 0 ) P (x, t; s) -

x ′ ,t ′ G(x, t; x ′ , t ′ )g(x ′ , t ′ ; s) = 0, leading to P (x, t; s) = x ′ ,t ′ G(x, t; x ′ , t ′ )g(x ′ , t ′ ; s). ( 29
)
Equation 29 represents the solution to equation 25 in terms of the Green's function G(x, t; x ′ , t ′ ) (i.e., the solution to equation 28). The modeled seismic traces are obtained by sampling P (x, t; s) at receivers r; that is,

P (r, t; s) = x δ(x -r)P (x, t; s) = x δ(x -r) x ′ ,t ′ G(x, t; x ′ , t ′ )g(x ′ , t ′ ; s) = x ′ ,t ′ G(r, t; x ′ , t ′ )g(x ′ , t ′ ; s) = x ′ ,t ′ G(r, t -t ′ ; x ′ )g(x ′ , t ′ ; s) (30) = x ′ ,t ′ G(x ′ , t -t ′ ; r)g(x ′ , t ′ ; s) (31) = x,t ′ 2 m(x) m 3 0 (x) G(x, t -t ′ ; r) ∂ 2 ∂t ′2 P 0 (x, t ′ ; s) = - x,t ′ 2 m(x) m 3 0 (x) ∂ ∂t ′ G(x, t -t ′ ; r) ∂ ∂t ′ P 0 (x, t ′ ; s) (32) = - x,t ′ 2 m(x) m 3 0 (x) ∂ ∂t ′ τ ∂ ∂τ G(x, τ ; r) ∂ ∂t ′ P 0 (x, t ′ ; s), with τ = t -t ′ (33) = x,t ′ 2 m(x) m 3 0 (x) ∂ ∂τ G(x, τ ; r) ∂ ∂t ′ P 0 (x, t ′ ; s), with τ = t -t ′ (34) = 2 x,t ′ Ġ(x, t -t ′ ; r) Ṗ0 (x, t ′ ; s) m(x) m 3 0 (x) , (35)
which demonstrates equation 19.

In equation 30, we have used

G(x, t; x ′ , t ′ ) = G(x, t -t ′ ; x ′ , 0) = G(x, t -t ′ ; x ′ ). ( 36
)
This follows from the definition of G(x, t; x ′ , t ′ ) in equation 28; in particular,

δ(x -x ′ )δ(t) = δ(x -x ′ )δ(t -0) = L(m 0 )G(x, t; x ′ , 0) = δ(x -x ′ )δ(t + t ′ -t ′ ) = L(m 0 )G(x, t; x ′ , 0) using l = t + t ′ → t = l -t ′ , it follows that = δ(x -x ′ )δ(l -t ′ ) = L(m 0 )G(x, l -t ′ ; x ′ , 0) L(m 0 )G(x, l; x ′ , t ′ ) = δ(x -x ′ )δ(l -t ′ ) = L(m 0 )G(x, l -t ′ ; x ′ , 0)
In equation 31, we have used

G(x, t; x ′ , t ′ ) = G(x ′ , t; x, t ′ ). ( 37
)
This is follows from the definition of G(x, t; x ′ , t ′ ) in equation 28; in particular,

L(m 0 )G(x, t; x ′ , t ′ ) = δ(x -x ′ )δ(t -t ′ ) = δ(x ′ -x)δ(t -t ′ ) = L(m 0 )G(x ′ , t; x, t ′ )
The equivalence between equation 18 and equation 20

In this paragraph, we demonstrate the equivalence between equation 18 and equation 20. First, let us define g(r, x, t, s) = 2 r δ(xr)D o (x, t, s), so that equation 18 can simplify to

   L * (m 0 )λ = g(r, x, t, s) m = -λ, 2 m 3 0 ∂ 2 ∂t 2 P 0 s,t . (38)
By looking at equation 26 and its solution, equation 29, it follows that the solution to the first equation in 38 is

λ = x ′ ,t ′ G(x, t; x ′ , t ′ )g(r, x ′ , t ′ ; s) = 2 x ′ ,t ′ G(x, t; x ′ , t ′ ) r δ(x ′ -r)D o (x ′ , t ′ , s) = 2 x ′ ,t ′ ,r G(x, t; x ′ , t ′ )δ(x ′ -r)D o (x ′ , t ′ , s) = 2 t ′ ,r G(x, t; r, t ′ )D o (r, t ′ , s) = 2 t ′ ,r G(x, t ′ ; r, t)D o (r, t ′ , s) (39) = 2 t ′ ,r G(x, t ′ -t; r)D o (r, t ′ , s), (40)
where we have used equation 37 in equation 39, and equation 36 in equation 40. Inserting the expression of λ, equation 40, in the second equation in 38 leads to

m(x) = -2 s,t 2 m 3 0 (x) ∂ 2 ∂t 2 P 0 (x, t; s) r,t ′ G(x, t ′ -t; r)D o (r, t ′ , s) = -2 r,s,t,t ′ G(x, t ′ -t; r) ∂ 2 ∂t 2 P 0 (x, t; s) 2 m 3 0 (x) D o (r, t ′ , s),
which demonstrates the equivalence between equation 18 and equation 20.

Chapter I.2

Reverse time migration as a preconditioner for waveform inversion

I.2.

Summary

Waveform inversion has gained recognition as a powerful technique for imaging the subsurface quantitatively by fitting synthetic data to observed data. On the one hand, linearized waveform inversion estimates the subsurface reflectivity model and, on the other hand, nonlinear waveform inversion estimates the velocity model in the subsurface. Waveform inversion problems are usually solved using iterative methods due to their large sizes. The high computational cost of waveform inversion is a challenge that has motivated geoscientists to develop preconditioning techniques to accelerate these iterative methods. In this chapter, I propose to use the reverse time migration (RTM) image as a preconditioner. The premise of the proposed method is that because the RTM image is obtained by a zero-lag cross-correlation of the source and receiver wavefields, the RTM image is a measure of both the subsurface reflectivity and illumination for a given acquisition geometry. I argue that the low frequency background signal traditionally regarded as noise in conventional RTM imaging is an approximation of the diagonal elements of the Hessian matrix of the objective function and can therefore be exploited to accelerate waveform inversion.

I.2.2 Introduction

Linearized waveform inversion (Tarantola, 1984a;[START_REF] Zhang | A stable and practical implementation of least-squares reverse time migration[END_REF][START_REF] Rocha | 3D acoustic least-squares reverse time migration using the energy norm[END_REF] and nonlinear waveform inversion (Tarantola, 1984b;[START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF] allow to generate optimal subsurface images by fitting synthetic data to observed data. Linearized inversion estimates the subsurface reflectivity, while nonlinear inversion seeks to estimate the subsurface velocity. Both linear and nonlinear waveform inversion problems are usually solved using iterative methods. Popular solvers for waveform inversion problems include gradient methods and quasi-Newton methods (Virieux and Operto, 2009). Gradient solvers generally show a low convergence rate, requiring preconditioning in order to make them competitive with quasi-Newton methods. The objective of this chapter is to propose a simple and efficient preconditioner for gradient methods.

A popular technique for preconditioning gradient based waveform inversion is to design a preconditioner that approximates the diagonal of the Hessian matrix of the L2-norm objective function (e.g. [START_REF] Plessix | Frequency-domain finite-difference amplitude-preserving migration[END_REF])

H(x, y) = s,r,t τ Ġr (t -τ, x) Ġs (τ, x) τ Ġr (t -τ, y) Ġs (τ, y) ,
where Ġs and Ġr represent time derivatives of the source and receiver wavefields; x and y are two points in space. [START_REF] Shin | Improved amplitude preservation for prestack depth migration by inverse scattering theory[END_REF] proposed to make the assumption that the receiver wavefields' contribution to the Hessian is negligible, resulting into a diagonal of the Hessian matrix obtained by a product of source wavefields; [START_REF] Guitton | Smoothing imaging condition for shot-profile migration[END_REF] further proposed to smooth the product of the source wavefields for more stability. [START_REF] Plessix | Frequency-domain finite-difference amplitude-preserving migration[END_REF] show that the product of source wavefields approximation is valid under infinite and uniform receiver coverage; they therefore proposed to approximate the receiver wavefields by analytical formulae, derived for a constant velocity model, for a finite receiver coverage acquisition. [START_REF] Rickett | Illuminationbased normalization for waveequation depth migration[END_REF] and [START_REF] Guitton | Amplitude and kinematic corrections of migrated images for nonunitary imaging operators[END_REF] use a first migrated image to generate synthetic data, which are migrated to produce a second image; they note that the two images are related by an approximation of the Hessian matrix, thus allowing them to estimate the Hessian. [START_REF] Hou | Accelerating extended least-squares migration with weighted conjugate gradient iteration[END_REF] use ray theory approximation to derive a pseudo inverse of the (extended) Born modeling operator; they show that the inverse of the Born modeling operator can be approximated by a weighted version of RTM. I reach a similar conclusion here. In particular, I argue that, under some crude assumptions, the diagonal elements of the pseudo Hessian matrix can be approximated by a low pass filter of the RTM image. This implies that the low frequency background signal in conventional RTM imaging, traditionally regarded as noise [START_REF] Zhang | Practical issues in reverse time migration: true amplitude gathers, noise removal and harmonic source encoding[END_REF], can be exploited to accelerate waveform inversion.

In section I.2.3.1, I present a heuristic motivation for the proposed method based on the fact that an RTM image is a good measure of the illumination distribution in the subsurface for a given acquisition geometry. In section I.2.3.2, I attempt to provide a more formal justification for the proposed method; the objective of this part is to reveal the underlying assumptions of the proposed method. In section I.2.4, I challenge the proposed method with numerical experiments on the Marmousi model [START_REF] Versteeg | The marmousi experience: Velocity model determination on a synthetic complex data set[END_REF]; both linearized inversion and nonlinear inversion are tested.

I.2.3 Theory I.2.3.1 Heuristics

The RTM image, under the Born approximation, can be expressed as (e.g. Symes, 2008, see also the last equation in section I.1.5)

I(x) = 1 m 0 (x) 3 s,r,t,τ Ġ(t -τ, x, r) Ġ(τ, x, s)D obs (s, r, t) = s,r,t B(x, s, r, t)D obs (s, r, t) = B * D obs (s, r, t) , (41)
where

B[•] : I(x) → D obs (s, r, t)
is the Born modeling operator (see for example section I.1.1), the square brackets indicate a linear relationship and the asterisk symbol indicates the adjoint operation; Ġ(t, x, s) denotes the time derivative of the Green's function of the acoustic wave equation for a source positioned at s, and D obs (s, r, t) are the input data recorded at r; m 0 (x) is the smooth background velocity model, and I(x) is the output RTM image. According to equation 41, if we ignore the data term D obs (s, r, t) for the sake of argument, the RTM image is a measure of the illumination distribution in the subsurface as a function of acquisition geometry and the background velocity model. The effect of the data term can be thought of as multiplying the illumination map by a highly oscillatory signal, whose effect can be removed by smoothing the absolute value of the output RTM image.

As an example, consider the synthetic model in Figure 28.a, discretized on a regular grid (dx=20m, dz=10m). For a survey with 39 sources, 200 receivers evenly distributed at the surface, and a Ricker source wavelet with a dominant frequency of 20Hz, the raw RTM image (using non muted data) is shown in Figure 28.d. Smoothing the absolute value of RTM image removes its oscillatory component, resulting in the illumination map in Figure 28.e; this illumination map can be used to normalize the raw RTM image (Figure 28.f). Numerical experiments suggest that this illumination map can be used to accelerate gradient methods for waveform inversion. Figure 29 compares the convergence rate for the linearized inversion (least squares RTM) problem using a conjugate gradient solver (e.g. [START_REF] Shewchuk | An Introduction to the Conjugate Gradient Method Without the Agonizing Pain[END_REF] with and without the proposed preconditioner; the preconditioned version shows a convergence rate about an order of magnitude higher. Figure 30 compares the convergence rate for the nonlinear waveform inversion (full waveform inversion) problem using a nonlinear conjugate gradient solver with and without the proposed preconditioner; the preconditioned version shows a convergence rate about two orders of magnitude higher.

I.2.3.2 Underlying assumptions

I now attempt to formally justify the proposition that, under suitable assumptions, the smoothed RTM image approximates the diagonal of the Hessian matrix. Assuming that the recorded data D obs (s, r, t) can be approximated by Born modeling, for some unknown reference reflectivity image I ref (x), it follows from equation 41 that [START_REF] Rickett | Illuminationbased normalization for waveequation depth migration[END_REF][START_REF] Guitton | Amplitude and kinematic corrections of migrated images for nonunitary imaging operators[END_REF])

I(x) = s,r,t B(x, s, r, t)D obs (s, r, t) = s,r,t B(x, s, r, t) y B(y, s, r, t)I ref (y) = y H(x, y)I ref (y) = H I ref (y) , ( 42
)
where H is the approximate Hessian matrix. First, let us assume that the Hessian H can be approximated by its diagonal elements, denoted by h(x); this reduces equation 42 to

I(x) = h(x)I ref (x). ( 43
)
According to equation 43, the RTM image is the product of the diagonal elements of the Hessian and the reference reflectivity. The idea will be to smooth out the oscillatory reference reflectivity. Let N (x) denote the set of points "near" x, and let f (x) denote a low-pass filter.

Define the filtered version of equation 43 as

Ī(x) = x i ∈N (x) f (x i )I(x i ) = x i ∈N (x) f (x i )h(x i )I ref (x i ). ( 44
)
If the diagonal elements h(x) vary much slower than the (unknown) reference reflectivity I ref (x), h(x) can be approximated by a quasi-constant function in N (x), leading to I(x) ≈ ĥ(x)

x i ∈N (x) f (x i )I ref (x i ), (45)
where ĥ(x) is the average value of h(x) in the neighborhood N (x). Finally, let us suppose that the reference reflectivity I ref (x) oscillates around some (unknown) average value, denoted by i ref , and that the smoothed reference reflectivity Īref (x) does not deviate too much from i ref ; it follows from 45 that

I(x) ≈ ĥ(x) Īref (x), ≈ i ref ĥ(x). ( 46
)
The average value of I ref (x), i ref , is likely to be zero in some cases; this is circumvented by first taking the absolute value of I(x) before filtering. The unknown i ref has no effect on image normalization since it is a constant. In summary, the smoothed RTM image approximates the Hessian under the following assumptions:

1. The Hessian H can be approximated by its diagonal elements, denoted by h(x).

2. The diagonal elements h(x) vary much slower than the (unknown) reference reflectivity

I ref (x).
3. The reference reflectivity I ref (x) oscillates around some (unknown) average value, denoted by i ref ; furthermore, the smoothed reference reflectivity does not deviate too much from i ref .

The first assumption is valid under a high frequency approximation (e.g. Miller et al., 1987); this may be a crude approximation in some cases but it has proven useful in practice [START_REF] Plessix | Frequency-domain finite-difference amplitude-preserving migration[END_REF]. The third assumption is clearly reasonable under the Born approximation.

I.2.4 Numerical examples

I now test the proposed preconditioner on the portion of the Marmousi velocity model [START_REF] Versteeg | The marmousi experience: Velocity model determination on a synthetic complex data set[END_REF] shown in Figure 31a. The model is discretized on a grid of size 200 × 215 with dx = 30m and dz = 15m. The model is used to generate synthetic data by finite difference time domain acoustic wave simulation (Virieux, 1984) with perfectly matched absorbing boundaries (Collino and Tsogka, 2001); the source is a Ricker wavelet with a dominant frequency of 15Hz. The input velocity model is also used to simulate the imaging velocity model shown in Figure 31b by smoothing. Figure 32 compares the results obtained by a multi-scale nonlinear waveform inversion [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF] without preconditioning (Figure 32a) and with preconditioning (Figure 32b). The preconditioned inversion in Figure 32b clearly converges much faster than the non-preconditioned version in Figure 32a. The difference between the results obtained for linearized inversion without preconditioning (Figure 33a) and with preconditioning (Figure 33b) is less striking than that previously observed for nonlinear inversion (Figure 32). However, a look at the residual curves in Figure 33c reveals that the preconditioned linearized inversion converges about an order of magnitude faster than the non-preconditioned version.

I.2.5 Conclusion

Waveform inversion has proven to be a promising technology for quantitative seismic depth imaging. One of the challenges of waveform inversion is its high computational cost. A number of techniques have been proposed to accelerate the convergence rate of waveform inversion problems. In this paper, I propose another method to accelerate gradient based waveform inversion by using the low frequency background signal in reverse time migration (RTM) images as a preconditioner. Numerical experiments suggest that the proposed preconditioner accelarates notably both linearized and nonlinear waveform inversion. The success of the preconditioner is justified by showing that, under appropriate assumptions, the low frequency background RTM signal can approximate the diagonal elements of the Hessian matrix of the objective function.

I.2.6 Acknowledgements

A lot of what I know about seismic imaging comes from discussions with Daniel Rocha, Esteban Diaz, Ludovic Métivier, Paul Cupillard, and Paul Sava. My knowledge of seismic imaging has also greatly benefited from exchanges within the HIWAI project led by Yann Capdeville. I am thankful to Guillaume Caumon for having sent me to many EAGE and SEG seismic imaging conferences and workshops over the course of my PhD. I would also like to thank Melchior Schuh-senlis for many discussions about the algebraic manipulations in I.2.3.2.

Chapter I.3

Extended Kirchhoff Imaging for Migration Velocity Analysis

This chapter resulted from an ongoing technical report co-authored with Paul Sava and Thomas Beraud.

I.3.1 Summary

Extended images are a substantial component of velocity analysis. We investigate extensions of the Kirchhoff imaging operator that capture image distortions caused by inaccurate velocity. We propose a general formula that describes all data-domain extensions (e.g. offset and reflection angle extensions), as well as image-domain extensions (e.g. space lag and time lag extensions). Both single-arrival and multi-arrival imaging are discussed. We propose a method to handle multi-arrivals without shooting rays from the subsurface, thus reducing the computational cost. One practical advantage of our generalized extension formula is that it allows a unified implementation for all the extensions mentioned herein. Another advantage of the proposed generalized extension formula is that it leads to a unified differential semblance migration velocity analysis formulation for all the data-domain extensions mentioned herein.

I.3.2 Introduction

Kirchhoff extended images, also known as common image gathers, find applications in velocity analysis (Stork, 1992), imaging in complex media (Xu et al., 2001), and reservoir characterization [START_REF] Fatti | Detection of gas in sandstone reservoirs using AVO analysis: A 3-D seismic case history using the Geostack technique[END_REF]. In this chapter, extended images are grouped into two classes: data-domain extensions, and image-domain extensions. Data-domain extensions are those for which the extension parameter depends explicitly on data parameters, such as source-receiver geometry; these include shot extension, receiver extension, offset extension and reflection angle extension, and have been extensively used in ray based imaging and velocity analysis (Xu et al., 2001;Woodward et al., 2008). Image-domain extensions are those for which the extension parameters do not depend explicitly on data parameters; these include time-lag and space-lag extensions extensively used in wave-equation based imaging and velocity analysis [START_REF] Sava | Time-shift imaging condition in seismic migration[END_REF]Symes, 2008;[START_REF] Sava | Extended imaging conditions for wave-equation migration[END_REF].

In this chapter, we propose a general Kirchhoff extension formula that describes all the extensions mentioned above. The generalized extension formula naturally leads to a unified implementation for constructing Kirchhoff extended images, and to a unified differential semblance [START_REF] Symes | Velocity inversion by differential semblance optimization[END_REF][START_REF] Mulder | Automatic velocity analysis by differential semblance optimization[END_REF] migration velocity analysis (MVA) formulation. First, we propose the general Kirchhoff extension formula for singlearrival imaging, then we propose a strategy to handle multi-arrivals without shooting rays from subsurface. We then apply the proposed extension formula within a typical ray-based MVA workflow (Stork, 1992;Woodward et al., 2008) to measure its suitability for velocity analysis.

Finally, the proposed extension formula is used to formulate a differential semblance MVA that is simultaneously valid for multiple data-domain extensions.

I.3.3 Imaging with single arrivals I.3.3.1 Data-domain extensions

The standard Kirchhoff imaging integral (e.g. Etgen et al., 2009) can be expressed as

I(x) = s,r,t W (s, r, x)δ t -T (s, r, x) ∂ ∂t D(s, r, t) = s,r W (s, r, x) ∂ ∂t D s, r, T (s, r, x) = s,r F 1 (s, r, x). (47)
In equation 47, x is an image point in the subsurface; s and r are source and receiver coordinates, assumed to be at the free surface; T (s, r, x) is the traveltime of a ray shot from s, scattered off x, then recorded at r; D(s, r, t) are the input data; W (s, r, x) are integration weights, whose physical meaning will become apparent when we discuss about image-domain extensions further below; I(x) is the output image. In this chapter, we define data-domain extensions as those for which the extension parameter depends explicitly on data parameters (s, r, x). Given this definition, all the data-domain extensions of equation 47 can be constructed using the generalized Radon transform

I(x, β) = s,r δ β -ζ(s, r, x) F 1 (s, r, x), (48)
where I(x, β) is the output extended image, and β is the extension parameter whose physical meaning depends on the user defined function ζ(s, r, x). β can be a scalar or a vector depending on the definition of ζ. The advantage of equation 48 is that it allows multiple extensions to be constructed using the same kernel of code. For example, the shot extended image is obtained by defining ζ(s, r, x) = s;

(49) the surface offset extended image is obtained by defining

ζ(s, r, x) = |s -r| 2 ; ( 50
)
and the reflection angle extended image is obtained by defining

ζ(s, r, x) = acos ∇T (x, s) • ∇T (x, r) ∇T (x, s) • ∇T (x, s) , (51)
where T (x, s) denotes first arrival traveltime maps from a source positioned at s.

I.3.3.2 Image-domain extensions

The standard wavefield imaging integral (e.g. Symes, 2008) can be expressed as

I(x) = s,r,t,t ′ G(t ′ , x, s)G(t -t ′ , x, r) ∂ ∂t D(s, r, t), (52)
where G(t, x, s) denotes the Green's function of the acoustic wave equation with a source positioned at s, and D(s, r, t) are the input data. Following [START_REF] Sava | Extended imaging conditions for wave-equation migration[END_REF], we introduce the time-lag (scalar) parameter τ , and space-lag (vector) parameter λ, and define an extension of equation 52 as

I(x, λ, τ ) = s,r,t,t ′ G (t ′ -τ ), x -λ, s × G t -(t ′ + τ ), x + λ, r ∂ ∂t D(s, r, t) = s,r,t,t ′ G(t ′ -τ, x -, s) × G(t -t ′ -τ, x + , r) ∂ ∂t D(s, r, t). (53)
Using the ray approximation

G(t, x, s) = A(x, s)δ t -T (x, s) , (54)
where A(x, s) and T (x, s) represent the amplitude and traveltime of a ray from a source at s, equation 53 leads to Let x = (x, y, z) denote a point in 3D; the traveltime map T (x, s), from a shot positioned at s, is a cube that can be multi-valued when the velocity field has strong spatial variations.

I(x, λ, τ ) = s,r,t,t ′ A(x -, s)δ t ′ -τ -T (x -, s) × A(x + , r)δ t -t ′ -τ -T (x + , r) ∂ ∂t D(s, r, t) = s,r,t A(x -, s)A(x + , r) × δ t -2τ -T (x -, s) -T (x + , r) ∂ ∂t D(s, r, t) = s,r W (s, r, x, λ) ∂ ∂t D s, r, 2τ + T (s, r, x, λ) = s,r F 2 (s, r, x, λ, τ ). ( 55
Let us denote the two ray take-off angles (dip and strike) by φ = (φ 1 , φ 2 ), and introduce the 5D coordinates x * = (x, φ). x * will be referred to as phase-space coordinates (e.g. [START_REF] Cerveny | Seismic ray theory[END_REF][START_REF] Bashkardin | Phase-Space Imaging of Reflection Seismic Data[END_REF], and x as physical-space coordinates. The phase-space traveltime map T (x * , s) is always single-valued; this concept is easier to illustrate in 2D as shown in Figure 34. Instead of projecting the single-valued phase-space traveltime map T (x * , s) to the lower dimensional multi-valued physical-space traveltime map T (x, s), which is later re-interpolated during the imaging loop (e.g. [START_REF] Xu | Fast migration/inversion with multivalued rayfields: Part I-Method, validation test, and application in 2D to Marmousi[END_REF], we propose to interpolate ray information from T (x * , s) directly during the imaging loop. That is, we propose the multi-valued version of equation 56 as

I(x, β, λ, τ ) = s,r,φ δ β -ζ(s, r, x, λ, φ) W (s, r, x, λ, φ) × ∂ ∂t D s, r, 2τ + T (s, r, x, λ, φ) . ( 57
)
In practice, the traveltime map T (x * , s) is discretized into a set of cells {C i (x * , s)}. We loop over the cells in {C i (x * , s)} implicitly through grid points x. We are allowed to do this because the order in which the cells are visited is irrelevant. A pseudo-code to implement equation 57 is given in algorithm 1. for each source s do for each receiver r do for each image-point x do x -= xλ;

x + = x + λ; for C s ∈ {C(x * j , s) | x * j = (x -, φ j ) } do for C r ∈ {C(x * k , r) | x * k = (x + , φ k ) } do β = ζ(C s , C r ); t = T (C s ) + T (C r ) + 2τ ; I(x, β, λ, τ )+ = A(C s )A(C r ) ∂ ∂t D(s, r, t); end end end end end Algorithm 1: A general pseudo-code to construct Kirchhoff extensions. {C(x * i , s) | x * = (x, φ i )}
is the set of cells that belong to a phase-space traveltime map T (x * , s) from a shot positioned at s and whose projection in the physical-space contains the point x. The notation f (C) means that the value of f (x, s) is estimated by interpolation from the cell C(x, φ, s). A 2D example comparing single-arrival with multi-arrival imaging is shown in Figure 35.

I.3.5 Velocity model building via ray tracing

We now test the ability of the proposed extended imaging operators (equations 56 and 57) to capture image distortion caused by inaccurate velocity by an application to ray-based MVA (Stork, 1992;Woodward et al., 2008). We use the simple velocity model in Figure 36a as the true model, in which we generate "observed" data by finite-difference acoustic modeling (Virieux, 1984) with perfectly matched absorbing boundaries (Collino and Tsogka, 2001). Sources and receivers cover the free surface. We then attempt to build a migration velocity model starting from the initial velocity model in Figure 36b. Figure 37a shows the output of equation 56 using the initial model (Figure 36b), with τ = 0, λ = 0, and ζ as defined in equation 50. For this simple example, residual moveout curves can be picked automatically by parabola fitting as shown in Figure 37b. Rays are then traced from each residual moveout curve towards the surface (Figure 37c), making it possible to build a ray tomographic system which can then be inverted to update the velocity model. Figure 37d shows the velocity model obtained after 15 iterations; the tomographic system is regularized using the smoothing operators proposed in chapter II.1 and in Irakarama et al. (2018). (Stork, 1992;Woodward et al., 2008). a) Offset extended image using the initial velocity model (Figure 36b). 36b) with common image gathers computed in the updated model (Figure 37d). We clearly note that the initial gathers (Figure 38a) are flattened (Figure 38b) by the updated velocity model. Figure 39 compares the migrated image obtained using the initial velocity model (Figure 39a) with the migrated image obtained using the updated velocity model (Figure 39b). The migrated image resulting from the updated model is more focused and has less artifacts. The analysis done in this section confirms that the proposed extended imaging operators can indeed be used for velocity modeling.

Figure 38: a) Common image gathers obtained using the initial velocity model (Figure 36b). b) Common image gathers obtained using the updated velocity model (Figure 37d).

Figure 39: a) Migrated image obtained using the initial velocity model (Figure 36b). b) Migrated image obtained using the updated velocity model (Figure 37d).

I.3.6 Velocity model building via the adjoint-state method

In this section, we investigate using the proposed generalized extended imaging operators to develop a general gradient for gradient-based MVA. Velocity modeling can also be formulated as an optimization problem (e.g. [START_REF] Mulder | Automatic velocity analysis by differential semblance optimization[END_REF]; we consider such a formulation in this section using a differential semblance objective function similar to that of Plessix (2006). For simplicity, we only investigate velocity modeling using the single arrival data-domain extended imaging equation 48. The velocity modeling problem in this framework can be stated as

                     m(x) = argmin m(x) Φ(m) = argmin m(x) x,β I(x, β, m) -I(x, β + 1, m) 2 subject to I(x, β, m) = s,r δ β -ζ(s, r, x, m) Ḋ s, r, T (s, r, x) T (s, r, x, m) = T s (x, m) + T r (x, m) ∇T s (x, m) 2 = 1 m(x) 2 , T s (s, m) ∇T r (x, m) 2 = 1 m(x) 2 , T r (r, m) , ( 58
)
where the notation I(x, β, m)-I(x, β+1, m) denotes a difference of two neihgbouring (partial) images along the β-axis. The problem in equation 58 can be solved using the adjoint-state method (e.g. Plessix, 2006). Other authors that have used the adjoint-state method and an eikonal solver for velocity modeling include Taillandier et al. ( 2009) who used the adjoint-state method to solve a first-arrival traveltime tomography problem, and Tavakoli et al. ( 2017) who used the adjoint-state method to solve a slope tomography problem. Here we consider using the adjoint-state method to solve a differential semblance ray-based MVA problem as done in [START_REF] Mulder | Automatic velocity analysis by differential semblance optimization[END_REF]. Our approach differs from that of [START_REF] Mulder | Automatic velocity analysis by differential semblance optimization[END_REF] in two ways. First, we use an eikonal solver (Sethian and Popovici, 1999) to compute traveltimes instead of a ray-tracer [START_REF] Vinje | Traveltime and amplitude estimation using wavefront construction[END_REF]. Second, because we use the generalized data-domain extended imaging equation 48, our result is simultaneously valid for multiple datadomain extensions. In what follows, we omit some arguments in expressions in equation 58 to simplify notation. The Lagrangian of the problem in equation 58 is given by

L = 1, I(β) -I(β + 1) 2 x,β + λ 1 , I(β) - s,r δ β -ζ(s, r) Ḋ(s, r) x,β + λ 2 , ∇T s 2 - 1 m 2 x,s + λ 3 , ∇T r 2 - 1 m 2 x,r . ( 59
)
The derivative of the Lagrangian with respect to the model is given by (see the appendix, section I.3.10, for more details)

∂ ∂ m L = 2 -I(β + 1) + 2I(β) -I(β -1) , ∂ ∂ m I(β) x,β + λ 1 , ∂ ∂ m I(β) x,β -λ 1 , s,r δ β -ζ D ∂ ∂ m T s x,β -λ 1 , s,r δ β -ζ D ∂ ∂ m T r x,β + λ 2 , 2∇T s • ∇ ∂ ∂ m T s x,s + λ 2 , 2 m 3 x,s + λ 3 , 2∇T r • ∇ ∂ ∂ m T r x,r + λ 3 , 2 m 3 x,r . (60)
The divergence theorem implies that

F (x), ∇g(x) x = ∇ • F (x), g(x) x , (61)
for some vector valued function F and a scalar valued function g (we have assumed that g or F is zero at boundaries). Using equation 61, equation 60 can be written as

∂ ∂ m L = λ 1 + 2 -I(β + 1) + 2I(β) -I(β -1) , ∂ ∂ m I(β) x,β + 2∇ • λ 2 ∇T s - r,β λ 1 δ β -ζ D, ∂ ∂ m T s x,s + λ 2 , 2 m 3 x,s + 2∇ • λ 3 ∇T r - s,β λ 1 δ β -ζ D, ∂ ∂ m T r x,r + λ 3 , 2 m 3 x,r , leading to                λ 1 = 2 I(β + 1) -2I(β) + I(β -1) 2∇ • λ 2 ∇T s = r,β λ 1 δ β -ζ D 2∇ • λ 3 ∇T r = s,β λ 1 δ β -ζ D ∂ ∂m L = λ 2 , 2 m 3 x,s + λ 3 , 2 m 3 x,r . ( 62
)
The gradient in equation 62 allows to solve the problem in equation 58 using an iterative algorithm of the form

                                         ∇T n s 2 = 1 m 2 n (a) ∇T n r 2 = 1 m 2 n (b) I n (β) = s,r δ β -ζ Ḋ(s, r, T n s + T n r ) (c) λ n 1 = 2 I n (β + 1) -2I n (β) + I n (β -1) (d) 2∇ • λ n 2 ∇T n s = r,β λ n 1 δ β -ζ D (e) 2∇ • λ n 3 ∇T n r = s,β λ n 1 δ β -ζ D (f ) ∆m n = λ n 2 , 2 m 3 n s + λ n 3 , 2 m 3 n r (g) m n+1 = m n + ∆m n (h) . ( 63

I.3.7 Discussion

The multivalued imaging equation 57 is included merely for completeness. In practice, the single valued version of the same equation is much easier to implement and gives comparable results. We found the multivalued imaging formula to perform considerably better only for idealized examples such as that in Figure 35. For more realistic models, like the Marmousi [START_REF] Versteeg | The marmousi experience: Velocity model determination on a synthetic complex data set[END_REF], for example, the multivalued formula did not perform considerably better than its single valued version. The challenge in implementing the multivalued imaging formula 57 comes from implementing efficiently the expression

{C(x * i , s) | x * = (x, φ i )}
in algorithm 1; that is, given a point x, find the set of ray-cells that belong to a phase-space traveltime map T (x * , s) from a shot positioned at s, and such that the projection of each raycell in the physical-space contains the point x.

The discussion in section I.3.6, particularly the algorithm in equation 63, has not yet been tested1 . Therefore, at the time of writing this manuscript, we are unable to comment on how well the algorithm in equation 63 would work in practice. There are, however, some familiar expressions in equation 63: the expression in equation 63d is similar to that obtained by Plessix (2006) for shot extended wave-based differential semblance MVA, and the expression in equation 63g is similar to that obtained by [START_REF] Tavakoli | Slope tomography based on eikonal solvers and the adjoint-state method[END_REF]. What is new in equation 63 are, to our knowledge, the expressions on the right-hand side of equations 63c, 63e, and 63f.

I.3.8 Conclusion

Image extension is a key component of modern velocity analysis and imaging. In wave-based velocity analysis, extended imaging is mainly dominated by image-domain extensions such as time-lag and space-lag extensions. In ray-based velocity analysis, extended imaging is mainly dominated by data-domain extensions such as offset and reflection angle extensions. In this paper, we propose a general Kirchhoff extension formula that handles all the extensions in a similar unified way, and both single-arrival and multi-arrival imaging are discussed. A practical advantage of the proposed method is that all the extensions are constructed using the same piece of code, only input parameters change for different extensions. Similarly, we propose a differential semblance MVA formulation that is simultaneously valid for multiple data-domain extensions.

I.3.9 Acknowledgements

Modeste Irakarama would like to thank Guillaume Caumon and Paul Cupillard for their encouragement and the permission to work on this project. Modeste Irakarama would also like to thank Melchior Schuh-senlis and Mustapha Zakari for many discussions about the algebraic manipulations presented in the appendix. Modeste Irakarama and Thomas Beraud acknowledge the industrial and academic sponsors of the RING-GOCAD Consortium managed by ASGA for their support.

I.3.10 Appendix

The first term in equation 59 is familiar from Plessix (2006). The derivative of the second term in equation 59 leads to

∂ ∂ m s,r δ β -ζ Ḋ = s,r Ḋ ∂ ∂ m δ β -ζ + s,r δ β -ζ ∂ ∂ m Ḋ = s,r Ḋ ∂ ∂ m η ∂ ∂ η δ(η) + s,r δ β -ζ D ∂ ∂ m T , with η = β -ζ = s,r Ḋ ∂ ∂ m η ∂ ∂ η lim σ→0 1 σ √ 2π e -η 2 2σ 2 + s,r δ β -ζ D ∂ ∂ m T = s,r Ḋ ∂ ∂ m η lim σ→0 -η σ 3 √ 2π e -η 2 2σ 2 + s,r δ β -ζ D ∂ ∂ m T = s,r Ḋ ∂ ∂ m ζ lim σ→0 (β -ζ) σ 3 √ 2π e -(β-ζ) 2 2σ 2 + s,r δ β -ζ D ∂ ∂ m T. ( 64
)
We

Part II Notes on Structural Modeling

This second Part of the thesis is devoted to structural modeling, particularly implicit structural modeling (e.g. Chilès et al., 2004;Caumon et al., 2013). The global implicit structural modeling workflow is illustrated in Figure 40: first, input data are interpolated to build a stratigraphic function; then, horizons are extracted from isovalues of the stratigraphic function. In this manuscript, we only discuss the interpolation part, referred to as structural interpolation hereafter. Structural interpolation is usually more challenging than the standard scattered data interpolation (e.g. [START_REF] Turk | Modelling with Implicit Surfaces that Interpolate[END_REF]Anjyo et al., 2014). First of all, structural interpolation has to account for the presence of discontinuities (faults, unconformities, ...) in the most general case. But even in the absence of discontinuities, structural interpolation has additional challenges compared to scattered data interpolation. This is because the output of structural interpolation is a stratigraphic function and is therefore subjected to some geological constraints. The two benchmark models of Renaudeau et al. (2017a) and [START_REF] Renaudeau | The weighted minimization criterion as a physical way to handle thickness variation in implicit structural modeling[END_REF] in Figure 41 are good examples to illustrate some of these challenges. Figures 42a and43a show the results of the scattered data interpolation problem

φ * (x) = argmin φ(x) x Rφ(x) 2 + δ(x -x i ) φ(x) -φ o (x i ) 2 ( 65
)
where R is a linear operator that measures roughness, and φ o (x i ) denotes data observed at the point x i . These results are not suited for implicit structural modeling. The stratigraphic function represents geological time, so it can not have saddle points; this condition is not satisfied by the function obtained in Figure 42a. Furthermore, geological time should increase/decrease monotonically between two successive horizons; this condition is not satisfied by the function obtained in Figure 43a. The more suitable structural interpolation problem may be stated as .

                 φ * (x) = argmin φ(x) x Rφ(x) 2 + δ(x -x i ) φ(x) -φ o (x i )
(66)

The constraint in equation 66b follows from geological field mapping practices. Some geological constraints may be nonlinear. In practice however, we often try to enforce the nonlinear geological constraints using some heuristic linear constraints. For example, the constraint in equation 66c (equation 66b as well) may be satisfied by using an appropriate roughness operator R; this was the case for the result shown in Figure 42b where we use a roughness operator proposed in Irakarama et al. (2018) and in chapter II.1) instead of the Laplacian roughness operator (e.g. [START_REF] Briggs | Machine countouring using minimum curvature[END_REF][START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF] used in Figure 42a. As for the constraint in equation 66d, [START_REF] Renaudeau | The weighted minimization criterion as a physical way to handle thickness variation in implicit structural modeling[END_REF] propose to use a weighted roughness operator where the roughness operator has small weights near data points to allow for rapid variation of the stratigraphic function near data points. Here we explore two new alternative solutions that explicitly take into account the orientation of geological structures via linear constraints, as detailed in chapter II.3; such methods can be quite successful as illustrated in Figure 43b. Chapter II.1

Finite difference implicit structural modeling of geological structures

This chapter resulted from a technical report co-authored with Laurent Gautier, Julien Renaudeau, and Guillaume Caumon.

II.1.1 Summary

We introduce a new method for implicit structural modeling. The main development in this chapter are the new regularization operators we propose by extending inherent properties of the classic 1D discrete second derivative operator to higher dimensions. The proposed regularization operators discretize very naturally on the Cartesian grid using finite differences, owing to the highly symmetric nature of the Cartesian grid. Furthermore, the proposed regularization operators do not require any special treatment on boundary nodes, and their generalization to higher dimensions is straightforward. As a result, the proposed method has the advantage of being simple to implement. Numerical examples show that the proposed method is robust and numerically efficient. Structural implicit modeling, in this chapter, is defined as the interpolation of randomly distributed, and possibly sparse, structural data. For simplicity, structural data, also referred to as structural constraints, will be limited to the following: value constraints : φ(x j ) = a j , (67) orientation constraints : ∇φ(x j ) = ||∇φ(x j )||u j ;

II.1.2 Introduction

(68) where φ(x) is the unknown function to be interpolated, x is a point in space, a is some given scalar, and u is some given unit vector. We refer to equation 68 as the normal form of an orientation constraint, as opposed to the tangential form. Let u 1 be a given unit vector, then there are N -1 unit vectors {u l } l=N l=2 such that {u l } l=N l=1 forms an orthogonal basis in N ≥ 2 dimensions. In that case, the tangential form of the normal orientation constraint ∇φ

(x) = ||∇φ(x)||u 1 is u l • ∇φ(x) = 0 , for l = 2, ..., N. (69)
In implicit structural modeling, the object being modeled is obtained by extracting a hypersurface along an iso-value of the interpolated function φ(x).

Implicit structural modeling has extensively been used in geosciences for contour mapping [START_REF] Briggs | Machine countouring using minimum curvature[END_REF]Mallet, 1984;[START_REF] Smith | Gridding with continous curvature splines in tension[END_REF][START_REF] Wessel | Interpolation with splines in tension: a Green's function approach[END_REF], and for 3D geological modeling [START_REF] Lajaunie | Foliation fields and 3D cartography in geology: Principles of a method based on potential interpolation[END_REF][START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF]. The method presented here was developed mainly for geological modeling applications, but it is suited for contour mapping as well. Most of the groundwork of modern structural implicit modeling was developed in the 2000's [START_REF] Cowan | Rapid geological modeling[END_REF][START_REF] Ledez | Modélisation d'objets naturels par formulation implicite[END_REF]Chilès et al., 2004;[START_REF] Moyen | Paramétrisation 3D de l'espace en géologie sédimentaire: le modèle GeoChron[END_REF][START_REF] Tertois | Création et modification de modéles géologiques par champs de potentiel. Application au modèle GeoChron[END_REF], starting with the pioneering work of [START_REF] Lajaunie | Foliation fields and 3D cartography in geology: Principles of a method based on potential interpolation[END_REF], to the theoretical framework developed by [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF], all the way to the more application-ready works of Chilès et al. (2004) and [START_REF] Frank | 3D reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF]. Structural implicit modeling can be separated into two main classes (Caumon et al., 2013): (1) methods based on dual Kriging and radial basis interpolation [START_REF] Lajaunie | Foliation fields and 3D cartography in geology: Principles of a method based on potential interpolation[END_REF][START_REF] Cowan | Rapid geological modeling[END_REF]Chilès et al., 2004;Calcagno et al., 2008), ( 2) and methods based on domain discretization [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF][START_REF] Frank | 3D reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF]Caumon et al., 2013;[START_REF] Souche | Volume Based Modeling -Automated construction of complex structural models[END_REF]. Domain discretization methods are collectively referred to as Discrete Smooth Interpolation (DSI) methods hereafter [START_REF] Mallet | Discrete smooth interpolation[END_REF][START_REF] Mallet | Discrete smooth interpolation in geometric modeling[END_REF][START_REF] Mallet | Discrete Modeling for Natural Objects[END_REF]. The principle of DSI methods is to discretize all structural constraints on a discrete domain, and assemble them into a least-squares system of linear equations supplemented with smoothing regularization constraints [START_REF] Mallet | Discrete smooth interpolation[END_REF][START_REF] Mallet | Discrete smooth interpolation in geometric modeling[END_REF][START_REF] Mallet | Discrete Modeling for Natural Objects[END_REF].

Implicit geological modeling is still an active area of research. Some authors have shown that existing methods still present limitations caused by the smoothing approaches and the discretization schemes currently used [START_REF] Laurent | Iterative thickness regularization of stratigraphic layers in discrete implicit modeling[END_REF][START_REF] Renaudeau | The weighted minimization criterion as a physical way to handle thickness variation in implicit structural modeling[END_REF]. In addition, recent work has started investigating stochastic modeling approaches of geological structures geometry [START_REF] Jessell | Next generation three-dimensional geologic modeling and inversion[END_REF]Cherpeau and Caumon, 2015;[START_REF] Godefroy | Stochastic association of fault evidences using graph theory and geological rules[END_REF] for which structural interpolation is becoming a bottleneck. Computation times that used to be acceptable for building a single "best model" are becoming far too long when considering stochastic approaches for sampling uncertainty. Other recent research and advances in implicit modeling include better modeling of folds [START_REF] Laurent | Implicit modeling of folds and overprinting deformation[END_REF][START_REF] Grose | Structural data constraints for implicit modeling of folds[END_REF], automated building of models that conform to seismic data (Wu, 2017), more numerically efficient discretization schemes [START_REF] Renaudeau | Implicit modeling of geological structures: A Cartesian grid method handling discontinuities with ghost points[END_REF], and many more (Mallet, 2014;[START_REF] Hillier | Three-Dimensional Modeling of Geological Surfaces Using Generelized Interpolation with Radial Basis Functions[END_REF][START_REF] Gonçalves | A machine learning approach to the potentialfield method for implicit modeling of geological structures[END_REF][START_REF] Martin | Iterative refinement of implicit boundary models for improved geological feature reproduction[END_REF][START_REF] Renaudeau | Implicit Structural Modeling by Minimization of the Bending Energy with Moving Least Squares Functions[END_REF]. As we move towards an era of multi-realization structural modeling [START_REF] Caumon | Towards stochastic time-varying geological modeling[END_REF], new challenges continuously emerge and motivate the quest for more robust and more efficient structural implicit modeling schemes.

In this chapter, we introduce a new approach based on finite differences; this method belongs to the DSI class of methods. Most DSI methods, as applied in geological modeling, rely on simplices (triangles in 2D, tetrahedra in 3D) for domain discretization due to their geometrical flexibility. Because simplicial meshes are unstructured and irregular, it is customary to discretize structural constraints and smoothing constraints using a piecewise linear approach [START_REF] Frank | 3D reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF]Caumon et al., 2013). The method we introduce is based on discretizing the domain on Cartesian grids, and on discretizing structural constraints and smoothing constraints using finite differences. Cartesian gridding and finite differences are arguably the simplest discretizations possible for structural interpolation; as a result, the proposed method is also easy to implement. Modern computing capabilities make the method numerically efficient as well. A 3D implicit structural modeling example using the proposed method is shown in Figure 44.

Implicit structural interpolation on Cartesian grids using finite differences is not new. [START_REF] Briggs | Machine countouring using minimum curvature[END_REF] solves the biharmoinc equation in 2D directly on the Carstesian grid using finite differences. In Briggs (1974)'s formulation, assignment constraints (equation 67) are treated as internal boundary conditions. Briggs also shows the equivalence between the biharmonic equation and the minimization of the global curvature of the implicit function. Wu and Hale (2015) use finite differences operators on Cartesian grids to create implicit functions on seismic images.

Assuming that (1) the normal of reflectors can be estimated at each point of the seismic image, and that (2) the implicit function increases monotonically with seismic traveltime (or depth), they use the normal form of orientation constraints (equation 68) to assemble an overdetermined system solved using least-squares. However, it is usually challenging to make an accurate estimation of reflector dips at each grid point due to coherent noise in seismic data [START_REF] Chauris | Migration velocity analysis from locally coherent events in 2-D laterally heterogeneous media, Part I: Theoretical aspects[END_REF][START_REF] Fehmers | Fast structural intepretation with structure-oriented filtering[END_REF]. As for the second assumption, it is not always valid; the section from the Ribaute model (Caumon et al., 2009) in Figure 45 is an example of where this assumption does not hold. Wu (2017) addresses this problem by first unfaulting and unfolding the seismic image before computing the implicit function. [START_REF] Mallet | Discrete smooth interpolation[END_REF] also proposed an interpolation method on Cartesian grids using finite differences; the method proposed here resolves some issues raised by [START_REF] Mallet | Discrete smooth interpolation[END_REF], as discussed in the section about discontinuities (section II.1.7). Furthermore, unlike [START_REF] Mallet | Discrete smooth interpolation[END_REF], the method proposed here does not require input data to be located at grid points. We propose an alternative approach. [START_REF] Briggs | Machine countouring using minimum curvature[END_REF] noted that the 1D interpolation problem had well defined properties; his strategy was then to extend these properties in 2D by solving the equivalent continuous problem (partial differential equation). Here, we notice that the 1D discrete interpolation problem has inherent properties that we wish to extend to high dimensional discrete interpolation problems. The key new elements of our method are the regularization operators we introduce. These operators discretize naturally on Cartesian grids: they do not require any special treatment at boundary nodes, and their generalization to higher dimensions is straightforward. Let Ω be the domain of definition of φ(x). Ω is discretized on a Cartesian grid, therefore Ω can be seen as a collection of N grid points {x l } l=N l=1 . We distinguish three subdomains in Ω,

Ω = {Ω E ∪ Ω B ∪ Ω I },
where Ω E is the set of external grid points, Ω B is the set of boundary grid points, and Ω I is the set of internal grid points. An example of this discretization is illustrated in Figure 46. A more formal definition of these subdomains will be given shortly. The computational domain Ω C is defined as

Ω C = {Ω B ∪ Ω I }.
(70) Each point in x j in Ω I has neighbors in Ω C , these neighbors define the neighborhood N (x j ) of x j . Figure 47a and Figure 47b illustrate the notion of neighborhoods in 1D and 2D; this figure also introduces directional vectors. For example, the 1D point x j has two neighbors

N (x j ) = {x j-1 , x j+1 } := {x j-d 1 , x j+d 1 };
that is, by definition, the 1D directional vector d 1 is [1]. In general, a point x j always has two neighbors along a given direction: one neighbor in front, and one neighbor behind. For example in 2D the two neighbors of the point x i,j = x k along the direction

d 4 = [1 1] are {x i-1,j-1 , x i+1,j+1 } := {x k-d 4 , x k+d 4 }. ( 71
)
In N dimensions, there are M n = 3 N -1 points in the neighborhood N (x j ) of x j ∈ Ω I , and there are M d = Mn 2 directions such that

N (x j ) = x j±d k for k = 1, ..., M d , (72)
where d k is the k th directional vector, as illustrated for example in Figure 47. Let us come back to the discretization Ω = {Ω E ∪ Ω B ∪ Ω I }. Given a set of points in Ω E , we formally define

Ω I = x j | N (x j ) ∩ Ω E = ∅ . ( 73
) Ω B = x j | N (x j ) ∩ Ω E = ∅ ∧ N (x j ) ∩ Ω I = ∅ . ( 74
)
That is, a grid point that does not have a neighbor in Ω E is by definition a point in Ω I , and a grid point that has at least one neighbor in Ω I and at least one neighbor in Ω E is by definition a point in Ω B . Points in Ω E will usually be specified as inputs. It should be noted that an internal grid point is never in contact with an external grid point: there is always at least one boundary point between them.

II.1.3.2 Problem discretization

The implicit function φ(x) is evaluated on the discrete computational domain Ω C by expansion into the series φ(x) =

x j ∈Ω C B(x, x j )φ(x j ) = x j ∈Ω C B(x, x j )φ j , (75)
where B(x, x j ) are some basis functions, and φ j are the unknown coefficients at grid points.

The basis functions B(x, x j ) are assumed to satisfy the Kronecker delta property, that is for

x i , x j ∈ Ω C , B(x i , x j ) = 1 if i = j 0 if i = j .
Let B(x, x j ) be 1 st -order Lagrange basis functions, simplifying equation 75 to

φ(x) = x j ∈N ′ (x i ) B(x, x j )φ j , (76)
where x i ∈ Ω I is the closest internal grid point to x, and

N ′ (x i ) := {N (x i ) ∪ x i };
it follows from definitions 73 and 74 that x j ∈ Ω C . Using equation 76, the first derivative along the k th direction is given by

∂ k φ(x) = x j ∈N ′ (x i ) ∂ k B(x, x j )φ j . (77)
The constraints 67 and 69 can therefore be discretized in N ≥ 2 dimensions as

φ(x) = x j ∈N ′ (x i ) B(x, x j )φ j = a. ( 78
)
u l • ∇φ(x) = N k u l k ∂ k φ(x) = N k u l k x j ∈N ′ (x i ) ∂ k B(x, x j )φ j = 0, for l = 2, ..., N. (79)
Some constraints may involve second derivatives; this is typically the case for smoothing regularization constraints. However, second derivatives of equation 76 cannot be computed directly, as was the case for first derivatives, because of the use of 1 st -order Lagrange basis functions.

The second derivative at a given point x is therefore approximated by the finite-difference second derivative at x i ∈ Ω I , the closest internal grid point to x. In particular, let d k denote the distance between two neighboring points along the directional vector d k ; then, using the notation in definition 71, the second derivative along the k th direction is given by

∂ 2 k φ(x) ≈ ∂ 2 k φ(x i ) = -2φ i + φ i+d k + φ i-d k d 2 k , (80)
and the mixed derivative along two orthogonal directions k, l is given by

∂ k ∂ l φ(x) ≈ ∂ k ∂ l φ(x i ) = φ i+d k +d l + φ i-d k -d l -φ i-d k +d l -φ i+d k -d l 4d k d l . (81)
In this chapter, all constraints involving second derivatives are imposed strictly on grid points; therefore, the approximation symbol ≈ in equations 80 and 81 can be replaced by the equality symbol =.

II.1.3.3 Solving the discrete problem

The system of equations 78-79, can be written in the more compact matrix form Ā Φ = fa .

(82) Ū Φ = 0 .

(83)

In practice, the system above is usually underdetermined and we have to regularize it with some smoothing constraints of the form R Φ = 0. That is, the solution Φ is obtained by solving

  Ā Ū R  Φ =   fa 0 0   ( 84
)
in a least-squares sense. Convenient choices for R will be proposed later in the chapter. The matrices involved in the system 84 are very sparse and straightforward to build. The structural interpolation problem can be described as a problem of finding a smooth function φ(x) that interpolates a given set of structural data. The smoothness is typically achieved by minimizing the function's roughness R(φ) [START_REF] Mallet | Discrete smooth interpolation[END_REF]. It is unclear what the definition of roughness is, precisely; roughness is usually understood to be related to the notion of curvature.

In 1D, a widely accepted definition of roughness is the second derivate operator [START_REF] Mallet | Discrete Modeling for Natural Objects[END_REF])

R(φ) := ∂ 2 x φ(x); (85)
Minimizing the second derivative in 1D minimizes the curvature of the function; this statement is unambiguous because there is only one way to define curvature in 1D. The 1D smoothing constraint, defined only on grid points, can therefore be discretized using equation 80 to give

R(φ j ) = -2φ j + φ j+1 + φ j-1 |dx| 2 = 0, (86)
for every grid point x j ∈ Ω I .

II.1.4.2 Problem formulation

In higher dimensions, minimizing the roughness of a function becomes ambiguous: there are infinitely many candidate roughness operators in higher dimensions that reduce to equation 85 in 1D. The standard way of overcoming this ambiguity is to explicitly give a physical meaning to equation 85. For example, [START_REF] Briggs | Machine countouring using minimum curvature[END_REF] and [START_REF] Levy | Topologie Algorithmique, Combinatoire et Plongement[END_REF] explicitly look for a function φ(x) that minimizes the global (mean) curvature; in that case, it becomes clear that the extension of equation 85 in higher dimensions is

R(φ) = ∆φ(x),
which is unambiguous since the Laplacian operator ∆ has a precise definition in all dimensions. Another common choice is to seek for a function φ(x) that minimizes the bending energy (see Renaudeau et al., 2017a, and references therein).

In this chapter, we take a different approach to determine the equivalent of equation 85 in higher dimensions. We find inherent properties of the 1D discrete equation 86 that we consider to be practical, from an implementation point of view, then we look for an equivalent discrete equation in higher dimensions which has the same properties. The resulting discrete operator in higher dimensions can then be transformed back to the continuous version if needed.

A very practical property of sparse data interpolation using equation 86 is that it does not require boundary conditions. We define a boundary condition as any constraint imposed on a grid point x ∈ Ω B . In the absence of discontinuities, a regularization matrix R built by imposing equation 86 on all grid points x ∈ Ω I , has a rank( R) = M -2, for a 1D problem with M degrees of freedom (i.e. M grid points in Ω C ). The property rank( R) = M -2 highlights the fact that a minimum of two independent points are required to define a function in 1D. In general, R(φ j ) will be said to satisfy the maximum-rank property if, in the absence of discontinuities, rank( R) = M -(N + 1); where R is the regularization matrix built by imposing R(φ j ) = 0 for all grid points x j ∈ Ω I , for a problem with M degrees of freedom in N dimensions. Our aim is to look for roughness operators R(φ j ) that satisfy the maximum-rank property in higher dimensions.

The maximum-rank property can be extended to higher dimensions if the high dimensional equivalent of equation ( 86) has the following properties:

Property 1. R(φ j ) should include all grid points in the neighborhood of x j , and x j itself; Property 2. R(φ j ) should smooth independently along all directional vectors at x j .

The two properties above are inherent to the 1D discrete problem, as it can be observed by looking at equation ( 86) and Figure 47.a. Therefore, we propose to generalize these two properties to higher dimensions. As we increase dimensions, the first property will guarantee that every degree of freedom is taken into account, and the second property will guarantee that there are "enough" equations to constrain all the degrees of freedom.

II.1.4.3 Towards interpolation in higher dimensions

We propose two strategies for extending equation 86 to higher dimensions. The first one is to generalize the 1D smoothing constraint

R(φ j ) = ∂ 2 x φ j = 0
as a minimization of second directional derivatives of φ j along all directional vectors, that is:

R(φ j ) = ∂ 2 k φ j = 0, for all directional vectors d k at x j . ( 87
)
The second approach is to generalize the 1D smoothing constraint

R(φ j ) = ∂ 2 x φ j = ∂ x (∂ x φ j ) = 0
as a minimization of mixed derivatives of φ j as follows:

R(φ j ) = ∂ k (∂ l φ j ) = 0, for all axis-aligned directional vectors d k , and for all directional vectors d l such that

d l = d k or d l ⊥ d k . ( 88
)
Both equation 87 and equation 88 satisfy the properties imposed in the last paragraph; these two equations lead to different roughness operators. Let us introduce the 2D Cartesian-axes matrix C2 and the 2D diagonal-axes matrix D2 defined as

II.1.5 Interpolation in 2D

C2 = 1 0 0 1 , D2 = 1 √ 2 1 -1 1 1 . ( 89
)
Each row of these matrices is a directional vector on the 2D grid as defined in Figure 47.b. The 2D full direction matrix is defined as

F2 = C2 D2 , (90)
its rows give all the directions available, locally, on a 2D grid.

II.1.5.1 Directional second derivatives

From the directional second derivative point of view (equation 87), the 2D version of equation 86 is

R(φ j ) = ∂ 2 k φ j = -2φ j + φ j+d k + φ j-d k d 2 k = 0, for all d k rows of F2 . (91)
Each of these directional derivatives ∂ 2 k φ j can be written analytically using the 2D Hessian matrix

H2 = ∂ x (∂ x φ) ∂ x (∂ z φ) ∂ z (∂ x φ) ∂ z (∂ z φ) as ∂ 2 k φ = d k H2 d T k .
It follows that a 2D version of the roughness operator 85 that satisfies the maximum-rank property is

R(φ) =                              ∂ 2 [1 0] φ = 1 0 H2 1 0 = ∂ 2 x φ ∂ 2 [0 1] φ = 0 1 H2 0 1 = ∂ 2 z φ ∂ 2 [1 -1] φ = 1 2 1 -1 H2 1 -1 = 1 2 (∂ 2 x + ∂ 2 z -2∂ x ∂ z )φ ∂ 2 [1 1] φ = 1 2 1 1 H2 1 1 = 1 2 (∂ 2 x + ∂ 2 z + 2∂ x ∂ z )φ . ( 92
)
The roughness operator 92 is discretized exactly as the 1D version. Let R denote the regularization matrix resulting from equation 91/92, then the j th row of the least-squares matrix Rt R is a finite-difference approximation of

∂ 2 [1 0] 2 + ∂ 2 [0 1] 2 + ∂ 2 [1 -1] 2 + ∂ 2 [1 1] 2 φ j . (93)
It is interesting to note that equation 93 equates to

∂ 2 x 2 + ∂ 2 z 2 + 2 ∂ x ∂ z ) 2 φ j + 1 2 ∂ 2 x 2 + ∂ 2 z 2 + 2∂ 2 x ∂ 2 z φ j , (94)
where the first term on the right hand side is the bending energy of [START_REF] Enriquez | Applications of bidimensional spline functions to geophysics[END_REF]; [START_REF] Turk | Modelling with Implicit Surfaces that Interpolate[END_REF]; [START_REF] Renaudeau | Implicit Structural Modeling by Minimization of the Bending Energy with Moving Least Squares Functions[END_REF], and the second term is half the squared (mean) curvature of [START_REF] Briggs | Machine countouring using minimum curvature[END_REF], also known as the squared Laplacian [START_REF] Levy | Topologie Algorithmique, Combinatoire et Plongement[END_REF]. Figure 48 shows an example using the roughness operator 91/92 to a modified version of a 2D benchmark model proposed by Renaudeau et al. (2017a).

II.1.5.2 Mixed derivatives

From the mixed derivative point of view (equation 88), another roughness operator R(φ) that satisfies the maximum-rank property is

R(φ) =                                ∂ [1 0] ∂ [1 0] φ = 1 0 H2 1 0 = ∂ x ∂ x φ ∂ [0 1] ∂ [0 1] φ = 0 1 H2 0 1 = ∂ z ∂ z φ ∂ [1 0] ∂ [0 1] φ = 1 0 H2 0 1 = ∂ x ∂ z φ ∂ [0 1] ∂ [1 0] φ = 0 1 H2 1 0 = ∂ z ∂ x φ , (95)
which, using

∂ x (∂ z ) = ∂ z (∂ x ), becomes R(φ) =      ∂ [1 0] ∂ [1 0] φ = ∂ 2 x φ ∂ [0 1] ∂ [0 1] φ = ∂ 2 z φ √ 2∂ [1 0] ∂ [0 1] φ = √ 2∂ x ∂ z φ i . ( 96
)
The factor √ 2 comes from noticing that

∂ [1 0] ∂ [0 1] 2 φ j + ∂ [0 1] ∂ [1 0] 2 φ j = √ 2∂ [1 0] ∂ [0 1] 2 φ j ;
as in the previous paragraph, the reason for equating squares of operators is because the problem is solved by least-squares. The complete sum of squares of terms in this operator equates to the first term on the right hand side of equation 94. The first two terms of equation 95 are discretized using equation 80 and the remaining terms are discretized using equation 81. The 3D Cartesian-axes matrix C3 and diagonal-axes matrix D3 are

II.1.6 Interpolation in 3D

C3 =   1 0 0 0 1 0 0 0 1   , D3 = 1 √ 3     1 -1 -1 1 -1 1 1 1 -1 1 1 1     . (97)
We also introduce the 3D extended-diagonal-axes matrix Ē3 defined as

Ē3 = 1 √ 2        0 1 -1 0 1 1 1 0 -1 1 0 1 1 -1 0 1 1 0        . ( 98
)
The 3D extended-diagonal matrix Ē3 is obtained by inserting an extra zero at different positions of each row of the 2D diagonal matrix D2 . The 3D full direction matrix is defined as

F3 =   C3 D3 Ē3   ; ( 99
)
its rows give all the directions available, locally, on a 3D grid.

Following the same reasoning as in section II.1.5, we use the direction matrix F3 to propose the following choices for the roughness operator:

R(φ) =                                                              ∂ 2 [1 0 0] φ ∂ 2 [0 1 0] φ ∂ 2 [0 0 1] φ ∂ 2 [0 1 -1] φ ∂ 2 [0 1 1] φ ∂ 2 [1 0 -1] φ ∂ 2 [1 0 1] φ ∂ 2 [1 -1 0] φ ∂ 2 [1 1 0] φ ∂ 2 [1 -1 -1] φ ∂ 2 [1 -1 1] φ ∂ 2 [1 1 -1] φ ∂ 2 [1 1 1] φ , and R(φ) =                                                          ∂ 2 [1 0 0] φ i ∂ 2 [0 1 0] φ i ∂ 2 [0 0 1] φ i √ 2∂ [1 0 0] ∂ [0 1 0] φ √ 2∂ [1 0 0] ∂ [0 0 1] φ √ 2∂ [0 1 0] ∂ [0 0 1] φ √ 2∂ [1 0 0] ∂ [0 1 1] φ √ 2∂ [1 0 0] ∂ [0 1 -1] φ √ 2∂ [0 1 0] ∂ [1 0 1] φ √ 2∂ [0 1 0] ∂ [1 0 -1] φ √ 2∂ [0 0 1] ∂ [1 1 0] φ √ 2∂ [0 0 1] ∂ [1 -1 0] φ . ( 100
)
The first operator comes from the directional second derivative approach (equation 87), while the second one comes the mixed derivatives approach (equation 88). It is also possible to combine these two approaches to get a third operator that satisfies the maximum-rank property (section II.1.4.2) ; in particular

R(φ) =                                              ∂ 2 [1 0 0] φ ∂ 2 [0 1 0] φ ∂ 2 [0 0 1] φ √ 2∂ [1 0 0] ∂ [0 1 0] φ √ 2∂ [1 0 0] ∂ [0 0 1] φ √ 2∂ [0 1 0] ∂ [0 0 1] φ ∂ 2 [1 -1 -1] φ ∂ 2 [1 -1 1] φ ∂ 2 [1 1 -1] φ ∂ 2 [1 1 1] φ . ( 101
)
The analytic expressions of each term in these operators are given in the Appendix. While these operators may look complex as we increase dimensions, their discretization remains straightforward: all the second derivative terms are discretized using 80, and all the mixed derivatives terms are discretized using 81. This means that first operator in equation 100 is implemented exactly as its 2D (equation 91), and its 1D (equation 86) version. The only thing that changes are the number of directional vectors and their dimension. Figure 49 shows an example using this smoothing operator to the Balzes model of Ramón et al. (2015). Discontinuities are often encountered in geological structural modeling. They usually consist of faults, unconformities and intrusive bodies. We propose a simple way to handle discontinuities: grid points affected by a discontinuity, after its rasterization, are tagged as external grid points (grid points in Ω E ). Consider for example the illustration in Figure 50: the discontinuity represented by the red line is rasterized to obtain the red external grid points, then boundary grid points are obtained from definition 74. In fact, handling discontinuities in this manner was the main motivation for most of the formalism presented earlier such as the definitions 72, 73, 74, and the maximum-rank property of Section II.1.4.2. This formalism resolves some issues raised by [START_REF] Mallet | Discrete smooth interpolation[END_REF]; [START_REF] Mallet | Discrete smooth interpolation[END_REF] cautioned that the discretization of the roughness operator R(φ) may require special attention at model boundaries and near discontinuities. No special treatment is required at boundaries here, since the roughness operator is only discretized on internal grid points. For the illustration in Figure 50, the solution would be obtained everywhere except at the red grid points. Values at the red grid points may be estimated by extrapolation as a post-processing step. As an example, tagging grid points rasterized by the fault in Figure 51.a as points in Ω E , changes the stratigraphic field from Figure 48.b into the stratigraphic field in Figure 51.b. This straightforward approach has proven capable of handling quite complex fault networks as shown in Figure 52. It also generalizes to higher dimensions as demonstrated by the 3D example in Figure 53. When a layer exhibits large thickness variations, the resulting implicit field can have strong artifacts. This is a problem common to most implicit modeling techniques. According to [START_REF] Smith | Gridding with continous curvature splines in tension[END_REF], this problem seems to come from our need to look for functions with continuous second derivatives everywhere. Current research in implicit structural modeling includes finding a solution to this problem suited for geological modeling applications [START_REF] Laurent | Iterative thickness regularization of stratigraphic layers in discrete implicit modeling[END_REF][START_REF] Renaudeau | The weighted minimization criterion as a physical way to handle thickness variation in implicit structural modeling[END_REF].

II.1.7 Handling discontinuities

II.1.8.2 Discontinuities

The method proposed to handle discontinuities in Section II.1.7 can be numerically inefficient in some cases. If two faults are very close to one another, and there are some data between them, we are forced to use a fine grid in order to take into account those points. This can dramatically increase both the computational time and the computational memory required to solve the system depending on how close the faults are to each other. Current research in finite difference implicit modeling involves finding ways to handle arbitrarily complex fault networks without having to resort to finer grids [START_REF] Renaudeau | Implicit modeling of geological structures: A Cartesian grid method handling discontinuities with ghost points[END_REF].

II.1.8.3 Performance

We use a multi-grid conjugate gradient solver running on GPU on a computer equipped with a 3.5 GHz core and a Quadro M4000 GPU. The 2D problem takes about a second for a problem with 250K grid points, and 3-4 seconds for a problem with 1M grid points. The 3D problem takes about 7-8 seconds for a problem with 1M grid points, and about 40 seconds for a problem with 5M grid points. Future research may include trying to improve performance by preconditioning our solver; some preconditioners for such a problem can be found in Wu and Hale (2015). One limitation of the proposed method is the high memory requirement for 3D problems; if the first roughness operator from equation 100 is used, each grid point has a memory cost of at least 13 equations × 3 nonzeros × (1 row id + 1 column id + 1 value ) × 4 bytes = 468 bytes.

II.1.8.4 Beyond 3D

Sections II.1.5 and II.1.6 gave particular applications of equations 87 and 88 in 2D and in 3D.

While the formulas for the roughness operator R(φ) appeared to become more complex as we went from 1D, to 2D, to 3D, the discrete version of the roughness operator did not change.

The only change was in the size of the full directional matrix F. It is reasonable to expect that this remains true for all higher dimensions N > 3. That is, all we need for an ND interpolation problem is to determine the ND directional matrix FN . In general, the directional matrix is given by

FN =   CN DN ĒN   ,
with the understanding that the diagonal-axis matrix DN only exists for N > 1, and the extented-diagonal-axis matrix ĒN only exists for N > 2. The Cartesian-axis and diagonal-axis matrices are given by

CN (i, j) = δ ij . DN (i, j) = (-1) f(i,j) , with f(i, j) = (2 • i) (j+1-N) .
The extended-diagonal-axis matrix ĒN is obtained by inserting a zero at different positions of each row of the matrix DN-1 ĒN-1 , in a similar way that we obtained Ē3 from D2 in the previous section. FN has 3 N -1 2 rows and N columns.

II.1.9 Conclusion

We have introduced a new technique for implicit modeling on Cartesian grids. We identified inherent practical properties of the well behaved 1D discrete second derivative operator classically used to regularize interpolation in 1D, and then we designed higher dimensional discrete regularization operators with the same properties. In doing so, we obtained regularization operators that discretize very naturally on the Cartesian grid: the operators do not require special treatment on boundary nodes, and they generalize to higher dimensions easily. Discarding boundary conditions allowed us to handle discontinuities by simply surrounding them with boundary nodes. As a result, our method is easy to implement, even in the presence of discontinuities. Numerical experiments demonstrate the robustness and numerical efficiency of the proposed method.

II.1.11 Appendix

The analytic expressions of terms appearing in the formulas for R(φ), as proposed in Section II.1.6, are given below:

∂ 2 [1 0 0] φ = ∂ 2 x φ, ∂ 2 [0 1 0] φ = ∂ 2 y φ, ∂ 2 [0 0 1] φ = ∂ 2 z φ, √ 2∂ [1 0 0] ∂ [0 1 0] φ = √ 2∂ x ∂ y φ , √ 2∂ [1 0 0] ∂ [0 0 1] φ = √ 2∂ x ∂ z φ , √ 2∂ [0 1 0] ∂ [0 0 1] φ = √ 2∂ y ∂ z φ , √ 2∂ [1 0 0] ∂ [0 1 1] φ = ∂ x ∂ y φ + ∂ z φ , √ 2∂ [1 0 0] ∂ [0 1 -1] φ = ∂ x ∂ y φ -∂ z φ , √ 2∂ [0 1 0] ∂ [1 0 1] φ = ∂ y ∂ x φ + ∂ z φ , √ 2∂ [0 1 0] ∂ [1 0 -1] φ = ∂ y ∂ x φ -∂ z φ , √ 2∂ [0 0 1] ∂ [1 1 0] φ = ∂ z ∂ x φ + ∂ y φ , √ 2∂ [0 0 1] ∂ [1 -1 0] φ = ∂ z ∂ x φ -∂ y φ , ∂ 2 [0 1 -1] φ = 1 2 ∂ 2 y + ∂ 2 z -2∂ y ∂ z φ, ∂ 2 [0 1 1] φ = 1 2 ∂ 2 y + ∂ 2 z + 2∂ y ∂ z φ, ∂ 2 [1 0 -1] φ = 1 2 ∂ 2 x + ∂ 2 z -2∂ x ∂ z φ, ∂ 2 [1 0 1] φ = 1 2 ∂ 2 x + ∂ 2 z + 2∂ x ∂ z φ, ∂ 2 [1 -1 0] φ = 1 2 ∂ 2 x + ∂ 2 y -2∂ x ∂ y φ, ∂ 2 [1 1 0] φ = 1 2 ∂ 2 x + ∂ 2 y + 2∂ x ∂ y φ, ∂ 2 [1 -1 -1] φ = 1 3 ∂ 2 x + ∂ 2 y + ∂ 2 z -2∂ x ∂ y -2∂ x ∂ z + 2∂ y ∂ z φ, ∂ 2 [1 -1 1] φ = 1 3 ∂ 2 x + ∂ 2 y + ∂ 2 z -2∂ x ∂ y + 2∂ x ∂ z -2∂ y ∂ z φ, ∂ 2 [1 1 -1] φ = 1 3 ∂ 2 x + ∂ 2 y + ∂ 2 z + 2∂ x ∂ y -2∂ x ∂ z -2∂ y ∂ z φ, ∂ 2 [1 1 1] φ = 1 3 ∂ 2 x + ∂ 2 y + ∂ 2 z + 2∂ x ∂ y + 2∂ x ∂ z + 2∂ y ∂ z φ.

Finite element implicit structural modeling of geological structures

This chapter resulted from an ongoing technical report co-authored with Morgan Thierry-Coudon, Mustapha Zakari, Paul Cupillard, and Guillaume Caumon.

II.2.1 Summary

We introduce a new method for implicit structural modeling. The method is based on a finite element discretization of recently proposed regularization operators for implicit modeling. A finite element implicit modeling scheme is believed to offer some geometrical flexibility as it is readily implemented on both structured and unstructured grids. While implicit modeling on unstructured grids is not new, the method proposed here is, to our knowledge, the first implementation based on finite elements. The finite element method is routinely used to solve boundary value problems. However, because boundary conditions are usually unknown in implicit structural modeling, the traditional finite element method requires minor adjustments in order to be suitable for implicit modeling.

The work presented in this chapter is preliminary and is part of an ongoing master's project by Morgan Thierry-Coudon.

II.2.2 Introduction

Implicit structural modeling is one of the main techniques used in geomodeling [START_REF] Lajaunie | Foliation fields and 3D cartography in geology: Principles of a method based on potential interpolation[END_REF][START_REF] Cowan | Rapid geological modeling[END_REF][START_REF] Frank | 3D reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF]Calcagno et al., 2008;Caumon et al., 2013;[START_REF] Souche | Volume Based Modeling -Automated construction of complex structural models[END_REF]Mallet, 2014). Structural implicit modeling has traditionally been separated into two main classes (Caumon et al., 2013;Wellmann and Caumon, 2018): (1) meshfree methods [START_REF] Lajaunie | Foliation fields and 3D cartography in geology: Principles of a method based on potential interpolation[END_REF][START_REF] Cowan | Rapid geological modeling[END_REF]Chilès et al., 2004;Calcagno et al., 2008;[START_REF] Renaudeau | Implicit Structural Modeling by Minimization of the Bending Energy with Moving Least Squares Functions[END_REF], (2) and mesh-based methods [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF][START_REF] Frank | 3D reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF]Caumon et al., 2013;[START_REF] Souche | Volume Based Modeling -Automated construction of complex structural models[END_REF][START_REF] Renaudeau | Implicit modeling of geological structures: A Cartesian grid method handling discontinuities with ghost points[END_REF]. Mesh-based methods are collectively referred to as Discrete Smooth Interpolation (DSI) methods hereafter [START_REF] Mallet | Discrete smooth interpolation[END_REF][START_REF] Mallet | Discrete smooth interpolation in geometric modeling[END_REF][START_REF] Mallet | Discrete Modeling for Natural Objects[END_REF]. The method presented here belongs to the DSI class of methods. The principle of DSI methods is to discretize all structural constraints on a discrete domain, and assemble them into a leastsquares system of linear equations supplemented with smoothing regularization constraints [START_REF] Mallet | Discrete smooth interpolation[END_REF][START_REF] Mallet | Discrete smooth interpolation in geometric modeling[END_REF][START_REF] Mallet | Discrete Modeling for Natural Objects[END_REF]. In particular, the DSI problem can be stated as that of finding a scalar function φ(x) that satisfies a system of linear equations of the form

D R Φ = fa 0 , (102)
where the data system D Φ = fa is assembled from geological data collected from the field or from seismic images. Geological data are often sparse; as a result, the data system is usually underdetermined and has to be regularized by the regularization system R Φ = 0. We only focus on the regularization system in this chapter.

The Laplacian has been found to be an appropriate regularization operator for some inverse problems (e.g. [START_REF] Clapp | Incorporating geologic information into reflection tomography[END_REF][START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF]. Laplacian regularization is not appropriate for implicit modeling however. This is because the implicit function φ(x) represents stratigraphic time and is therefore subjected to some physical constraints. For example, the gradient of φ(x) should be non-zero everywhere in the domain of study. This constraint is usually not satisfied by Laplacian regularization when interpolating sparse data as the resulting φ(x) may exhibit saddle points. Another issue with Laplacian regularization is that it requires boundary conditions -in this chapter, boundary conditions are defined as any constraint defined on a boundary point. Often, one sets the derivatives of φ(x) to zero at boundaries. As a result, the iso-values of the function φ(x) tend to be perpendicular to boundaries, which is also not desirable in implicit modeling: the iso-values of the function φ(x) represent geological horizons so it is not reasonable to have them artificially terminate perpendicularly on boundaries and geological discontinuities. To avoid the aforementioned shortcomings of Laplacian regularization in implicit modeling, alternative regularization operators have been proposed. For example, [START_REF] Frank | 3D reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF]; Caumon et al. (2013) present an implicit modeling method on unstructured meshes where they assemble a regularization system by minimizing second directional derivatives across tetrahedron faces. In particular, let ∇φ(x) 1 , ∇φ(x) 2 be the gradients of two neighboring tetrahedrons T 1 , T 2 , and n the normal vector of the face shared by T 1 , T 2 , the regularization system is built by imposing

∇φ 1 -∇φ 2 • n = 0 = ∂ nφ 1 -∂ nφ 2 ⇒ ∂ 2 nφ = 0 (103)
across all faces shared by tetrahedrons. The notation ∂ n denotes the directional derivative along the direction vector n. The regularization operator 103 therefore minimizes second directional derivatives along all the normals of tetrahedrons' faces; this is conceptually similar to the regularization operator proposed in chapter II.1, which minimizes second directional derivatives along all the directions on a Cartesian grid. For the general unstructured mesh, the directions of derivation in equation 103 may change from one tetrahedron to another. The directions of derivation on Cartesian grids do not change on the other hand, making it easier to propose an analytical formula for the resulting regularization operator. For example, one of the regularization operators proposed in chapter II.1 in 2D is

R(φ) =            ∂ 2 [1 0] φ = ∂ 2 x φ ∂ 2 [0 1] φ = ∂ 2 z φ ∂ 2 [1 -1] φ = 1 2 (∂ 2 x + ∂ 2 z -2∂ x ∂ z )φ ∂ 2 [1 1] φ = 1 2 (∂ 2 x + ∂ 2 z + 2∂ x ∂ z )φ . ( 104
)
In this chapter, we present a finite element discretization of the regularization operator 104. In practice, both equation 103 and 104 work equally well, so we are merely proposing an alternative method. The finite element method is familiar from solving boundary value problems. However, because boundary conditions are usually unknown in implicit modeling, the traditional finite element method requires minor adjustments in order to be suitable for implicit modeling. The subject of this chapter is to discuss these adjustments.

II.2.3 Finite elements for boundary value problems

In this section, we consider first the case where the regularization system is built by imposing the Laplacian

R(φ) = ∆φ = 0 in Ω I ∇φ • n = 0 on Ω B , (105)
where Ω I is the interior of the domain of interest Ω, Ω B its boundary, such that

Ω = {Ω I ∪ Ω B } , {Ω I ∩ Ω B } = ∅;
n is the normal on Ω B . The short review of finite elements presented here follows the text of [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], to which the reader is referred to for more details. One reduces the order of derivation on φ in equation 105 using integration by parts and the divergence theorem. In particular, let w be a test function, then equation 105 can be transformed to

Ω w(x)∆φ(x) = - Ω ∇w(x)∇φ(x) + Ω B w(x)∇φ(x) • n = 0 (106) = - Ω ∇w(x)∇φ(x) = 0. ( 107
)
In practice, the domain of interest Ω is discretized; we will occasionally think of Ω as a set of discrete points x i . Let N (x, x i ) denote the standard piecewise linear finite element basis function such that

N (x j , x i ) = δ ji , (108)
and express w(x) and φ(x) as

w(x) = x i ∈Ω N (x, x i )w i , (109)
φ(x) = x j ∈Ω N (x, x j )φ j , (110)
allowing equation 107 to be expressed as 0 =

x i ∈Ω w i

x j ∈Ω Ω ∇N (x, x i )∇N (x, x j )φ j = x i ∈Ω w i x j ∈Ω K ij φ j = x i ∈Ω w i G i . (111)
Equation 111 should be valid for all w i ; because the coefficients w i are arbitrary, it follows that

G i = x j ∈Ω K ij φ j = x j ∈Ω Ω ∇N (x, x i )∇N (x, x j )φ j = 0, for all x i ∈ Ω. ( 112
)
The matrix K ij in equation 112 can then be assembled by integrating over each element independently as is customary in finite elements [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF].

II.2.4 Finite elements for structural interpolation

We now consider the finite element discretization of equation 104. We detail only the first term; the rest of the terms are discretized in a similar way. Consider the integral

Ω w(x)∂ 2 x φ(x) = z x w(x)∂ 2 x φ(x) = z - x ∂ x w(x)∂ x φ(x) + w(x)∂ x φ(x) x 2 ∈Ω B x 1 ∈Ω B , (113)
for some two boundary points x 1 , x 2 . In implicit modeling, the derivatives of φ(x) are not known at Ω B , so the boundary term in equation 113 is undesirable. A straightforward way to get rid of that boundary term is to set the derivatives of φ(x) to zero as was done in section II.2.3 for Laplacian regularization. This leads to the matrix

K 1 ij φ j = Ω ∂ x N (x, x i )∂ x N (x, x j )φ j = 0, for all x i , x j ∈ Ω, (114)
where the superscript 1 in K 1 ij refers to the 1 st term of the operator 104. The complete regularization matrix is obtained by gathering the four matrices, one for each term; particularly

K =      K1 K2 K3 K4      . ( 115
)
However, the regularization system built this way is not appropriate for implicit modeling.

Setting the derivatives of φ(x) to zero implies that geological horizons are perpendicular to the model boundary and geological discontinuities; this is too strong of an assumption.

We propose a simple strategy to get rid of the boundary term in equation 113. Instead of imposing heuristic boundary conditions on φ(x), which is subjected to physical constraints, we propose to impose that w(x) vanishes at boundaries. In traditional finite elements, a similar argument is used to get rid of the boundary term on Dirichlet boundaries [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF].

Here we simply extend the argument to the entire boundary. In practice, because the basis functions satisfy the Kronecker delta property 108, imposing w(x) to vanish at boundaries means replacing equation 109 by

w(x) = x i ∈Ω I N (x, x i )w i , (116)
leading to

K 1 ij φ j = Ω ∂ x N (x, x i )∂ x N (x, x j )φ j = 0, for all x i ∈ Ω I , x j ∈ Ω. ( 117
)
Note that the only difference between equations 109 and 116 is the domain of summation. Furthermore, the matrix in equation 117 is not square and symmetric, unlike the matrix in equation 114, since |Ω| > |Ω I |. Note also that using equation 116 with Laplacian regularization would not work for sparse data interpolation because the resulting regularization matrix would not have enough equations (rows); it works for the operator 104 because this operator has multiple linearly independent terms for each grid point. Figure 54 compares interpolation results obtained using the regularization operator 112 with results obtained using the regularization operator 117. Results obtained using the regularization operator 114 are not shown as they are similar to the results shown in Figure 54a. The input data, the white curves, represent two different iso-values. a) Result obtained on a triangle mesh using Laplacian regularization (eq. 112). b) Result obtained on a triangle mesh using the proposed regularization operator (eq. 117).

Equation 116 imposes a minor restriction on mesh generation. Consider for example the illustrative square domain in Figure 55a discretized with nine nodes and eight elements. There is only one internal node: node 5. Because there is no edge connecting the boundary nodes 3 and 7 to node 5, it follows from equation 117 that the 3 rd and 7 th columns of the matrix K (equation 115) will contain only zeros. In order to avoid such cases, each element of the mesh should have at least one internal node as shown in Figure 55b.

II.2.5 Numerical experiments

We now test the different regularization strategies on a section of the Ribaute model of Caumon et al. (2009). The input data is made of three faults, and iso-values picked along three horizons (Figure 56). Faults, and other discontinuities, are handled by duplicating nodes on fault-surfaces as proposed by Caumon et al. (2013). Figure 56a shows results obtained using Laplacian regularization 112 on a triangle mesh; the stratigraphic field tends to be perpendicular to boundaries as may be observed clearly on the bottom boundary. Figure 56b shows results obtained using the regularization operator 114 on a triangle mesh; the artifacts, due to imposing the derivative of the stratigraphic field along boundary normals to zero, are more severe compared to Laplacian regularization. Figures 56c and56d show results obtained using the proposed regularization operator 117 respectively on a triangle mesh and a quadrilateral mesh. The quadrilateral mesh has fewer elements than the triangle mesh for the same resolution. As a result, interpolation on quadrilateral meshes has shown to converge slightly faster. However, for some examples it was not possible to generate a quadrilateral mesh of the entire domain of interest. In particular, the mesh generator was unable to generate quadrilaterals at some fault intersections. a) Result obtained on a triangle mesh using Laplacian regularization (eq. 112). b) Result obtained on a triangle mesh using the regularization operator (eq. 114). c) Result obtained on a triangle mesh using the proposed regularization operator (eq. 117). d) Result obtained on a quadrilateral mesh using the proposed regularization operator (eq. 117).

II.2.6 Discussion & Perspectives

One major challenge of the proposed method is being able to handle complex fault intersections.

The current standard practice of handling faults in implicit modeling on irregular meshes is by duplicating nodes on fault-surfaces. This imposes an explicit computation of fault-intersection points. Furthermore, it imposes nodes on fault-surfaces to carry some topological information, such as which faults does a given node belong to. Such restrictions do not exist in finite difference implicit structural modeling (chapter II.1). When interpolating on unstructured meshes, a preprocessing step is usually required to address those restrictions. The preprocessing step can be long and tedious; ongoing research includes automating this preprocessing step [START_REF] Anquez | A Graph-Based Method to Detect and Correct Invalid Features in Subsurface Structural Models[END_REF][START_REF] Anquez | Automatic correction and simplification of geological maps and cross-sections for numerical simulations[END_REF].

Interpolation on quadrilateral meshes has shown to converge slightly faster than interpolation of triangle meshes; this is in part because a quadrilateral mesh typically has fewer elements than a triangle mesh for the same resolution. However, in some cases, it has proven challenging to generate an entirely quadrilateral mesh in the presence of intersecting faults. One way to take advantage of quadrilateral-elements may be to use a mixed-elements mesh (eg. [START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF][START_REF] Botella | Génération de maillages non structurés volumiques de modèles géologiques pour la simulation de phénom'enes physiques[END_REF]; the method proposed here is readily implemented on mixed-elements meshes.

II.2.7 Conclusion

We present a new implementation of implicit structural modeling based on finite elements of the recently proposed regularization operators (chapter II.1). A finite element implementation is believed to offer some geometrical flexibility as it is readily available for regular and iregular meshes. We show that finite elements for implicit structural modeling requires a minor modification of the standard finite elements for boundary value problems. These modifactions are imposed by the fact that we rarely have boundary conditions in structural modeling; this also imposes a constraint on mesh generation: every boundary node should have at least one edge connecting it to an internal node.

II.2.8 Acknowledgements

Modeste Irakarama's knowledge of implicit modeling has benefited greatly from discussions with Gautier Laurent and Julien Renaudeau. This work also benefited from discussions with François Bonneau, Pierre Anquez and Melchior Schuh-Senlis. Part of the code used in this chapter was written in collaboration with Maxim Torgonskiy and Olivier Giles as a school project while we were still students at École Nationale Supérieure de Géologie. Mesh generation was performed using the open source software Gmsh [START_REF] Geuzaine | Gmsh: A 3d finite element mesh generator with builtin preand postprocessing facilities[END_REF].

This work was done in the frame of the RING project at Université de Lorraine. We would therefore like to thank the sponsors of the RING-GOCAD Consortium managed by ASGA for their support.

Discussion

Figure 57: Example of a highly faulted model. A closer look at fault intersections shows that branching faults may overshot beyond the intersection point in practice. This eventually makes it difficult to generate a mesh that is suitable for implicit modeling. Synthetic data courtesy of ExxonMobil.

In this second Part of the manuscript, we presented new perspectives on implicit structural modeling. In addition to proposing new regularization operators, we introduced Finite Difference Implicit Structural Modeling and Finite Element Implicit Structural Modeling. The main motivation for the work presented in this Part was to develop a simple method for implicit modeling: a method based on finite differences on Cartesian grids. The finite element method was developed later following the interest of some colleagues in implementing the new regularization operators on irregular meshes. Current standard grid-based implicit modeling methods rely on irregular meshes. The conceptual advantage of irregular meshes over regular meshes is that irregular meshes require much less elements, and therefore less computer memory, to accurately model geometries with curves. The disadvantage of irregular meshes is that it becomes a challenge to generate a mesh in complex models with discontinuities. Consider for example the highly faulted model in Figure 57. A closer look at fault-fault intersections shows that some branching faults overshoot beyond the branch point. Because current practice for handling faults in implicit modeling on irregular meshes involves duplicating fault nodes (e.g. Caumon et al., 2013), one has to explicitly compute the fault-fault intersection point and discard the overshooting part in order to generate a suitable mesh. This preprocessing step usually involves a lot of manual manipulations and can be quite fastidious, especially in 3D. Furthermore, computing fault branch lines in 3D is a challenging task on its own. The preprocessing step is not needed for the finite difference implicit modeling proposed in chapter II.1 (see Figure 52 for an example using the model in Figure 57). Current implicit modeling methods perform poorly in the presence of large thickness variations (e.g. [START_REF] Collon | 3D Geomodelling in Structurally Complex Areas -Implicit vs. Explicit representations[END_REF][START_REF] Renaudeau | The weighted minimization criterion as a physical way to handle thickness variation in implicit structural modeling[END_REF]. In particular the constraint in equation 66d is usually violated in the presence of large thickness variation as illustrated by the results in Figure 58a. We now discuss two new strategies to address this problem; both methods use the orientation of horizons to guide the interpolation. The first step is common to both methods, it consists in estimating the orientation of the stratigraphic field in the entire domain from input data as shown in Figure 58b. The first method we propose is simply to constrain the stratigraphic field to follow the estimated vector field using equation 69; this results in the stratigraphic field in Figure 58c. The second method we propose is to use the estimated vector field to build a regularization operator that smooths preferentially along the estimated vector field; in particular, the regularization operator in equation 92 is modified to

Rφ(x) =              w 1 (x)∂ 2 d1 φ = w 1 (x)∂ 2 [1 0] φ = w 1 (x)∂ 2 x φ w 2 (x)∂ 2 d2 φ = w 2 (x)∂ 2 [0 1] φ = w 2 (x)∂ 2 z φ w 3 (x)∂ 2 d3 φ = w 3 (x)∂ 2 [1 -1] φ = w 3 (x) 1 2 (∂ 2 x + ∂ 2 z -2∂ x ∂ z )φ w 4 (x)∂ 2 d4 φ = w 4 (x)∂ 2 [1 1] φ = w 4 (x) 1 2 (∂ 2 x + ∂ 2 z + 2∂ x ∂ z )φ , (118)
with

w i (x) = di • s(x) ||di||||s(x)|| α + ǫ, ( 119
)
where di is any of the directions of derivation in equation 118, s(x) is the value of the estimated vector field (Figure 58b) at the point x, ǫ is a small value to avoid zero weights, and α is a parameter that determines the degree of anisotropic smoothing. The results in Figure 58d were obtained using α = 10.

The first method proposed (Figure 58c) is conceptually similar to that of [START_REF] Laurent | Iterative thickness regularization of stratigraphic layers in discrete implicit modeling[END_REF] who also proposed to estimate a vector field from input data and use the vector field to constrain the interpolation. Our approach differs from that of [START_REF] Laurent | Iterative thickness regularization of stratigraphic layers in discrete implicit modeling[END_REF] at least in two ways:

1. Laurent (2016) estimates the stratigraphic orientation vector field (e.g. Figure 58b) from the gradient of an initial scalar field estimated from input data (e.g. Figure 58a) and then updates both the vector field and the scalar field iteratively; this approach is not guaranteed to converge and is computationally intensive since it is iterative. Here, we use every two consecutive points along input isovalue data to define vectors. The resulting vectors are then interpolated componentwise in the entire domain.

2. Laurent (2016) uses the estimated vector field to constrain the norm of the gradient of the stratigraphic field while we use estimated vector field to constrain the orientation of the stratigraphic field. Constraining the norm of the gradient is a nonlinear problem whereas constraining the orientation is a linear problem.

The differences between our method and that of Laurent ( 2016) can be justified by the difference in the problems we are seeking to solve: Laurent (2016) seeks to enforce the constraint in equation 66b while we seek to enforce the constraint in equation 66d.

The second method proposed (Figure 58d) reduces the artifacts observed in Figure 58a considerably but it does not eliminate them entirely: a closer look shows that artifacts are still present in Figure 58d (e.g. near the point 170,1100). It should also be noted that equation 119 is a specific case of the more general

w i (x) = M di, s(x) ,
where M (•, •) is an operator that measures the difference between two vectors. Current research includes defining an appropriate operator M that may eliminate the remaining artifacts in Fig- ure 58d.

One disadvantage common to the two methods proposed is the need to specify the orientation of each input isovalue curve used to estimate the stratigraphic vector field (e.g. Figure 58b); this may not always be straightforward, especially in 3D.

Part III

Reducing Structural Uncertainties using Seismic Data This third Part of the manuscript is devoted to reducing structural uncertainties by appraising structural models using seismic data. This problem is illustrated in Figure 59: a single seismic image supports multiple structural interpretations; can we use seismic data to determine which models are more likely than others?

Let M denote a structural model, and let M denote an arbitrary set of candidate structural models. Consider the optimization problem

         M * = argmin M Φ d (M), M ∈ M (a)
subject to M * is geologically consistent (b) M * is the closest model in M to the "true subsurface geology" (c)

, ( 120
)
where Φ d (M) denotes a seismic data misfit value associated to the structural model M. A solution to problem 120 implies a solution to the appraising structural models using seismic data problem: let Φ m (M), the model misfit, be a measure of how close the structural model M is to the true subsurface geology; then, if problem 120 has a solution, it follows that

∞ > Φ d (M j ) -Φ d (M i ) Φ m (M j ) -Φ m (M i ) > 0 when i = j, (121)
for any pair of models M i , M j ∈ M (this is certainly true for a set that has only two candidate models {M 1 , M 2 }; because these two models are arbitrary, it follows that it's also true for any two pairs of models M i , M j ∈ M ). Given a definition of Φ m , the problem of appraising structural models using seismic data will consist in defining seismic data misfit functions Φ d that satisfy equation 121 for all pairs of models M i , M j ∈ M . Equation 121implies that ranking models using Φ m is equivalent to ranking models using Φ d ; therefore, if we can satisfy that equation we can rank models using seismic data (Φ d ) without any knowledge of the true subsurface geology needed to define Φ m .

Problem 120 has several challenges:

1. The data misfit Φ d in equation 120a implies that there exists a forward modeling operator that generates synthetic seismic data from a structural model. In this work, structural models are populated with velocities estimated from the migration velocity model, which are then used to run numerical simulations. Populating structural models with velocities from a smooth migration velocity model is not a trivial problem.

2. The constraint in equation 120c implies that we can define a distance between two models, denoted Φ m ; this is also not a trivial problem.

3. Given a definition of Φ d and Φ m , we still have the problem of satisfying the constraint in equation 120c, or equivalently satisfying equation 121 for any pair of models M i , M j ∈ M ; this is also not a trivial problem.

In this thesis, we will only propose straightforward solutions to the first and second challenge so that we can focus on the third challenge. The constraint in equation 120b will be satisfied by requiring that all structural models in M be constructed from structural interpretations proposed by geologists.

Problem 120 is conceptually similar to methods that have been proposed in reservoir modeling to determine reservoir models that fit production data (e.g. [START_REF] Suzuki | Dynamic data integration for structural modeling: Model screening approach using a distance-based model parameterization[END_REF][START_REF] Cherpeau | Method for Stochastic Inverse Modeling of Fault Geometry and Connectivity Using Flow Data[END_REF]. What is new in problem 120, other than the use of seismic data, is the constraint in equation 120c. It is usually implicitly assumed that a model that minimizes data misfit is also the "most correct" model. However, this assumption is not necessarily true, especially in geology where we usually solve nonlinear physical problems in complex models. For example, [START_REF] Carter | Our calibrated model has poor predictive value: An example from the petroleum industry[END_REF] show that reservoir models that match available production data may have poor predictive values.

Chapter III.1

Structural Interpretation Uncertainties

We speak of structural uncertainties when the position and/or the presence of a structural surface is uncertain on a seismic image (Thore et al., 2002). Sources of uncertainties in seismic imaging include data uncertainty and velocity uncertainty (e.g. [START_REF] Li | Assessing seismic uncertainty via geostatistical velocity-model perturbation and image registration: An application to subsalt imaging[END_REF][START_REF] Pawelec | UNCERTAINTY QUANTIFICATION IN SEISMIC IMAGING[END_REF], acquisition geometry (e.g. [START_REF] Lecomte | 2(3)D convolution modeling of complex geological targets -Beyond 1D convolution[END_REF], and geological complexity (e.g. [START_REF] Suiter | Reducing structural uncertainties through anisotropic pre-stack depth imaging: examples from the Elgin/Franklin/Glenelg HP/HT fields area, Central North Sea[END_REF]. As an example of uncertainties due to geological complexity, consider the Caumon salt model in Figure 60a. The model is relatively simple from a structural geology point of view but the presence of salt and thin layer structures makes the model challenging for seismic imaging.

The model was built by adding thin layer structures in the macrolayered velocity model in Figure 60b. The seismic image computed from the thinlayered model (Figure 60c) has a more realistic texture than the seismic image computed in the macrolayered model (Figure 60d); therefore, the true subsurface geology is expected to look more like the model in Figure 60a than the model in Figure 60b [START_REF] Landa | Realistic Finite Differences Modeling -A Case Study[END_REF]. When building structural models however, we are usually only interested in major layers (see for example the structural model in Figure 62). The problem of structural interpretation is therefore that of identifying the structural surfaces (faults, horizons, ...) observed in the macrolayered model in Figure 60b from the seismic image in Figure 60c. This is not entirely straightforward even for this experiment where we used the exact velocity model for seismic imaging. Consider for example the well illuminated portion of the model between the distances x = 2000m and x = 6000m in Figure 60c. While it is fairly straightforward to identify the fault, it is less obvious how to identify the major horizons and to correlate them across the fault.

When the positions and/or the presence of structural surfaces are uncertain on a given seismic image, it is usually possible to propose multiple structural interpretations from the same seismic image (e.g. Bond et al., 2007). Consider for example the seismic image in Figure 61; this image was generated from the sandbox model of Colletta et al. (1991) (see Figure 70 in chapter III.2) using the fine layering technique described in the previous paragraph. The image was then given to different interpreters who provided 7 different structural interpretations, shown in Figure 63 to Figure 69. The structural interpretation in Figure 62 was performed directly on the reference velocity model (Figure 70) in order to have one structural interpretation that we are sure is the most accurate, for benchmarking purposes. The different interpretations are geologically consistent, but because they are different, we can expect some interpretations to be more consistent than others. In this third Part of the manuscript, we are going to study the problem of determining which of those structural interpretations are more consistent than others using seismic data. In particular, we will generate synthetic data for each structural interpretation, and then compare the synthetic data with observed data in order to assign a misfit value to each interpretation. Chapter III.2

Appraising Structural Models using Seismic Data -Theory

This chapter resulted from a technical report co-authored with Paul Cupillard, Guillaume Caumon, Paul Sava, and Jonathan Edwards.

III.2.1 Summary

Structural interpretation of seismic images can be highly subjective, especially in complex geological settings. A single seismic image will often support multiple geologically valid interpretations. However, it is usually difficult to determine which of those interpretations are more likely than others. We refer to this problem as structural model appraisal herein. In this chapter, we use misfit functions to rank and appraise multiple interpretations of a given seismic image. Given a set of possible interpretations, we compute synthetic data for each structural interpretation, then compare these synthetic data against observed seismic data; this allows us to assign a data-misfit value to each structural interpretation. Our aim is to find data-misfit functions which enable a ranking of interpretations. To do so, we formalize the problem of appraising structural interpretations using seismic data and derive a set of conditions to be satisfied by the data-misfit function for a successful appraisal. We investigate both vertical seismic profiling and surface seismic configurations. An application of the proposed method to a realistic synthetic model shows promising results for appraising structural interpretations using vertical seismic profiling data, provided the target region is well illuminated. However, we find appraising structural interpretations using surface seismic data to be more challenging, mainly due to the difficulty of computing phase-shift data-misfits.

III.2.2 Introduction

Uncertainties in structural interpretations of seismic images have been of interest to geoscientists for a long time. Earlier studies were mainly focused on migration velocity errors and seismic resolution; some examples can be found in [START_REF] Hajnal | Maximum uncertainty of interval velocity estimates[END_REF] and [START_REF] Beylkin | Spatial resolution of migration algorithms[END_REF]. Structural uncertainties propagated from velocity errors usually translate to subtle displacement or flexing of horizons and faults away from the reference position of those horizons and faults [START_REF] Thore | A practical formulation of migration errors due velocity uncertainties[END_REF][START_REF] Bube | The influence of stacking velocity uncertainties on structural uncertainties[END_REF][START_REF] Pon | Sensitivity analysis of seismic depth migrations[END_REF][START_REF] Fomel | Structural uncertainty of time-migrated seismic images[END_REF]. By the reference position of a given horizon or fault, one often means the position obtained from the tomographic velocity model (e.g. Messud et al., 2017). Geological realizations obtained by such perturbations of a reference structural model often underestimate structural uncertainties, since interpretation uncertainties are not taken into account. Interpretation uncertainties arise where reflectors cannot be tracked deterministically; this can be a result of an inaccurate imaging velocity model [START_REF] Li | Assessing seismic uncertainty via geostatistical velocity-model perturbation and image registration: An application to subsalt imaging[END_REF], poor illumination and/or poor resolution [START_REF] Lecomte | 2(3)D convolution modeling of complex geological targets -Beyond 1D convolution[END_REF]. In these situations, multiple geologically possible structural models can be interpreted from the same seismic image (Bond et al., 2007;Bond, 2015;[START_REF] Alcalde | Impact of seismic image quality on fault interpretation uncertainty[END_REF].

Generating multiple structural models from a single data set was one of the earliest methods proposed to evaluate structural uncertainties in geomodeling (Thore et al., 2002). This method opened ways to new challenges, two of which are :

1. Sampling structural uncertainties. This is essentially a structural modeling problem [START_REF] Abrahamsen | Bayesian kriging for seismic depth conversion of a multi-layer reservoir[END_REF][START_REF] Lecour | Modelling of stochastic faults and fault networks in a structural uncertainty study[END_REF][START_REF] Holden | Stochastic Structural Modeling[END_REF]Wellmann et al., 2010[START_REF] Wellmann | Validating 3-D structural models with geological knowledge for improved uncertainty evaluations[END_REF]Cherpeau and Caumon, 2015;[START_REF] Julio | Sampling uncertainty about segmented normal fault interpretation using a stochastic downscaling method[END_REF]. In particular, how can we generate multiple structural models respecting both uncertainty bounds and geological rules efficiently?

2. Appraising structural models. This usually amounts to evaluating the model's response for some observed physical phenomenon. In most studies on structural model appraisal, the physical phenomenon that has been used to evaluate the model is fluid flow [START_REF] Suzuki | Dynamic data integration for structural modeling: Model screening approach using a distance-based model parameterization[END_REF][START_REF] Cherpeau | Method for Stochastic Inverse Modeling of Fault Geometry and Connectivity Using Flow Data[END_REF]. Potential fields (gravity and magnetic) have also been used to validate model scenarios [START_REF] Foss | Geologically constrained seismic imaging -Workflow[END_REF][START_REF] Fullagar | Constrained inversion of geologic surfaces -pushing the boundaries[END_REF].

With standard geomodeling workflows, we are able to generate multiple geological scenarios from a seismic image, but we rarely go back to check if these models are consistent with the initial data; doing so would be a way to reduce structural uncertainty. Here, initial data refers to the recorded pre-migrated seismic traces. An interesting approach proposed by [START_REF] Lecomte | Simulated Prestack Local Imaging: a robust and efficient interpretation tool to control illumination, resolution, and time-lapse properties of reservoirs[END_REF] to address this issue is to generate synthetic seismic images from different geological scenarios and compare them to the initial seismic image; a similar approach was adopted by [START_REF] Lallier | Relevance of the stochastic stratigraphic well correlation approach for the study of complex carbonate settings: application to the Malampaya buildup (Offshore Palawan, Philippines)[END_REF]. A major challenge of this approach is that it is not clear how to objectively compare seismic images. An alternative way for reducing interpretation uncertainties in structural modeling is to use "data-driven" interpretation approaches. Here, data-driven interpretation includes all the methods that have been proposed to automatically track horizons and faults, and to build a relative-geologic-time function from a seismic image [START_REF] Stark | Relative geologic time (age) volumes-Relating every seismic sample to a geologically reasonable horizon[END_REF][START_REF] Pauget | A global approach in seismic interpretation based on cost function minimization[END_REF]Wu and Hale, 2015;Wu, 2017). Data-driven techniques, however, require good quality seismic data [START_REF] Hoyes | A review of global interpretation methods for automated 3d horizon picking[END_REF]; where the geology is complex or lacks reflectivity contrasts, manual interpretation is still necessary.

In this chapter, we aim to reduce interpretation uncertainties by appraising structural interpretations using seismic data. Our approach is conceptually similar to that proposed by [START_REF] Lecomte | Simulated Prestack Local Imaging: a robust and efficient interpretation tool to control illumination, resolution, and time-lapse properties of reservoirs[END_REF]: we generate synthetic data from a set of candidate geological models, then we compare the synthetic data against observed data. The candidate geological models are obtained from different interpretations of the same seismic image. There are two main differences between the method proposed here and the method proposed by [START_REF] Lecomte | Simulated Prestack Local Imaging: a robust and efficient interpretation tool to control illumination, resolution, and time-lapse properties of reservoirs[END_REF]. First, we compare initial pre-migrated data instead of seismic images; second, we propose a quantitative comparison instead of a qualitative comparison. The objective of this chapter is to investigate how exactly we should compare synthetic data against observed data. We mainly focus on theoretical aspects of the problem and defer the more practical issues for further investigations.

In what follows, we will introduce the problem of appraising structural interpretations using seismic data, and propose a mathematical formulation to describe the process. We then propose some strategies to build a velocity model with structural discontinuities from a structural interpretation and a (smooth) migration velocity model; the resulting velocity model with structural discontinuities is then used to generate synthetic data for the given structural interpretation.

Finally, we propose a workflow to design data misfit functions in order to evaluate and rank the different structural interpretations. The proposed method is then applied on a realistic synthetic vertical seismic profiling (VSP) case.

III.2.3 The structural interpretation appraisal problem

We use the model in Figure 70 to illustrate an example of interpretation uncertainty. The model was built by assigning constant velocity values in each major layer of the sandbox model of Colletta et al. (1991). Additional thin-layering structures were added inside each major layer in order to generate more realistic synthetic seismic data [START_REF] Landa | Realistic Finite Differences Modeling -A Case Study[END_REF]. The resulting velocity model was then used to simulate "observed" seismic data using an acoustic, constant density, staggered finite-difference scheme (Virieux, 1984) with perfectly matched absorbing boundaries (Collino and Tsogka, 2001). All subsequent seismic modeling mentioned in this chapter were performed with the same code. The reference velocity model in Figure 70 was also used to obtain the migration velocity model in Figure 71a by Gaussian smoothing of the slowness; this migration velocity model was then used to migrate the "observed" data to obtain the seismic image in Figure 71b using Kirchhoff depth migration (e.g. Etgen et al., 2009).

Figures 72a and 72b

show two different structural interpretations of the seismic image in Figure 71b from two different interpreters. Both interpretations are geologically possible, but they are different in some regions. Assuming that one is more accurate than the other in the regions where they differ, can we use seismic data to determine which of the two models is more accurate in those regions? To answer this question, we rely on macro-layered velocity models. Macro-layered velocity models are defined as velocity models with discontinuities from structural models and velocity values from the migration model (e.g. Figure 72c and Figure 72d). The migration velocity model is assumed to be kinematically accurate throughout this chapter.

We use these macro-layered velocity models to compute synthetic data, which are compared against the observed data to appraise the different geological scenarios. Structural interpretations are typically segments or curves picked by interpreters along struc-tural discontinuities on a seismic image. A structural model is a set of consistent structural surfaces, such as horizons and faults, that represent a geological model (Caumon et al., 2009) (e.g. Figure 72a and Figure 72b). In this chapter, structural models are built from structural interpretations and there is a one-to-one relation from structural interpretations to structural models. Therefore, appraising structural interpretations is synonymous with appraising structural models herein. The relation from structural models to macro-layered velocity models is assumed to be one-to-one as well. Let us define a model-misfit as a number that quantifies the mismatch between the reference structural model and a candidate structural model. The reference structural model is the "true" structural model; it is unknown in practice. However, for the sake of argument, we assume the reference model to be known for now. Let us also assume that we can determine this modelmisfit value for each model in a given set of candidate structural models. Furthermore, we assume that different models will have different model-misfit values, allowing us to rank those models from the best model to the worst model. This ranking of structural models using model-misfit values will be referred to as model-space ranking of structural models. Objects in the model-space are structural interpretations, or alternatively structural models. We also define a data-misfit as a number that quantifies the mismatch between observed data and the synthetic data generated from candidate structural models. Data-misfit values will allow us to rank structural interpretations in the data-space; we refer to this as data-space ranking. The data-space here refers to the pre-migration domain; i.e., objects in the data-space are seismic traces. We assume that the observed data are the wave equation's response of the reference model (true earth) subjected to some acquisition geometry, just like synthetic data are the wave equation's response of a given candidate (synthetic) model. If the ranking in data-space is the same as the ranking in model-space, then we do not need to know the reference model to determine which among our candidate models are more probable than others. For this reasoning to be applicable, the following conditions have to be satisfied:

1. Different candidate structural models should, overall, have different data-misfit values.

2. Data-misfit values should, overall, be rank-correlated to model-space misfit values.

In practice, a pair of different structural models may exhibit similar seismic response at receivers; in this case, the first condition can be satisfied by ignoring one of the models. The second condition means that if we plot model-misfits and data-misfits on a scatter plot, we should be able to find a monotonic curve that fits the data.

In general, it is straightforward to design data-misfit functions and compute data misfit-values. However, it is much harder to design data-misfits that satisfy the conditions mentioned above. The objective of this chapter is to investigate what it takes to define data-misfit functions that honor the conditions mentioned above as much as possible, and thus run the process with real data, where the reference (true) model truly is unknown.

III.2.4 Mathematical formulation

Let M i denote a candidate structural model (e.g. Figure 72a and Figure 72b), the index i taking different values for different interpretation scenarios. Each structural model M i can be used, together with the migration velocity model, to build a macro-layered velocity model, also denoted by M i . Examples of macro-layered velocity models will be presented in the next section (see also Figure 72c and Figure 72d).

Let us introduce the model-space misfit function Φ m i,j that measures the difference between two structural models M i and M j , and denote the mismatch between the reference model and a candidate model

M i as Φ m i = Φ m i,ref .
Similarly, we introduce a data-space misfit Φ d i,j that measures the difference between data computed in two macro-layered models M i and M j , and denote the mismatch between data computed in the macro-layered model M i and the observed (reference) data as

Φ d i = Φ d i,ref . Two models M i and M j are said to be data-different if |Φ d i -Φ d j | > ǫ,
for some "small" user-defined threshold ǫ > 0.

Let M be a set of data-different structural models. If Φ d i and Φ m i are perfectly rank-correlated for all models in M , then ranking models in M using Φ d i is the same as ranking models in M using Φ m i . In that case, by definition, the structural models in M can be appraised (i.e. ranked) using seismic data. Therefore, given a set of macro-layered velocity models M , we consider the problem of appraising structural models using seismic data to be solvable in M if the following conditions are satisfied:

Condition 1. |Φ d j -Φ d i | > ǫ when i = j. Condition 2. ∞ > Φ d j -Φ d i Φ m j -Φ m i > 0 when i = j.
Given a set of random structural models R, we can build a set of data-different structural models M by progressively moving each model M i from R to M , when M i is data-different from all the models already in M . Therefore, condition 1 can be satisfied in practice. We focus on condition 2 hereafter. Let us define Φ d i (Ω) and Φ m i (Ω) respectively as the data-space and model-space misfit values corresponding to a specific region Ω in the candidate model M i . The reason for localizing misfits in space will be justified later. In that case, two necessary conditions, although not sufficient, for (the localized version of) condition 2 to be satisfied are:

Condition 3. A change of Φ d i (Ω) ⇒ a change of Φ m i (Ω). Condition 4. A change of Φ m i (Ω) ⇒ a change of Φ d i (Ω)
. The symbol ⇒ stands for "implies". Furthermore, suppose that Φ m satisfies the following condition:

Condition 5. A change of Φ m i (Ω) ⇔ a change of the structural model M i in Ω.
Then conditions 3 and 4 become:

Condition 6. A change of Φ d i (Ω) ⇒ a change of the structural model M i in Ω. Condition 7. A change of the structural model M i in Ω ⇒ a change of Φ d i (Ω)
. At this point, we have transformed the problem of appraising structural models using seismic data to that of satisfying one condition on Φ m (condition 5) and three conditions on Φ d (conditions 1, 6, and 7). The conditions imposed on Φ d can be checked in practice as they do not depend on the (unknown) reference model. Note that condition 2 implies that we assume the relation between the model and data to be weakly nonlinear1 , i.e., we assume the perturbations in conditions 3 and 4 to be "small". In other words, we assume, to some extent, that all the candidate structural models are "close" to the (unknown) reference model. Our goal will be to design Φ d such that conditions 6 and 7 are honored as much as possible. We will argue that for VSP data-misfits, condition 7 is relatively easy to satisfy whereas condition 6 is practically impossible to satisfy in general. As for surface seismic data-misfit functions, we will argue that if all the candidate macro-layered velocity models are kinematically equivalent to the migration velocity model, itself assumed to be kinematically accurate, both conditions 6 and 7 can be satisfied in theory for stratified (i.e. layered) models; however, it remains challenging to design and compute misfit function satisfying these conditions. These difficulties make the ideal condition 2 too strict in practice; we replace it with the less ambitious condition:

Condition 8. Φ d is statistically rank-correlated with Φ m .
Condition 2 means that the rank correlation coefficient of the variables Φ m and Φ d is 1, while condition 8 means that it is "high enough".

III.2.5 Macro-layered velocity models

We propose several options for building macro-layered velocity models from a structural interpretation and a migration model. In this section, we use the illustrative model in Figure 73a as the reference model, and the model in Figure 73b as the migration model; the migration model was obtained by smoothing the slowness of the reference model.

III.2.5.1 Block macro-layered velocity model

The first type of macro-layered model we propose is obtained by averaging migration slowness values (Figure 73b) in each layer of a structural model. Averaging slowness values preserves traveltimes, as opposed to averaging velocities. The velocity of the i th layer is therefore given by

v i = mean 1 M mig (x)
x in layer i .

This results in a block velocity model, denoted M b hereafter, as illustrated in Figure 74a.

Figure 74b shows velocity model residuals obtained by subtracting the reference model from the block macro-layered model. Well-log information, when available, can also be used to constrain velocity values assigned in each layer of block macro-layered models. Block macrolayered velocity models are not in general kinematically equivalent to the migration velocity model.

III.2.5.2 Wire macro-layered velocity model

The second type of macro-layered model is the wire velocity model, M w , illustrated in Figure 75. This model is obtained by

M w (x) = M mig (x) + ∆M b (x), (123)
where ∆M b (x) denotes the velocity contrast at point x; it is equal to zero everywhere except at block boundaries where it is equal to the velocity jump between the velocities in adjacent blocks. Wire macro-layered velocity models are kinematically equivalent to the migration velocity model everywhere expect at geological interfaces.

III.2.5.3 Reflectivity macro-layered velocity model

We also propose the reflectivity velocity model, M r , illustrated in Figure 76 and obtained by

M r (x) = M mig (x) + M b (x) -Mb (x) , (124)
where

III.2.5.4 Density macro-layered velocity models

Finally, we propose a variable density macro-layered model, M d , given by the pair (M mig , M * ) where M * is a density model estimated from M b , by Gardner's law for example [START_REF] Gardner | Formation velocity and density-the diagnostic basics for stratigraphic traps[END_REF]. The variable density model allows us to propose different macro-layered velocity models that are kinematically equivalent to the migration model by proposing different impedance models.

III.2.6 Data-space misfit functions

Seismic data can only contain structural information about subsurface regions that have been illuminated by the wavefield recorded at receivers. Therefore, a data-misfit value computed for a given structural interpretation will only be representative of the illuminated region, and not of the entire model. Furthermore, a misfit value is a scalar, so it summarizes all the interpretation errors in the model in just one value. As a consequence, the best information a global misfit can offer is to determine which of two models is better than the other, but not why (i.e., where the errors are coming from). This motivates the need to look for misfit functions that are localized in space, i.e., data-misfit functions that are a function of space. It is important that we know where the errors are coming from in our model because we might be able to update our structural interpretation in that region to lower the misfit. With this in mind, we propose a data-misfit function of the form

Φ d i (Ω) = s,r x∈Ω ω(s, r, x)χ f o (s, r, t), f i (s, r, t) , ( 125
)
where Ω is the region in the model we would like to evaluate, ω is a weighting function depending on the illumination in Ω, f o (s, r) is an observed seismogram from a shot fired at x s = s and recorded at a receiver at x r = r, f i (s, r) is the corresponding synthetic seismogram computed in the candidate macro-layered velocity model M i , and χ is a data residual function computing the mismatch between observed and synthetic data. We will consider one region of interest Ω at a time; this region can have an arbitrary shape, provided it is illuminated by at least one source-receiver raypath. The data-misfit 125 can be written in the more compact form

Φ d i (Ω) = x∈Ω W χ f o (s, r, t), f i (s, r, t) , (126)
where it is now more apparent that the data-misfit function 126 acts in two steps: first, the residual function χ maps the data plane D×D to a reduced data-space D ′ ; second, the weighting operator W maps the reduced data-space D ′ to the image-space I. Here, the image-space refers to the range of any seismic imaging operator; it is usually referred to as the image-domain. The reduced data-space is, by definition, the range of the residual function χ; a concrete example of a reduced data-space will be presented in the following paragraph. Therefore, if we can find an appropriate linear relation L that maps I to D ′ , the weight function ω can be determined from L * : D ′ → I, where L * = W is either the adjoint or the pseudo-inverse of L. In this chapter, we limit ourselves to the case where L * is the adjoint and defer the possibility of L * being the pseudo-inverse for further investigations. In summary, we propose the following steps for using the localized data-misfit function 126 defined above:

1. Find an appropriate residual function χ :

D × D → D ′ .
2. Find an appropriate linear map L : I → D ′ .

3. Find the adjoint L * : D ′ → I.

4. Determine the weight function ω from L * = W.

III.2.6.1 Example for vertical seismic profiling data

A common choice of χ : D × D → D ′ for VSP data is the L 2 -norm of data windowed around the first-arrival pick [START_REF] Pratt | Seismic waveform inversion in the frequency domain, part 2: Fault delineation in sediments using crosshole data[END_REF]:

χ f o (s, r, t), f i (s, r, t) = t=τ (s,r)-∆t t=τ (s,r)+∆t f o (s, r, t) -f i (s, r, t) 2 = f (s, r) , (127)
where τ (s, r) is the traveltime from s to r, and ∆t is half the time-window size. In this case, an object in the reduced data-space D ′ is a function of the source and receiver coordinates (s, r); therefore, we need to look for L : f (x) → f (s, r). Let us choose L to be the integration of traveltime along a ray from the source to the receiver:

L[m(x)] = x dl(s, r, x)[m(x)] = τ (s, r) ,
where dl(s, r, x) is a ray-segment centered at x along the ray, and m(x) is the slowness model at x. The adjoint L * is readily available (e.g. Christensen, 2010, chapter 4):

L f (x) , f (s, r) D ′ = s,r x dl(s, r, x)f (x) f (s, r) = x f (x) s,r dl(s, r, x)f (s, r)
= f (x), L * f (s, r) I ; (128)
< •, • > I denotes an inner product defined in I. By comparing L * as defined in equation 128 to equation 125, we can identify ω(s, r, x) = dl(s, r, x). Because the rays can always be interpolated such that all ray-segments dl(s, r, x) = 1, we conclude that w(s, r, x) = 1 if x is along the ray connecting s and r, 0 otherwise .

This implies that, for VSP, we can also rewrite the data-misfit 125 as

Φ d i (Ω) = s,r ω(s, r, Ω)χ f o (s, r, t), f i (s, r, t) , (130)
where ω(s, r, Ω) is the length of the ray-segment that actually goes through Ω.

The cartoon in Figure 77 shows the geometrical interpretation of using the data-misfit 125 along with equation 127 and equation 129. The data-misfit is projecting data residuals along rays that illuminate the region of interest Ω; this projection (of data residuals into the image-space) is equivalent to a generalized Radon transform back-projection (e.g. Toft, 1996, chapter 4).

Figure 77a shows rays emanating from sources at the free surface to receivers in a well. The green raypaths illuminate the region of interest Ω, while the black raypaths do not; only the green raypaths will actually contribute to the data-misfit 125. However, the data-misfit is still not localized enough: errors are projected along entire rays since every point x along the green rays has ω(s, r, x) = 1, according to equation 129. Therefore, for χ and ω as given by equation 127 and equation 129, respectively, it is impossible to limit the contribution of the data-misfit 125 to points x ∈ Ω in the most general case. This is a violation of condition 6. Figure 77b shows how we can localize the misfit better by adding an additional well in order to illuminate the region of interest from a different direction. In this case, points x along the green rays have ω(s, r, x) = 2 if x ∈ Ω (i.e. where the two family of raypaths cross each other), otherwise ω(s, r, x) = 1.

III.2.6.2 Example for surface seismic data

A common choice of χ : D × D → D ′ for surface data is the L p (p ∈ [1, 2]) difference data residuals [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF], i.e.

χ f o (s, r, t), f i (s, r, t) = |f o (s, r, t) -f i (s, r, t)| p . (131)
Our numerical experiments suggest that this is not the best choice of χ for appraising structural interpretations using reflection data. Consider for example the reference and interpreted models in Figure 73a and Figure 74a, respectively; reflection shot-gathers from each of these models are shown in Figure 78a and Figure 78b. In Figure 74b, we distinguish interpretation errors from picking errors. Interpretation errors typically result from mispositioning geological interfaces due to poor resolution of a seismic image, while picking errors result from mispositioning geological interfaces due to the fact that it is rarely possible to perfectly track manually a reflector on a seismic image. Looking at the model-space residuals (Figure 74b), picking errors are practically negligible, while in the data-space (Figure 78c) picking errors are as strong as interpretation errors. This is problematic because picking errors are inevitable in practice.

In Figure 78c, the picking errors identified by the top green arrow are strong because of the high amplitude of early arrivals; the picking errors identified by the bottom green arrow are important because of the velocity error introduced by the interpretation error above them. This second point motivates the use of macro-layered velocity models that are kinematically equivalent to the migration velocity model when appraising structural interpretations using surface In addition to being very sensitive to picking errors, equation 131 also suffers from cycleskipping: if a horizon's position is progressively shifted away from its true position, there is a limit beyond which a data-misfit computed using equation 131 will no longer depend on the horizon's position [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF]. This is a violation of condition 7. For instance, the data-misfit r, t)] between the two traces in Figure 79a would not change significantly if the mispositioned horizon in Figure 74a was shifted further down because the corresponding reflection-pulse no longer overlaps with the reference one. We therefore conclude that the L p (p ∈ [1, 2]) seismogram difference is not the best choice for χ. In principle, phase data residuals can overcome the aforementioned problems. Here, a phase residual refers to any function of time that measures phase-shifts between two signals. Such residual functions include dynamic-warping-based techniques [START_REF] Hale | Dynamic warping of seismic images[END_REF], and optimaltransport-based techniques [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF]. The dynamic-time-warping (DTW) phase residual δφ(s, r, t) of f o (s, r, t) and f i (s, r, t) is defined as [START_REF] Hale | Dynamic warping of seismic images[END_REF]:

φ d i (s, r) = t χ[f o (s, r, t), f i (s,
δφ(s, r, t) such that t f o s, r, t -f i s, r, t + δφ(s, r, t) 2 is minimized. ( 132
)
Figure 79d shows the data residuals of the two signals in Figure 79a using this choice of χ. We note that such a misfit is more sensitive to horizon mispositioning errors, and less sensitive to picking errors. However, such a misfit can have some undesired effects as highlighted by the red ellipse. When using equation 131 or 132, an object in the reduced data-space D ′ is a function of time as well as source and receiver coordinates, that is D ′ = D; therefore we need to look for L : f (x) → f (s, r, t). For simplicity, let us choose L to be a Kirchhoff-type modeling operator

L m(x) = x α(s, x, r)S t -τ (s, x, r) m(x) = f (s, r, t) , (133)
where f (s, r, t) is the computed trace recorded at receiver r from a shot at s; S(t) is the source wavelet; α(s, x, r) are "appropriate" migration weights as may be found, for example, in [START_REF] Cohen | Three-dimensional Born inversion with an arbitrary reference[END_REF]; τ (s, x, r) is the traveltime for a ray originating from s scattered off x and recorded at r; m is the reflectivity model. The adjoint L * is readily available:

L f (x) , f (s, r, t) D ′ =D = s,r,t x α(s, x, r)S t -τ (s, x, r) f (x) f (s, r, t) = x f (x) s,r,t α(s, x, r)S t -τ (s, x, r) f (s, r, t) = f (x), L * f (s, r, t) I . (134)
L * as defined in equation 134 differs from the standard Kirchhoff depth migration (KDM) in that it does not have a derivative along the vertical axis. One often defines [START_REF] Santos | Seismic modeling by demigration[END_REF]; we use this latter definition of KDM hereafter.

KDM f (s, r, t) = ∂ ∂z L * f (s, r, t) (Schneider, 1978), or KDM f (s, r, t) = L * ∂ ∂t f (s, r, t) (135)
By comparing L * as defined in equation 134 to equation 125, we can identify

ω(s, r, x) = t α(s, x, r)S t -τ (s, x, r) . (136)
The interpretation of this result is that the localized data-misfit 125 merely amounts to migrating data residuals for surface data. Figure 80a shows the L * projection of L 1 data residuals of waveforms computed in the model in Figure 73a and waveforms computed in the model in Figure 74a; notable contribution of picking errors (green arrows) and velocity errors (red arrow) are observed. Figure 80b shows the L * projection of L 1 data residuals of waveforms computed in the model in Figure 73a and waveforms computed in the model in Figure 76; velocity errors are dramatically reduced since the model in Figure 76 is more kinematically accurate than the model in Figure 74a. Figure 81a shows the KDM of DTW data residuals of waveforms computed in the model in Figure 73a and waveforms computed in the model in Figure 76. We observe that high-frequency noise is present and that the residual map has two poles: positive contribution above the black line, negative contribution below the black line.

Both the high-frequency noise and the bipolarity of projected residuals are a footprint of the derivative along the vertical axis in KDM, and can be removed by using the L * projection as shown in Figure 81b.

III.2.7 Application & Discussion

Appraising structural interpretations quantitatively is a challenging problem. Consider the eight different structural models in Figure 82 interpreted from the seismic image in Figure 71b by different interpreters. We would like to rank those structural models using VSP data from the most likely to the least likely. For each structural model, we built a block macro-layered velocity model, allowing us to compute VSP synthetic data for each structural model. Observed VSP data were generated in the reference model in Figure 70, with sources covering the free surface and receivers positioned in a vertical well at x = 7500 m. First, we tried to rank the structural models using a data-misfit defined as the L 1 -norm of the entire data set, for each model. Figure 83a shows that the data-space ranking obtained this way is not consistent with the ranking expected in the model-space (the model-space misfits were computed as the L 1 -norm of block macro-layered velocity models and the reference velocity model). We then tried to rank the structural models using the method proposed in this chapter. Figure 83b shows the results obtained by limiting the model-misfit and the data-misfit to the region of interest, i.e. the black box in Figure 71b, and by relying on equations 125, 127 and 129 to define the data-misfit. We observe that the scatter plot in Figure 83b satisfies condition 8, while the scatter plot in Figure 83a does not. The model-misfit axis in Figures 83a and83b show that model 1 is the best model. The data-misfit in 83b was able to indentify model 1 as the best model, whereas the data-misfit in Figure 83a identified model 3 as the best model.

In fact, the data-misfit in Figure 83b was able to successfully rank all the models except model 7. The present chapter focuses on misfit functions with the aim of determining the possible reasons as to why the straightforward approach implemented in Figure 83a fails. First, we argue that one should use localized data-misfits Φ d (Ω). It follows that conditions 6 and 7 are necessary conditions for being able to rank structural models using Φ d (Ω), if we expect data-space ranking to be consistent with model-space ranking. These conditions are expected to be valid if the available candidate structural models are "close enough" to the (unknown) reference model.

Condition 6 mandates that a perturbation of the data-misfit value Φ d i (Ω) in the region Ω of the i th model should only occur if that model is modified in that specific region Ω. Condition 7 mandates that if the i th model is modified in the region Ω, the data misfit Φ d i (Ω) corresponding to that region should change as well.

For VSP data, condition 7 can readily be satisfied by computing Φ d using only source-receiver pairs that illuminate Ω (e.g. using equation 129) and windowing data, along the time axis, around the first arrival (e.g. using equation 127). The difficulty then comes from satisfying condition 6, which is impossible in the most general case as any perturbation of the structural model at any point along any raypath through Ω will lead to a perturbation of Φ d , not just points in Ω. This difficulty can be alleviated, in principle, by adding data that illuminate Ω from different directions as illustrated in Figure 77.

For surface seismic data, using L p (p ∈ [1, 2]) data residuals, we show that velocity errors introduced in macro-layered velocity models can lead to a perturbation of Φ d (Ω) for a region Ω below the velocity error even in the absence of any structural interpretation errors in Ω (e.g. Figure 80a), thereby violating condition 6. We also argue that if a reflector in the structural model is shifted away from its true position until its reflection-pulse no longer overlaps with true reflection-pulse along the time axis, shifting the reflector further in the model would not necessarily affect the data misfit, thereby violating condition 7. However, it is possible, in principle, to approximate condition 7 by using more appropriate data residuals like phase-shift residuals. The challenge then becomes how to compute phase-shifts for complex data sets; this is a subject of ongoing investigations.

The proposed method assumes that the region of interest is well illuminated. We can therefore expect the method to perform poorly for more challenging seismic imaging targets, such as subsalt areas for example. Moreover, appraising structural models using surface seismic data requires a kinematically accurate migration velocity model, which is not easy to obtain in practice. Challenging our approach with a realistic migration velocity model obtained by migration velocity analysis is a subject of ongoing investigations.

III.2.8 Conclusion

In this chapter, we propose a theoretical formulation and some general solutions to the problem of appraising structural interpretations using seismic data. Assuming that the different structural models are close to the unknown true model, we propose a set of conditions imposed on data-space misfit functions needed for a reliable model appraisal. We argue that misfit functions should be able to localize interpretation errors in the image-space. This localization of errors is achieved by back-projection of data-residuals into the image-space. It follows that, because it is not possible to truly localize errors using VSP data, one cannot predict which interpretations are more probable than others using VSP data in the most general case. However, it is expected that VSP data can always be used to statistically rank structural interpretations if the data illuminates the target from different directions. As for appraising structural models using surface seismic data, we expect a better localization of errors, compared to VSP, and therefore a better chance for ranking structural interpretations. The challenge when using surface seismic data is to define appropriate data residual functions. We argue that phase-shift residual functions, such as dynamic-warping and optimal-transport residuals, are good candidates; our current work involves finding robust ways to compute such phase-shift residuals for data acquired in complex geological settings.

III.3.2 Introduction

Structural models built from interpretation of seismic images inherit uncertainties propagated from interpretation uncertainties. As a result, it is often possible to propose multiple structural models from a single seismic image. In this chapter we consider the problem of determining which structural models, from a set of structural models, are more likely than others using seismic data. This problem is referred to as appraising structural models hereafter. A detailed introduction to the problem of appraising structural models using seismic data is given in chapter III.2. In chapter III.2, an application of appraising structural models using vertical seismic profiling data is presented. Here, we discuss the challenges of appraising structural models using surface seismic data. The main idea of the proposed method is to simulate synthetic data for each candidate structural model available, and then compare the resulting synthetic data to observed data. Candidate structural models are assumed to result from structural interpretations of a single seismic image. As discussed in chapter III.2, the datamisfit Φ i (Ω) of the i th structural model in the region of interest Ω can be expressed as

Φ i (Ω) = x∈Ω P R(d i , d o ) = x∈Ω E(x), (137)
where d o is the observed data, d i is the synthetic data generated from the i th structural model. The data residual operator R maps the data-plane D × D to a reduced-data-space D ′ , and the weighting operator P maps the reduced-data-space D ′ to the image-domain I. The interpretation error map E defined in the data misfit 137 is therefore computed in two steps: first, the residual operator R computes data residuals; second, the weight operator P projects the data residuals in the image-domain. We will investigate only residual operators R for which D ′ = D (i.e. an object in D ′ is a function of source and receiver coordinates and time) so that one may consider P = F † , where F † is either the adjoint (e.g. section III.3.3) or the inverse (e.g. section III.3.4) of any appropriate operator

F : I → D. (138)
In section III.3.3, we consider the case where R is the dynamic image warping (DIW) operator [START_REF] Hale | Dynamic warping of seismic images[END_REF] that computes phase-shift residuals of seismic data, and P is a linear imaging operator. In this section, we study only candidate interpreted structural models that have the same number of geological layers as the reference (observed) structural model. This is an ideal case where the assumptions of DIW, and therefore its validity within this framework, can be expected to hold. In practice, however, we can expect candidate models to have a different number of layers from the reference structural model, especially if the reference model has some fine layered structures (see for example chapter III.1); in this case, the assumptions of DIW do not hold and its validity within this framework is questionable. In section III.3.4, we consider the case where R is the L1-difference and P is a waveform inversion operator. This approach does not constrain the number of layers in candidate structural models; however, it inherits all the computational challenges of waveform inversion. We will only study macrolayered models in this chapter, that is, models without any fine layer structures, even though this is not necessary for the approach presented in section III.3.4.

III.3.3 Appraising models with the same number of layers

Let us consider the case where F in equation 138 is the Born modeling operator (see for example section I.1.1)

L(m)P i (x, t; s) = 1 m 2 (x) ∂ 2 ∂t 2 -∆ P i (x, t; s) = m i (x)P (x, t; s), (139)
where P i is the output synthetic wavefield, m i is the input model, P is the solution to

L(m)P (x, t; s) = δ(x -s)w(t), (140)
m is the smooth background velocity model (assumed to be kinematically accurate), and w is the source wavelet. Synthetic data d i (t, s, r) are obtained by sampling the synthetic wavefield

P i (x, t; s) = d i (x, s, t) = x ′ ,t ′ G(x, t; x ′ , t ′ )P (x ′ , t ′ ; s)m i (x ′ ) = Bm i (x) (141)
at receivers; that is

d i (t, s, r) = x δ(x -r)d i (x, s, t) = x δ(x -r) x ′ ,t ′ G(x, t; x ′ , t ′ )P (x ′ , t ′ ; s)m i (x ′ ) = x,t ′ G(x, t -t ′ ; r)P (x, t ′ ; s)m i (x) = Fm i (x), (142)
(see the appendices in section I.1.5 for more details) where G(x, t; r) is the solution to

L(m)G(x, t; r) = δ(x -r)δ(t).
For F as defined in equation 142, one may equate its adjoint to the weight operator P in equation 137. For example, let f (s, r, t) denote the DIW data residual

f (s, r, t) = R d o (s, r, t), d i (s, r, t) = argmin l(s,r,t) s,r,t d o (s, r, t -l) -d i (s, r, t) 2 , (143)
then the adjoint of equation 142 leads to P : f (s, r, t) → f (x) being the reverse time migration (RTM) operator

P[f (s, r, t)] = r,s,t,t ′ G(x, t -t ′ ; r)P (x, t ′ ; s)[f (s, r, t)], (144)
or more customarily P[f (s, r, t)] = s,t R(x, t; r)P (x, t; s)

L * (m)R(x, t; r) = f (s, r, t) , (145)
(see the appendices in section I.1.5 for more details) where the asterisk in L * (m) denotes adjoincy; that is, the equation is solved backwards in time (e.g. Fichtner et al., 2006).

In practice, the straightforward application of equation 145 to project DIW data residuals in the image space does not work well. Consider for example the illustrative reference model in Figure 84a, the corresponding candidate synthetic model is shown in Figure 84b; the difference between the two models is shown in Figure 84c. Figure 85a shows an "observed" shot gather computed in the reference model (Figure 84a), Figure 85b shows a synthetic shot gather computed in the synthetic model (Figure 84b), and Figure 85c shows the difference between the observed and synthetic data. Ideally, given a migration velocity model (Figure 86a), we would like to construct an interpretation error map E(x) similar to Figure 84c using the the observed data (Figure 85a) and synthetic data (Figure 85b). Figure 86b shows the DIW data residuals obtained from the observed and synthetic data in Figure 85 and Figure 86c shows the result of applying equation 145 to the DIW data residuals using the migration model in Figure 86a. The resulting error map is noisy and does not reveal the interpretation error observed in Figure 84c. An improved error map, shown in Figure 86d, is obtained by replacing G and P in equation 144 by their ray-based approximations (i.e. Kirchhoff depth migration). For this simple example, the error map obtained by Kirchhoff migration (KDM) reveals the interpretation error in Figure 84c. For more complex models however, the KDM-based error map does not work well either. Let us first consider the slightly more complex model in Figure 87. The reference and the candidate model from interpretation are shown in Figures 87a and87b; the difference between the two models is shown in Figure 87c. The "observed" and synthetic data are shown in Figures 88a and88b; the difference between the observed and synthetic data is shown in Figure 88c.

The DIW data residuals computed from the observed (Figure 88a) and synthetic (Figure 88b) are shown in Figure 89a. The KDM projection of DIW data residual in the image-space, using the migration model in Figure 89b, is shown in Figure 89c. The interpretation error in blue about 1000m deep in Figure 87c is detected, as clearly seen in Figure 89c. However, the rest of the interpretation errors observed beyond 1000m are not visible in Figure 87c. Now let us consider the more complex reference model in Figure 90a; the corresponding interpreted model is shown in Figure 90b, and the difference between the two is shown in Figure 90c.

The "observed" and synthetic data are shown in Figures 91a and91b; the difference between the observed and the reference data is shown in Figure 91c. The DIW data residuals computed from the observed (Figure 91a) and synthetic (Figure 91b) are shown in Figure 92a. The KDM projection of DIW data residual in the image-space, using the migration model in Figure 92b, is shown in Figure 92c. The resulting error map was clearly unable to detect the interpretation errors observed in Figure 90c. While it is still unclear why equation 145 does not work as expected, there are at least two issues that we think may be leading to the unsatisfactory results observed:

1. The shifts computed by DIW may not be accurate; this seems to be especially true towards the end of seismic traces dominated by multiples. Even in the ideal case where the reference and synthetic model have the same number of layers, as investigated in this section, the synthetic data will generally not satisfy the DIW assumption requiring the synthetic data to be a warped version of the observed data. Furthermore, DIW, as presented by [START_REF] Hale | Dynamic warping of seismic images[END_REF], has a couple of free parameters. DIW data-residuals have shown to strongly depend on these parameters and it has proven challenging to determine the optimal values for these parameters.

2. The output data residuals from equation 143 are usually very low frequency data, as can be observed in Figures 86b,89a, and 92a. Both our intuition and numerical experiments suggest that it is unlikely to reconstruct interpretation errors using equation 145 because interpretation errors usually contain both high and low frequency as observed in Figures 84c,87c, and 90c.

The first issue raised above, that of computing accurate phase-shift residuals, remains an open problem, at least within the scope of this work.

We now discuss the second issue raised above: the low frequency content of DIW phase-shift residuals. An alternative way to think of P in equation 144 is as the first iteration of the least-square inverse of equation 142 (i.e. least-squares RTM, see section I.1.2 for details). The advantage of this approach is that it automatically accounts for the low frequency nature of residuals obtained using 143. Consider the optimization problem

             m(x) = argmin m(x) s,r,t f (s, r, t, m) 2 = s,r,x,t δ(x -r)f (s, x, t, m) 2 subject to f (s, r, t, m) = argmin l(s,r,t,m) s,r,x,t δ(x -r) d o (s, x, t -l) -d i (s, x, t, m) 2 d i (s, x, t, m) = Bm i , (146)
resulting from equations 141 and 143. Ma and Hale (2013) note from equation 143 that

∂ ∂l R| f = 0 = -2 s,r,t d o (s, r, τ ) -d i (s, r, t) ∂ ∂τ d o (s, r, τ ), with τ = t -f, (147)
(see the appendix in section III.3.7 for details) transforming problem 146 to

           m(x) = argmin m(x) = s,r,x,t δ(x -r)f (s, x, t, m) 2 subject to 0 = s,r,x,t δ(x -r) d o (s, x, τ ) -d i (s, x, t, m) ∂ ∂τ d o (s, x, τ ), with τ = t -f (s, x, t, m) d i (s, x, t, m) = Bm i . (148)
The Lagrangian of problem 148 is given by

L = δ(x -r), f (s, x, t, m) 2 s,r,x,t + λ 1 δ(x -r), d o (s, x, τ ) -d i (s, x, t, m) ∂ ∂τ d o (s, x, τ ) s,r,x,t + λ 2 , d i (s, x, t, m) -Bm i s,x,t . ( 149
)
The derivative of the Lagrangian with respect to the model is given by (see the appendix in section III.3.7 for details)

∂ ∂m L = δ(x -r) ∂ ∂m f, 2f -λ 1 ḋo (s, x, τ ) 2 + d o (s, x, τ ) -d i (s, x, t, m) do (s, x, τ ) s,r,x,t + ∂ ∂m d i (s, x, t, m), λ 2 - r λ 1 δ(x -r) ḋo (s, x, τ ) s,x,t - B * λ 2 , 1 x , (150)
or more conveniently

         λ 1 = 2 f (s,x,t,m) ḋo(s,x,τ ) 2 + do(s,x,τ )-d i (s,x,t,m) do(s,x,τ ) λ 2 = r λ 1 δ(x -r) ḋo (s, x, τ ) ∂ ∂m L = -B * λ 2 , 1 x , ( 151
)
where B * is the adjoint of the Born modeling operator defined in equation 141 (i.e. B * is the RTM operator). It follows from equation 151 that the first iteration of least-squares RTM, using the DIW equation 143, leads to

       P[f (s, r, t)] = s,t R(x, t; r)P (x, t; s) (a) L * (m)R(x, t; r) = f (s,r,t,m) ḋo(s,r,τ ) ḋo(s,r,τ ) 2 + do(s,r,τ )-d i (s,r,t,m) do(s,r,τ ) (b) τ = t -f (s, r, t) (c)
.

The term on the right-hand side of equation 152b will be referred to as the adjoint-source.

According to equation 152b, the projection operator 152a overcomes the low-frequency issue by scaling the low-frequency DIW phase-shifts by the time derivative of observed traces. However, the error map obtained by the projection operator 152a will still be unreliable in general because it still strongly depends on our ability to compute accurate DIW phase-shifts. Consider for example adjoint-source in Figure 93a computed from the traces in Figure 85 and the phase-shifts in Figure 86b; the RTM projection, shown in Figure 93b, was able to detect the interpretation error in Figure 84c. For comparison, Figure 93c shows the RTM projection of the L 1 difference in Figure 85c. The L 1 data-difference error map misleadingly amplifies the picking error in Figure 86c around the depth of 2000m. Such shortcomings of L p (p ∈ [1, 2]) differences were among the motivations for investigating phase-shift residuals (see chapter III.2). These theoretical advantages of phase-shift residuals over data-difference residuals become less convincing in practice as models get more complex. For example, given the adjoint-source in Figure 94a computed from the traces in Figure 88 and the phase-shifts in Figure 89a, the superiority of the adjoint-source error map in Figure 94b over the L 1 difference error map in Figure 94c in detecting the interpretation errors in Figure 87c is less obvious. In fact, for the more complex example where the adjoint-source in Figure 95a is computed from the traces in Figure 91 and the phase-shifts in Figure 92a, the adjoint-source error map in Figure 95b appears to give a poorer estimate of the interpretation errors in Figure 90c than the L 1 data-difference error map in Figure 95c. We know from theoretical arguments and from numerical experiments that L p (p ∈ [1, 2]) datadifference residuals are not suited for appraising structural interpretations using a linear projection operator P (see chapter III.2). As we have just shown in this section, while DIW phase-shift residuals offer some theoretical advantages, they appear not to be suited for appraising structural interpretations using a linear projection operator P as well in practice. We investigate the use of a nonlinear projection operator P in the following section.

III.3.4 Appraising models with different number of layers

Consider the experiment shown in Figure 96. The velocity model in Figure 96a, obtained by adding the migration velocity model in Figure 86a to the error map in Figure 84c, was used to generate the synthetic data in Figure 96b. A waveform inversion (see section I.1.3) of this data set results in the velocity error map in Figure 96c, whose absolute value is presented in Figure 96d. In some sense, the error map in Figure 96d represents the best estimate of the error map we can expect to retrieve from data residuals using waveform inversion because it was obtained from the true error map (Figure 84c). It follows from this experiment that if we can use the observed data (Figure 85a) and synthetic data (e.g. Figure 85b) to compute data residuals that approximate the data in Figure 96b "well enough", then these data residuals can be used to estimate the error map by waveform inversion. Following the discussion in the previous paragraph, let

a o (s, r, t) = R(d i , d o )
denote data residuals computed using the observed data d o and synthetic data d i obtained from the i th structural model, and let us assume that a o (s, r, t) is a "good" approximation of data obtained by means of the experiment in Figure 96. Then the interpretation error map corresponding to the i th structural model can be estimated by solving the following waveform inversion problem

         δm i = argmin b i (x) s,r,t a o (s, r, t) -a i (s, r, t, b i ) 2 = argmin b i (x) s,r,x,t δ(x -r) a o (s, x, t) -a i (s, x, t, b i ) 2 subject to L(m + b i )a i (x, t, b i ; s) = δ(x -s)w(t) . ( 153
)
The Lagrangian of problem 153 is given by

L = δ(x -r), a o (s, x, t) -a i (s, x, t, b i ) 2 s,r,x,t + λ, L(m + b i )a i (x, t, b i ; s) -δ(x -s)w(t) s,x,t . (154)
The derivative of the Lagrangian with respect to the model b(x) is given by

∂ ∂b L = ∂ ∂b a i (s, x, t, b i ), L * (m + b i )λ -2 r δ(x -r) a o (s, x, t) -a i (s, x, t, b i ) s,x,t + ∂ ∂t λ, 2 (m + b i ) 3 ∂ ∂t a i (x, t, b i ; s) s,x,t , (155)
(see the appendix in section III.3.7 for details), or more conveniently

     L * (m + b i )λ = 2 r δ(x -r) a o (s, x, t) -a i (s, x, t, b i ) ∂ ∂b L = ∂ ∂t λ, 2 (m+b i ) 3 ∂ ∂t a i (x, t, b i ; s) s,x,t . (156)
The derivative of the Lagrangian, equation 156, can be used to solve problem 153 using an iterative algorithm of the form

           L(m + b i )a i (x, t, b i ; s) = δ(x -s)w(t) (a) L * (m + b i )R(x, t, b i ; r) = r δ(x -r) a o (s, x, t) -a i (s, x, t, b i ) (b) ∆b i = 1 (m+b i ) 3 s,t ∂ ∂t R(x, t, b i ; r) ∂ ∂t a i (x, t, b i ; s) (c) b n+1 i = b n i + ∆b n i (d) . ( 157
)
Let us consider the case where the data residuals are obtained using a simple data difference

a o (s, r, t) = d i -d o . (158)
Figure 97 shows an estimate of the error map in Figure 84c obtained by waveform inversion when the synthetic data d i in equation 158 is generated using a block macro-layered velocity model (section III.2.5.1). We note that the error map in Figure 97b differs from the error map in Figure 96d at least in two aspects:

1. The error map in Figure 97b shows an additional interpretation error around the depth of 2000m.

2. The error map in Figure 97b shows additional random oscillatory noise scattered throughout the error map.

The first issue comes from the fact that block macro-layered models are not kinematically accurate; this issue can be alleviated by using macro-layered models that are more kinematically accurate. For example, Figure 98 shows the error map obtained by waveform inversion when the synthetic data d i in equation 158 is generated using a ref lectivity macro-layered velocity model (section III.2.5.3). Clearly the error map in Figure 98b is a better estimate of the error map 96d than the error map in Figure 97b. The second issue raised above (random oscillatory noise) probably comes from the fact that data residuals obtained by the data difference 158 differ from data obtained by means of the experiment in Figure 96. One way to reduce this oscillatory noise is to use a regularization operator that penalizes oscillations such as the total variation (TV) regularization (e.g. [START_REF] Anagaw | Edge-preserving seismic imaging using the total variation method[END_REF][START_REF] Lopez | Speed up and regularization techniques for seismic full waveform inversion[END_REF][START_REF] Esser | Total-variation regularization strategies in full-waveform inversion[END_REF]. The TV regularized version of problem 153 is

   δm i = argmin b i (x) s,r,x,t δ(x -r) a o (s, x, t) -a i (s, x, t, b i ) 2 + β x ( ∂ ∂ x1 b i ) 2 + ( ∂ ∂ x2 b i ) 2 subject to L(m + b i )a i (x, t, b i ; s) = δ(x -s)w(t)
.

(159) where x = (x 1 , x 2 ), and β ∈ R is a user-defined regularization weight. It can be shown that problem 159 can be solved iteratively using an algorithm of the form

             L(m + b i )a i (x, t, b i ; s) = δ(x -s)w(t) (a) L * (m + b i )R(x, t, b i ; r) = r δ(x -r) a o (s, x, t) -a i (s, x, t, b i ) (b) ∆b i = 1 (m+b i ) 3 s,t ∂ ∂t R(x, t, b i ; r) ∂ ∂t a i (x, t, b i ; s) + β∇ • 1 ( ∂ ∂ x1 b i ) 2 +( ∂ ∂ x2 b i ) 2 ∇b i (c) b n+1 i = b n i + ∆b n i (d)
.

The algorithm in equation 160 is derived from problem 159 in the same way the algorithm in equation 157 was derived from problem 153; the derivation of the extra term in equation 160c is given in the appendices (section III.3.7). The TV regularized waveform inversion estimate of the error map in Figure 84c is shown in Figure 99, the random oscillatory noise has been considerably reduced. We now test the waveform inversion method on the slightly more complex model shown in Figure 87. Figure 100a shows the error map estimated by means of the experiment in Figure 96; this error map is the best we can hope to get without knowledge of the reference model by waveform inversion as described in this section. Figure 100b shows the error map estimated from synthetic data generated in a block macro-layered velocity model, while Figure 100c shows the error map estimated from synthetic data generated in a reflectivity macro-layered velocity model. The error map from the block macro-layered model (Figure 100b) visually looks closer to the target error map (Figure 100a) than the error map from the reflectivity macro-layered model (100c). However, because we know that the block macro-layered model may introduce velocity errors as discussed previously, we will always prefer the error map from the reflectivity macro-layered model and try to improve it by regularization. The TV regularized version of the error map in Figure 100c is shown in Figure 100d. We find that the error map in Figure 100d is a reasonable estimate of the target error map in Figure 100a, itself a reasonable estimate of the true error map in Figure 87c.

III.3.5 Conclusion & Discussion

In this chapter, we study the practical aspects of appraising structural interpretations using surface seismic data. The global workflow we use for appraising structural interpretations is to generate synthetic data for each interpretation, and then use the synthetic data to compute a misfit value for each interpretation. We use misfit functions that are localized in space (i.e. a misfit function that is a function of space). The localized misfit functions are composed of two main components: (1) a residual operator that computes the difference between synthetic and observed data, (2) and a projection operator that projects the data residuals into the image-domain. We present two strategies for computing the localized misfit functions:

• In the first strategy (section III.3.3), the residual operator does most of the work in the misfit function by computing phase-shifts between observed and synthetic data; these phase shifts are intended to be proportional to horizon shifts between the reference and interpreted structural models. A linear projection operator is then used to project the phase-shifts to the image-domain. The challenge in this first strategy is to propose ap-propriate residual operators that can compute phase-shifts accurately. We conclude that phase-shift residual operators based on dynamic warping are not appropriate mainly because they rely on the assumption that observed data are a warped version of synthetic data (e.g. [START_REF] Hale | Dynamic warping of seismic images[END_REF]. This assumption implies that synthetic structural models (models from seismic interpretation) should have the same number of layers as the true structural model, which is usually not the case. Furthermore, the assumption is usually also not valid for faulted complex structural models even when they have the same number of layers as the reference model. A natural perspective for this first strategy is to test other phaseshift residual operators, particularly operators based on optimal transport (e.g. [START_REF] Engquist | Application of the Wasserstein metric to seismic signals[END_REF][START_REF] Métivier | Measuring the misfit between seismograms using an optimal transport distance: Application to full waveform inversion[END_REF].

• In the second strategy (section III.3.4), the residual operator is based on a simple L1difference; in this case, the projection operator does most of the work in the misfit function by projecting the data residual in the image domain using an iterative process. Here, the projection operator is based on waveform inversion. This strategy therefore inherits all the computational challenges already familiar from the waveform inversion literature. The advantage of this strategy is that it does not require synthetic structural models to have the same number of layers as the reference model. Current research on this second strategy involves finding appropriate regularization schemes to stabilize the waveform inversion of data residuals. A promising perspective to investigate is the regularization scheme presented by [START_REF] Esser | Total-variation regularization strategies in full-waveform inversion[END_REF]; their regularization scheme seems to be able to recover low wavenumber velocity components from inversion. Such an ability would certainly be valuable for our method.

While we concluded that phase-shift residuals based on dynamic warping were not appropriate for the first strategy, we expect them to be useful in the second strategy. In particular, consider the problem

               a o (s, r, t) = R 1 d i (s, r, t), d o (s, r, t) (a) δm i = argmin b i (x) s,r,t R 2 a i (s, r, t), a o (s, r, t) 2 (b) = argmin b i (x) s,r,x,t δ(x -r)R 2 a i (s, x, t), a o (s, x, t) 2 (c) subject to L(m + b i )a i (x, t, b i ; s) = δ(x -s)w(t) (d)
;

(161) the method presented in section III.3.4 is a specific case of problem 161 where both R 1 and R 2 are L1-difference data residual operators. A promising perspective is to solve problem 161 where R 1 is an L1-difference data residual operator and R 2 is a phase-shift data residual operator (based on dynamic warping for example). It should be noted in problem 161 that R 1 is the residual operator R in equation 137 while R 2 is part of the projection operator P in equation 137.

III.3.7 Appendix

Derivation of equation 147

∂ ∂l R| f = s,r,t ∂ ∂f d o (s, r, t -f ) -d i (s, r, t) 2 = s,r,t ∂ ∂f d o (s, r, g) -d i (s, r, t) 2 , with g(f ) = t -f = s,r,t ∂ ∂f h(g) 2 , with h(g) = d o (s, r, g) -d i (s, r, t) = s,r,t ∂ ∂h h 2 ∂ ∂g h(g) ∂ ∂f g(f ) = -2 s,r,t d o (s, r, t) -d i (s, r, g) ∂ ∂g d o (s, r, g)

Derivation of equation 150

The derivative of the first term in equation 149 with respect to the model m is given by

∂ ∂m δ(x -r), f (s, x, t, m) 2 s,r,x,t = δ(x -r), 2f (s, x, t, m) ∂ ∂m f (s, x, t, m) s,r,x,t . The derivative ∂ ∂m d o (s, x, τ ) -d i (s, x, t, m) ḋo (s, x, τ ) = ḋo (s, x, τ ) ∂ ∂m d o (s, x, τ ) -d i (s, x, t, m) + d o (s, x, τ ) -d i (s, x, t, m) ∂ ∂m ḋo (s, x, τ ) = ḋo (s, x, τ ) ḋo (s, x, τ ) ∂ ∂m τ - ∂ ∂m d i (s, x, t, m) + d o (s, x, τ ) -d i (s, x, t, m) do (s, x, τ ) ∂ ∂m τ = ḋo (s, x, τ ) 2 + d o (s, x, τ ) -d i (s, x, t, m) do (s, x, τ ) ∂ ∂m τ -ḋo (s, x, τ ) ∂ ∂m d i (s, x, t, m) = -ḋo (s, x, τ ) 2 + d o (s, x, τ ) -d i (s, x, t, m) do (s, x, τ ) ∂ ∂m f -ḋo (s, x, τ ) ∂ ∂m d i (s, x, t, m)
is used to compute the derivative of the second term in equation 149. The derivative of the last term in equation 149 is given by

∂ ∂m λ 2 , d i (s, x, t, m) -Bm i (x) s,x,t = λ 2 , ∂ ∂m d i (s, x, t, m) s,x,t -λ 2 , B ∂ ∂m m i (x) s,x,t = λ 2 , ∂ ∂m d i (s, x, t, m) s,x,t -B * λ 2 , 1 x ; that f (s, x, t), Bf (x) s,x,t = B * f (s, x, t), f (x)
x follows from equation 19 and equation 20:

B is the Born modeling operator, B * is the reverse time migration operator, and B * f (s, x, t) = s,t Bf (s, x, t) .

Derivation of equation 155

The derivative of the first term in equation 154 with respect to the model parameter b is given by

∂ ∂b δ(x -r), a o (s, x, t) -a i (s, x, t, b i ) 2 s,r,x,t = δ(x -r), -2 a o (s, x, t) -a i (s, x, t, b i ) ∂ ∂b a i (s, x, t, b i ) s,r,x,t
.

The derivative of the second term in equation 154 is given by

∂ ∂b λ, L(m + b i )a i (x, t, b i ; s) -δ(x -s)w(t) s,x,t = λ, ∂ ∂b L(m + b i ) a i (x, t, b i ; s) s,x,t + λ, L(m + b i ) ∂ ∂b a i (x, t, b i ; s) s,x,t = λ, - 2 (m + b i ) 3 ∂ 2 ∂t 2 a i (x, t, b i ; s) s,x,t + L * (m + b i )λ, ∂ ∂b a i (x, t, b i ; s) s,x,t = ∂ ∂t λ, 2 (m + b i ) 3 ∂ ∂t a i (x, t, b i ; s) s,x,t + L * (m + b i )λ, ∂ ∂b a i (x, t, b i ; s) s,x,t
, where L * is the adjoint of L.

Derivation of equation 160c

The first term in equation 160c is already familiar from equation 157c. The second term in equation 160c comes from the derivative of the regularization term in problem 159, in particular

∂ ∂b 1, ( ∂ ∂ x1 b) 2 + ( ∂ ∂ x2 b) 2 x = 1, ∂ ∂b ( ∂ ∂ x1 b) 2 + ( ∂ ∂ x2 b) 2 x = 1, ∂ ∂b f x , with f = ( ∂ ∂ x1 b) 2 + ( ∂ ∂ x2 b) 2 = 1, ∂ ∂f f ∂ ∂b f x = 1 2 √ f , ∂ ∂b f x = 1 2 √ f , ∂ ∂b ( ∂ ∂ x1 b) 2 + ( ∂ ∂ x2 b) 2 x = 1 2 √ f , ∂ ∂b f 2 1 + ∂ ∂b f 2 2 x , with f 1 = ∂ ∂ x1 b and f 2 = ∂ ∂ x2 b = 1 2 √ f , ∂ ∂f 1 f 2 1 ∂ ∂b f 1 + ∂ ∂f 2 f 2 2 ∂ ∂b f 2 x = 1 √ f , f 1 ∂ ∂b f 1 + f 2 ∂ ∂b f 2 x = 1 √ f , ∂ ∂ x1 b ∂ ∂b ∂ ∂ x1 b + ∂ ∂ x2 b ∂ ∂b ∂ ∂ x2 b x = 1 √ f , ∇b • ∂ ∂b ∇b x = 1 √ f , ∇b • ∇ ∂ ∂b b x = 1, 1 √ f ∇b • ∇ 1 x = -∇ • 1 √ f ∇b , 1 x = -∇ • 1 ( ∂ ∂ x1 b) 2 + ( ∂ ∂ x2 b) 2 ∇b , 1 x
The last line follows from the divergence theorem and neglecting boundary contributions.

General Conclusion

Subsurface structural models find applications in natural resource estimation, risk management, and numerical simulations; it is therefore important that structural models represent the geometry of subsurface geological objects accurately. In this manuscript, we study the problem of reducing structural uncertainties using seismic data. In particular, we study the problem of using seismic data to determine which structural models are more likely than others in a set of geologically plausible structural models. We refer to this problem as appraising structural models using seismic data. Appraising structural models using seismic data is an integrative workflow that involves seismic imaging and structural modeling. Consequently, this manuscript also investigates some problems in seismic imaging and structural modeling.

The first Part of the manuscript is devoted to seismic imaging. Seismic imaging based on waveform inversion produces optimal quantitative seismic images of the subsurface by minimizing the difference between observed and synthetic data. Waveform inversion is however computationally intensive. In chapter I.2, I propose to use reverse-time migration (RTM) as a preconditioner for waveform inversion. Numerical experiments show that the proposed preconditioner accelerates both linearized waveform inversion (least squares reverse time migration) and nonlinear waveform inversion (full waveform inversion) by at least an order of magnitude. I justify the positive numerical performance of the proposed preconditioner by showing algebraically that a low pass filter of the RTM image can approximate the diagonal elements of the Hessian matrix of the objective function under appropriate assumptions. However, I am still unable to propose a physical meaning of the low pass filtering and how it relates the RTM image to the elements of the diagonal of the Hessian matrix; this is still a subject of ongoing investigation. All seismic imaging methods require an accurate imaging velocity model. In chapter I.3, we propose a generalized extended Kirchhoff imaging operator for velocity modeling; the operator is generalized in the sense that it describes multiple data-domain extensions (e.g. shot, offset, and angle extensions) and image-domain extensions (e.g. time-lag and space-lag extensions) simultaneously. The advantages of the proposed generalized extended operator are twofold: firstly, it allows a unified implementation for multiple extensions (i.e. a single implementation that is valid for multiple extensions); secondly, the operator leads to a unified gradient-based migration velocity analysis (MVA) scheme. We confirm the ability of the proposed generalized extended operator to capture image distortion caused by inaccurate velocity by applying it to a ray-based MVA experiment. We then use the proposed operator to derive a diferential semblance gradient-based MVA that is valid simultaneously for all data-domain extensions; the practical implementation of this unified gradient-based MVA has been postponed for further investigation.

The second Part of the manuscript is devoted to structural modeling, particularly implicit structural interpolation. In chapter II.1, we introduce Finite Difference Structural Implicit Modeling (FDSIM). The advantages of FDSIM are twofold: firstly, it is relatively easy to implement and to optimize since it is based on finite differences on regular grids; secondly, because it handles discontinuities by rasterization, FDSIM has shown to easily handle very complex fault-networks. The main disadvantage of the method is that it may require a very fine resolution depending on the complexity of the fault-network, sometimes leading to memory limits. We also propose new regularization operators in chapter II.1; the particularity of these operators is that they do not need to be implemented on boundary nodes, a property which is very appealing in implicit modeling where boundary conditions are usually unknown. In chapter II.2, we introduce Finite Element Structural Implicit Modeling (FESIM). FESIM is based on a finite element implementation of the regularization operators proposed in chapter II.1. We show that the conventional finite element familiar for solving boundary value problems has to be slightly modified for implicit modeling where boundary conditions are usually unknown. One advantage of FESIM over FDSIM is that it requires less memory for complexly faulted structural models. The disadvantage of FESIM over FDSIM is that it can be quite difficult to model structural models with complex fault-networks due to meshing challenges. Current research in FESIM includes finding more automated strategies to generate a mesh for structural models with complex fault-networks. In chapter II.3, we propose two new strategies for handling large thickness variations in implicit structural modeling: the first strategy constrains the output implicit function (stratigraphic field) to follow the orientation of geological layers estimated from input data, while the second strategy uses an anisotropic regularization operator that is consistent with the orientation of geological layers. These strategies are still under investigation and they show promising results.

The third Part of the manuscript is devoted to appraising structural models/interpretations using seismic data. In chapter III.1, we discuss about interpretation uncertainties and show how a single seismic image may support multiple geologically consistent structural models. In chapter III.2, we introduce and formalize the problem of appraising structural interpretations using seismic data. We propose to solve the problem by generating synthetic data for each structural interpretation and then compute misfit values for each interpretation; this allows us to rank the different structural interpretations. The main challenge of appraising structural models using seismic data is to propose appropriate data misfit functions. We derive a set of conditions that have to be satisfied by the data misfit function for a successful appraisal of structural models. We argue that since it is not possible to satisfy these conditions using vertical seismic profile (VSP) data, it is not possible to appraise structural interpretation using VSP data in the most general case. However, we show that it is always possible to get a better approximation of those conditions using VSP data that illuminate the target region from different directions. We present an application of appraising structural interpretations using VSP data to a realistic synthetic case. The conditions imposed on the data misfit function can in principle be satisfied for surface seismic data. In practice however, it remains a challenge to propose and compute data misfit functions that satisfy those conditions. In chapter III.3, we focus on the practical aspects of appraising structural interpretations using surface seismic data. We propose a general data misfit formula that is made of two main components: (1) a residual operator that computes data residuals, that is, the difference between observed and synthetic data, and (2) a projection operator that projects the data residuals to the imagedomain. This misfit function is therefore localized in space, as it outputs data misfit values in the image-domain. We present two strategies for implementing the localized data misfit function:

• In the first strategy, the residual operator does most of the work in the misfit function by computing phase shifts between synthetic and observed data; phase shifts are, in principle, correlated with horizon shifts between the true structural model and the interpreted structural model. A linear operator is then used to project the phase-shift residuals into the image-domain to give an indication of horizon shifts in the subsurface. We conclude that phase-shift based on dynamic-warping are not appropriate for appraising structural models, when using a linear projection operator, mainly because they require observed data to be a warped version of synthetic data, a requirement usually not satisfied for this particular problem.

• In the second strategy, the residual operator does most of the work in the misfit function by iteratively projecting to the image-domain data residuals based on L1-differences. In this case, the projection operator is based on nonlinear waveform inversion and therefore inherits all the computational challenges familiar from the full waveform inversion literature. Waveform inversion of data residuals has shown to be quite unstable and often requires regularization.

We believe that future research on appraising structural models using surface seismic data should include investigating other phase-shift residual operators (e.g. optimal transport) for the first strategy, and investigating robust regularization techniques and phase-shift based waveform inversion projection operators for the second strategy.

The work presented in this manuscript is a small part of the larger workflow shown in Figure 101 for reducing structural interpretation uncertainties using seismic data. The long term objective is to use seismic data to update a structural model built from seismic images. There are three main problems in this global workflow: (1) generating seismic data from structural models, (2) computing appropriate data misfit values, (3) and updating the structural model in a manner that reduces data misfit values. Each of these three problems has its own challenges.

The challenge in generating synthetic data from structural models is to create macro-layared velocity models that have structural discontinuities while keeping the velocity kinematically accurate as much as possible (section III.2.5); this is still a problem under investigation. The main focus of this manuscript was to study the second problem of proposing appropriate data misfit functions. The connection of the work presented here and the global workflow in Figure 101 can be seen by the following experiment: consider a sequence of data misfit minimizing structural models generated by the loop in Figure 101; if the data misfit satisfies equation 121, then that sequence converges to the "true structural model", at least in theory. As we have seen in this manuscript, the problem of designing such data misfit functions is also challenging and we still consider it unsolved. The third problem, which we have not yet begun to investigate, is also non-trivial: how can we update structural models from a scalar field error map? Introducing the SIGMA library

IV.1.1 Introduction & Motivation

The SIGMA (Seismic Imaging for GeoModel Analysis) library is a collection of programs that my students (students I supervised for a master's project) and I developed over the course of my thesis. The goal of the package is to develop a suite of programs, i.e. a numerical laboratory, to make my research easier. All the results presented in this thesis are reproducible and were obtained using SIGMA.

I started writing the library when I was learning about seismic imaging and back then I was often asked, "Why are you writing seismic imaging code when there are so many open-source seismic imaging packages? Why reinvent the wheel?". Here are some of the main reasons:

• Programming an algorithm is a reliable way for me to evaluate how well I have understood the algorithm. For example, I don't consider myself to know how seismic imaging works until I have written a seismic imaging code that works.

• Sometimes it takes me about the same amount of time to write a new code as to read and understand an existing code.

• Sometimes the existing code was written to solve a problem similar to mine but not exactly. So I find myself making ad hoc modifications to adapt the code to my problem. The ad hoc code can be quite substantial and usually results in a badly organized code.

• I like to code! Some chapters in this thesis come from my "reinventing the wheel". For example, the generalized Kirchhoff imaging operator proposed in chapter I.3 was motivated by the need to avoid code duplication when implementing different Kirchhoff extended operators (i.e. one code for shot imaging, one code for angle imaging, ...), and the need to avoid shooting rays from each point in the subsurface for Kirchhoff angle imaging, as is common practice. The regularization operators proposed in chapter II.1 was motivated by the need to avoid implementing boundary conditions as this led to a lot of "if statements" (remember that faults are also boundaries), resulting in an untidy code. Had I not "reinvented the wheel", chapters I.3, II.1, II.2, and II.3 would be missing from this manuscript.

IV.1.2 Hardware requirements

The current version of the library runs on Linux, so we'll need a Linux computer. The computer should preferably be equipped with a CUDA-enabled GPU, some programs in the tutorials presented in the following sections run on GPU1 . Finally, the computer should support OpenGL 4.5 or higher, needed by SIGMA's viewer.

IV.1.3 Getting started

There are a couple of ways to benefit from SIGMA's services. First, developers can include it directly in their own projects. However, for recurring practice-oriented tasks, it may be more convenient to compile a set of stand-alone programs to use as commands. We will follow the latter approach in this tutorial.

Let us start by opening a terminal and then launch the SIGMA viewer by invoking the command "S2GoViewer". We then open a second terminal and place it below the viewer. Our screen should look as shown in Figure 102. The terminal will be used not only to run stand-alone programs but also to interact with the viewer. That's right, the viewer can receive commands from the decoupled Linux terminal! This has proven to be very convenient, as it eliminates the need to transfer data into thirdparty viewers at the end of each computation to confirm that results are good.

To conclude this section, it should be pointed out that commands addressed to stand-alone programs are prefixed with "S2Go ",i.e. "S2Go CommandName", while commands addressed to the viewer are prefixed with "S2GoV ", i.e. "S2GoV CommandName". Commands addressed to the viewer are prefixed with "S2GoV ", i.e. "S2GoV CommandName".

IV.1.4 The RSF file format

We will use a Regularly Sampled File (RSF) file-format very similar to the one used by Madagascar [START_REF] Fomel | Madagascar: open-source software project for multidimensional data analysis and reproducible computational experiments[END_REF]. An RSF file is simply an nd-cube. For instance, a 2d RSF file is a section (2d regular grid), and a 3d RSF file is a cube (3d regular grid, a voxet). I will present the RSF file format as used by SIGMA. I distinguish three parts of an RSF file: (1) 2 nd order in time. One of such commands is "S2Go 2dAWECD gpu" (short for 2d Acoustic Wave Equation with Constant Density, runs on GPU); a CPU version of the command is also available: "S2Go 2dAWECD cpu". The result of this command using the velocity model and source in Figure 105.b, and the source wavelet created in the last section, are shown in Figure 106; the experiment can be reproduced using the bash script in Figure 107. Now that we have all the tools needed to run simulations, let us simulate a seismic survey. We start by creating an array of sources and receivers as shown in Figure 108 and in Figure 109 respectively. Then, using the source wavelet created earlier, we invoke the "S2Go 2dAWECD gpu" command discussed earlier. The shot gathers computed using the source-receiver configuration and command from this paragraph are shown in Figure 110.

IV.2.2.2 A seismic survey

IV.2.3 Seismic imaging

The first thing we need to do is to remove the direct arrivals in a our data. We do this by making a copy of our initial data and then using the command shown in Figure 111. Then, we smooth our velocity section to simulate an imaging velocity model. This is done by the command shown in Figure 112. In practice, the imaging velocity model is estimated from data by migration velocity analysis (Woodward et al., 2008;Beraud et al., 2018), but this is beyond the scope of this tutorial.

I(x) is the output image; x = (x, z) is a point in the image domain; x s = s and x r = r are source and receiver coordinates respectively; W (s, r, x) is a weight function; τ (s, r, x) is the traveltime from a ray from s, scattered off x, then recorded at r; and D(s, r, t) is the input data set. The traveltime maps needed to evaluate τ (s, r, x) are computed using the smooth velocity model by the command "S2Go 2dTTimes" as shown in Figure 113; this command solves the eikonal equation using the fast marching method (Sethian and Popovici, 1999). In particular, τ (s, r, x) = T (x, s) + T (x, r), where T is the traveltime cube in Figure 113.

IV.2.3.1 Common image gathers

In practice, it is useful to construct extended images I(x, λ) obtained by extending the image domain x to (x, λ), that is I(x) → I(x, λ) (Symes, 2008). Here, I present extended imaging as described in chapter I.3. Let us rewrite Eq.162 as

I(x) =
for some function ζ(s, r, x). The physical meaning of the parameter λ will depend on our definition of ζ(s, r, x). Extended images are also known as common image gathers. The advantage of image extension (e.g. Eq.164) is that it outputs a collection of structural images as a function of λ. Because a structural image is expected to be invariant if imaged using an accurate velocity, any depth variation of a reflector along the λ-axis suggests an error in the imaging velocity model; this is the basis of some migration velocity analysis methods.

Shot common image gather

The shot common image gather is obtained from Eq.164 by defining ζ(s, r, x) := s.

(165)

The shot common image gather of our data set (Figure 111) is shown in Figure 114 and was obtained by the command "S2Go 2dFAKDMSCIG gpu" (short for 2d First Arrival Kirchhoff Depth Migration, Shot CIG, runs on GPU).

The offset common image gather of our data set (Figure 111) is shown in Figure 115 and was obtained by the command "S2Go 2dFAKDMOCIG gpu" (short for 2d First Arrival Kirchhoff Depth Migration, Offset CIG, runs on GPU). The second example is the interpolation of geological data from the Ribaute model (Caumon et al., 2009). Figure 124a shows the input data. (Ramón et al., 2015). This experiment can be reproduced using the bash script in Je propose un ensemble de conditions qui doivent être satisfaites par la fonction-objectif pour une évaluation réussie des modèles structuraux à l'aide des données sismiques. Ces conditions imposées à la fonction-objectif peuvent, en principe, être satisfaites en utilisant les données sismiques de surface ("surface seismic data"). Cependant, en pratique il reste tout de même difficile de proposer et de calculer des fonctions-objectifs qui satisfassent ces conditions. Je termine le manuscrit en illustrant les difficultés rencontrées en pratique lorsque nous cherchons à évaluer les interprétations structurales à l'aide des donnés sismiques de surface. Je propose une fonction-objectif générale faite de deux composants principaux: (1) un opérateur de résidus qui calcule les résidus des données, et (2) un opérateur de projection qui projette les résidus de données depuis l'espace de données vers l'espace physique (le sous-sol). Cette fonction-objectif est donc localisée dans l'espace car elle génère des valeurs en fonction de l'espace. Cependant, je ne suis toujours pas en mesure de proposer une implémentation pratique de cette fonction-objectif qui satisfasse les conditions imposées pour une évaluation réussie des interprétations structurales; cela reste un sujet de recherche.

Mots-clés: incertitudes structurales, modélisation structurale, imagerie sismique, interprétation sismique, problèmes inverses Towards Reducing Structural Interpretation Uncertainties Using Seismic Data Abstract: Subsurface structural models are routinely used for resource estimation, numerical simulations, and risk management; it is therefore important that subsurface models represent the geometry of geological objects accurately. The first step in building a subsurface model is usually to interpret structural features, such as faults and horizons, from a seismic image; the identified structural features are then used to build a subsurface model using interpolation methods. Subsurface models built this way therefore inherit interpretation uncertainties since a single seismic image often supports multiple structural interpretations. In this manuscript, I study the problem of reducing interpretation uncertainties using seismic data. In particular, I study the problem of using seismic data to determine which structural models are more likely than others in an ensemble of geologically plausible structural models. I refer to this problem as appraising structural models using seismic data. I introduce and formalize the problem of appraising structural interpretations using seismic data. I propose to solve the problem by generating synthetic data for each structural interpretation and then to compute misfit values for each interpretation; this allows us to rank the different structural interpretations. The main challenge of appraising structural models using seismic data is to propose appropriate data misfit functions. I derive a set of conditions that have to be satisfied by the data misfit function for a successful appraisal of structural models. I argue that since it is not possible to satisfy these conditions using vertical seismic profile (VSP) data, it is not possible to appraise structural interpretations using VSP data in the most general case.

The conditions imposed on the data misfit function can in principle be satisfied for surface seismic data. In practice, however, it remains a challenge to propose and compute data misfit functions that satisfy those conditions. I conclude the manuscript by highlighting practical issues of appraising structural interpretations using surface seismic data. I propose a general data misfit function that is made of two main components: (1) a residual operator that computes data residuals, and (2) a projection operator that projects the data residuals from the data-space into the image-domain. This misfit function is therefore localized in space, as it outputs data misfit values in the imagedomain. However, I am still unable to propose a practical implementation of this misfit function that satisfies the conditions imposed for a successful appraisal of structural interpretations; this is a subject for further research.

Keywords: structural uncertainties, structural modeling, seismic imaging, seismic interpretation, inverse problems

Figure 1 :

 1 Figure 1: Exemple d'un modèle structural 3D du sous-sol. Le modèle est composé de deux surfaces d'horizon et dix surfaces de faille.

Figure 2 :

 2 Figure 2: Le concept de la modélisation structurale implicite. a) Données initiales (fournies par Total); les lignes noires ont été interprétées à partir d'une image sismique par un processus similaire à celui illustré dans la Figure 3b. b) Les données initiales sont ensuite interpolées pour obtenir un champ scalaire qui représente le temps stratigraphique. Le champs scalaire servira enfin pour construire le modèle dans la Figure 1 par extraction des isovaleurs: les surfaces d'horizons sont des isovaleurs du champs scalaire. Voir Caumon et al. (2009), par exemple, pour la construction des surfaces de failles.

Figure 3 :

 3 Figure 3: Illustration des incertitudes d'interprétation structurale. a) Image sismique 2D du sous-sol. b) Deux interprétations structurales par deux géologues différents. c) Les modèles géologiques 2D du sous-sol issus des interprétations structurales.

  plus d'informations sur la RTM) en tant que préconditionneur pour accélérer l'inversion des formes d'ondes (voir par exemple Virieux and Operto, 2009, pour plus d'informations sur l'inversion des formes d'ondes). L'idée de base repose sur l'observation que l'image RTM donne une bonne estimation de l'illumination sismique dans le sous-sol pour une géométrie d'acquisition donnée. Cette méthode a le mérite d'être simple à implémenter tout en restant numériquement efficace comme l'illustre la Figure 4. Deuxièmement, dans le chapitre I.3, je propose une formule généralisée pour faire du "Kirchhoff extended imaging" (voir par exemple Symes, 2008, pour plus d'informations sur le concept d'"extended imaging").

Figure 4 :

 4 Figure 4: Accélération de l'inversion des formes d'ondes par la RTM. a) Modèle de vitesses de reference. b) Modèle de vitesses initial. c) Modèle de vitesses obtenu après inversion sans préconditionneur. d) Modèle de vitesses obtenu après inversion en utilisant la RTM comme préconditionneur.

  plus d'informations sur le "Kirchhoff imaging") standard ζ(s, r, x) F (s, r, x),

Figure 5 :

 5 Figure 5: La modélisation structurale implicite par différences-finies. a) Données 2D initiales (fournies par ExxonMobil). b) Interpolation structurale 2D par différences-finies. c) Données 3D initiales (fournies par l'IFPEN). d) Interpolation structurale 3D par différences-finies.

Figure 6 :

 6 Figure 6: Stratégie globale pour l'évaluation des modèles structuraux par données sismiques. a) Données observées utilisées pour l'imagerie sismique. b) Modèles structuraux construits pour chaque interprétation structurale. c) Modèles de vitesses contruits pour chaque modèle structural. d) Données synthétiques calculées pour chaque modèle de vitesses, permettant ainsi de calculer la fonction-objectif pour chaque interprétation structurale.

Figure 7 :

 7 Figure 7: Example of a 3D subsurface structural model. The model is made of two horizon surfaces and ten fault surfaces.

Figure 8 :

 8 Figure 8: The concept of structural implicit modeling. a) Input data (courtesy of Total); the black lines have been interpreted from a seismic image by a process similar to that illustrated in Figure 9. b) The input data are then interpolated to obtain a scalar field which represents the stratigraphic time. The scalar field is then finally used to build the model in Figure 7 by extracting isovalues: horizon surfaces are isovalues of the scalar field. See Caumon et al. (2009), for example, on how to build fault surfaces.

Figure 9 :

 9 Figure 9: Illustration of structural interpretation uncertainties. a) A 2D seismic image of the subsurface. b) Two structural interpretations from two different geologists. c) Subsurface geological models resulting from the structural interpretations.

Figure 10 :

 10 Figure 10: Speeding up waveform inversion with RTM. a) Reference velocity model. b) Initial velocity model. c) Velocity model obtained after inversion without preconditioning. d) Velocity model obtained after inversion using RTM as a precontioner.

Figure 11 :

 11 Figure 11: Finite difference structural implicit modeling. a) 2D input data (courtesy of ExxonMobil). b) 2D finite difference structural interpolation. c) 3D input data (courtesy of IFPEN). d) 3D finite difference structural interpolation.

Figure 12 :

 12 Figure 12: Global workflow for appraising structural models using seismic data. a) Observed data is used for seismic imaging. b) Structural models are built for each interpretation of the seismic image. c) Velocity models are built for each structural model. d) Synthetic data are generated for each velocity model, allowing to compute a data misfit value for each interpretation.

  . Structural modeling has been around at least since the early nineties and it is now a well established field. The reader interested in learning more on structural modeling is referred to Mallet (2002); Caumon et al. (2004); Chilès et al. (2004); Calcagno et al. (2008); Caumon et al. (2009, 2013); Jackson et al. (2013); Mallet (2014); Wellmann and Caumon (2018).

Figure 13 :

 13 Figure 13: Example of a 3D structural model (data courtesy of Total). a) Input data consisting of 10 fault surfaces and horizon points picked along 3 horizons. b) Structural model built from the input data.

Figure 14 :

 14 Figure 14: Structural modeling from seismic interpretation. a) Input data picked on the seismic image. b) A 2D structural model built from the input data.

Figure 15 :

 15 Figure 15: Typical seismic imaging workflow. a) Reference velocity model, unknown in practice. b) "Observed" data generated from the reference model. c) Imaging (migration) velocity model estimated from the observed data. d) Seismic image of the reference model obtained by migrating the observed data using the imaging velocity model.

Figure 16 :

 16 Figure 16: Velocity modeling by ray-based migration velocity analysis. a) Initial velocity model. b) Final velocity model obtained using a low regularization weight. c) Final velocity model obtained using a moderate regularization weight. d) Final velocity model obtained using a high regularization weight.

Figure 17 :

 17 Figure 17: a) Seismic image and common-image-gathers obtained using the initial velocity model in Figure 16a. b) Seismic image and common-image-gathers obtained using the velocity model in Figure 16b. c) Seismic image and common-image-gathers obtained using the velocity model in Figure 16c. d) Seismic image and common-image-gathers obtained using the velocity model in Figure 16d.

Figure 18 :

 18 Figure 18: Structural uncertainties due to the non-uniqueness of the imaging velocity model. Red) Horizon picking on the seismic image in Figure 17d. Blue) Horizon picking on the seismic image in Figure 17c. Yellow) Horizon picking on the seismic image in Figure 17b.

Figure 19 :

 19 Figure 19: A different structural interpretation of the seismic image in Figure 14a. a) Input data picked on the seismic image. b) A 2D structural model built from the input data.

  2. Part II: Notes on Structural Modeling. 3. Part III: Reducing Structural Uncertainties using Seismic Data. Part I presents seismic imaging methods that are considered to be useful for Part III; it's therefore recommended to read Part I before reading Part III. Part II, on the other hand, is independent and can safely be skipped as it is not essential to Part III.There is a fourth Part, the last Part of the manuscript, where I present SIGMA, a software that my students (students I supervised for a master's project) and I developed over the course of my thesis. All the results presented in this thesis are reproducible and were obtained using SIGMA.

Figure 20 :

 20 Figure 20: Workflow for appraising structural models using seismic data. a) Observed data is used for seismic imaging. b) Structural models are built for each interpretation of the seismic image. c) Velocity models are built for each structural model. d) Synthetic data are generated for each velocity model, allowing to compute a data misfit value for each interpretation.

Figure 21 :

 21 Figure 21: Waveform inversion of synthetic data computed in a portion of the Marmousi model (Figure 22a). a) Input seismic traces. b) Subsurface image obtained after waveform inversion. The initial velocity model is shown in Figure 22b.

  The model is discretized on a grid of size 200 × 215 with dx = 30 m and dz = 15 m, and data are computed with a Ricker source wavelet with a dominant frequency of 15 Hz. The data set consists of 60 shot-gathers recorded at 200 receivers; both sources and receivers are evenly distributed along the free surface.

Figure 22 :

 22 Figure 22: Portion of the Marmousi velocity model (Versteeg, 1994) used for numerical experiments. a) The reference (true) velocity model. b) The migration velocity model obtained by smoothing the reference model.

Figure 23 :

 23 Figure 23: Linearized seismic imaging. a) Subsurface image obtained by linearized inversion after 10 iterations (least-squares reverse time migration). b) Subsurface image obtained by a single iteration of linearized inversion (reverse time migration).

Figure 24 :

 24 Figure 24: Nonlinear seismic imaging (waveform inversion). a) Subsurface image obtained using conjugate gradients with a preconditioner. b) Subsurface image obtained using conjugate gradients without a preconditioner.

Figure 25 :

 25 Figure 25: Linearized seismic imaging, 10 iterations. a) Subsurface image obtained by shotwise inversion followed by stacking. b) Subsurface image obtained by simultaneously inversion of all shot gathers.

Figure 26 :

 26 Figure 26: Some common post-processing steps in seismic imaging. a) Reverse-time migration image (Figure 23b before clipping). b) Automatic gain control (here, normalization by the smoothed absolute value of the initial image). c) Smoothing. d) Clipping.

Figure 27 :

 27 Figure 27: Ray-based seismic imaging. a) Subsurface image obtained by linearized inversion after 10 iterations (least-squares Kirchhoff migration). b) Subsurface image obtained by a single iteration of linearized inversion (Kirchhoff migration). These images have been post-processed as described in Figure 26.

Figure 28 :

 28 Figure 28: Synthetic model used for numerical experiments. a) Reference velocity model. b) Imaging model obtained by smoothing the reference model. c) Reference reflectivity obtained by taking the difference between the two previous models. d) Raw output of RTM. e) Illumination map estimated by smoothing the RTM image. f ) The RTM image, normalized using the illumination map.

Figure 29 :

 29 Figure 29: Convergence rate for linearized waveform inversion. a) Estimated reflectivity, nonpreconditioned conjugate gradient. b) Estimated reflectivity, preconditioned conjugate gradient. c) Blue: non-preconditioned; red: preconditioned.

Figure 30 :

 30 Figure 30: Convergence rate for nonlinear waveform inversion with data filtered (low pass) up to a 10Hz cutoff frequency. a) Estimated velocity, non-preconditioned nonlinear conjugate gradient. b) Estimated velocity, preconditioned nonlinear conjugate gradient. c) Blue: non-preconditioned; red: preconditioned.

Figure 31 :

 31 Figure 31: Marmousi velocity model (Versteeg, 1994). a) Reference velocity model. b) Imaging velocity model obtained by smoothing the reference model.

Figure 32 :

 32 Figure 32: Nonlinear waveform inversion with data filtered (low pass) up to 6, 8, 10, and a 12 Hz cutoff frequency. The final velocity model inverted for one frequency is used as the initial velocity model for the next frequency. a) Estimated velocity, non-preconditioned conjugate gradient. b) Estimated velocity, preconditioned conjugate gradient.

Figure 33 :

 33 Figure 33: Linearized waveform inversion (least squares reverse time migration). a) Image obtained with non-preconditioned conjugate gradient. b) Image obtained with preconditioned conjugate gradient. c) Residual curves; blue: non-preconditioned, red: preconditioned.

)

  When λ = 0 and τ = 0, equation 55 reduces to equation 47. Based on equation 48 and equation 55, we propose the more general Kirchhoff extension formula I(x, β, λ, τ ) = s,r δ βζ(s, r, x, λ) F 2 (s, r, x, λ, τ ). (56) I.3.4 Imaging with multiple arrivals

Figure 34 :

 34 Figure 34: Traveltimes computed in the model of Figure 35b. a) First arrival traveltimes section computed by fast marching(Sethian and Popovici, 1999). b) Multi-valued traveltimes section computed by wavefront construction[START_REF] Vinje | Traveltime and amplitude estimation using wavefront construction[END_REF]; only ray cells are displayed. c), d) Single-valued traveltimes surfaces obtained by expressing the multi-valued traveltimes section b) in (offset,depth,take-off angle) coordinates.

Figure 35 :

 35 Figure 35: Illustrative example for offset extension (i.e. using equation 50). a) Reference velocity model. b) Migration model obtained by smoothing the reference model. c) Synthetic data computed in the reference model; sources and receivers cover the free surface. d) Offset extension using first arrivals. e) Offset extension using multiple arrivals. f ) First arrival image. g) Multi-arrival image. An improved continuity of reflectors below the strong velocity anomaly is observed in the multi-arrival image. 57

Figure 36 :

 36 Figure 36: a) Reference velocity model. b) Initial velocity model for ray-based migration velocity analysis.

Figure 37 :

 37 Figure 37: Typical ray-based migration velocity analysis workflow(Stork, 1992; Woodward et al., 2008). a) Offset extended image using the initial velocity model (Figure36b). b) Residual moveout picking by parabola fitting. c) Ray-tracing from residual moveout curves towards the surface. d) Updated model after 15 iterations.

  Figure 37: Typical ray-based migration velocity analysis workflow(Stork, 1992; Woodward et al., 2008). a) Offset extended image using the initial velocity model (Figure36b). b) Residual moveout picking by parabola fitting. c) Ray-tracing from residual moveout curves towards the surface. d) Updated model after 15 iterations.

Figure 38

 38 Figure38compares common image gathers computed in the initial model (Figure36b) with common image gathers computed in the updated model (Figure37d). We clearly note that the initial gathers (Figure38a) are flattened (Figure38b) by the updated velocity model. Figure

  ) Different definitions of ζ in equation 63 give different types of MVA; for example, equation 49 leads to shot-based MVA, equation 50 leads to (surface) offset-based MVA, and equation 51 leads to (reflection) angle-based MVA.

Figure 40 :

 40 Figure 40: Implicit structural modeling workflow. a) Interpolation of input data to build a stratigraphic function. b) Extraction of horizons from the isovalues of the stratigraphic function.

  ∇φ * (x) away from observations (b) ||∇φ * (x)|| > 0 for all x (c) d dl φ * (x) is monotonic along any line joining two successive horizons (d)

Figure 41 :

 41 Figure 41: Benchmark models; the white points represent input data. a) Benchmark model proposed by Renaudeau et al. (2017a). b) Benchmark model proposed by Renaudeau et al. (2017b).

Figure 42 :

 42 Figure 42: Interpolation of the data in Figure 41a. a) Results obtained using the Laplacian roughness operator. b) Results obtained using a more appropriate roughness operator for structural modeling.

Figure 43 :

 43 Figure 43: Interpolation of the data in Figure 41b. a) Results obtained using an isotropic roughness operator. b) Results obtained when we explicitly take into account the orientation of horizons.

Figure 44 :

 44 Figure 44: Implicit stratigraphic field computed using the proposed method. a) Input data consisting of 10 fault surfaces and horizon points picked along 3 horizons. b) Solid view of the resulting stratigraphic field. c) Section view of the resulting stratigraphic field. Data courtesy of Total.

Figure 45 :

 45 Figure 45: Implicit stratigraphic modeling of a section from the Ribaute model of Caumon et al. (2009). a) Input data, three faults and three horizons. b) Resulting stratigraphic field.

  Figure 46: a) 1D domain discretization. b) 2D domain discretization. The red nodes are external nodes, the blue nodes are boundary nodes, and the black nodes are internal nodes.

Figure 47 :

 47 Figure 47: Notion of neighborhoods and directional vectors in 1D and 2D. a) 1D neighborhood: the neighborhood of the black node is made of the two (green) points that can be reached from it by taking one step (one forward and one backward) along the 1D directional vector. b) 2D neighborhood: the neighborhood of the black node is made of the eight (green) points that can be reached from it by taking one step (one forward and one backward) along the four 2D directional vectors.

  II.1.4 Interpolation in 1D: problem formulation II.1.4.1 Interpolation in 1D

Figure 48 :

 48 Figure 48: Implicit stratigraphic modeling of a modified section from one of the 2D benchmark models of Renaudeau et al. (2017a). a) Input data, two horizons. b) Resulting stratigraphic field.

Figure 49 :

 49 Figure 49: Implicit stratigraphic modeling of the Balzes fold model of Ramón et al. (2015). a) Input data consisting of horizon points picked along two horizons. b) Solid view of the resulting stratigraphic field. c) Section view of the resulting stratigraphic field.

Figure 50 :

 50 Figure 50: Domain discretization in the presence of a discontinuity. The red line and nodes represent a discontinuity and its rasterization. Red nodes are tagged as external nodes. Blue nodes are boundary nodes, and black nodes are internal nodes.

Figure 51 :

 51 Figure 51: Implicit stratigraphic modeling of model in Figure 48, in the presence a discontinuity. a) Input data, two horizons and one fault. b) Resulting stratigraphic field.

Figure 52 :

 52 Figure 52: Implicit stratigraphic modeling in a highly faulted model. a) Input data, 1400 data points and 50 faults. b) Resulting stratigraphic field. Synthetic data courtesy of ExxonMobil.

Figure 53 :

 53 Figure 53: Implicit stratigraphic modeling of the 3D sandbox model of Chauvin et al. (2018). a) Input data consisting of 27 fault surfaces and horizon points picked along 6 horizons. b) Solid view of the resulting stratigraphic field. c) Section view of the resulting stratigraphic field. Data courtesy of IFPEN and C&C Reservoirs.

Figure 54 :

 54 Figure 54: Implicit stratigraphic modeling of the 2D benchmark model of Renaudeau et al. (2017a).The input data, the white curves, represent two different iso-values. a) Result obtained on a triangle mesh using Laplacian regularization (eq. 112). b) Result obtained on a triangle mesh using the proposed regularization operator (eq. 117).

Figure 55 :

 55 Figure 55: Example of valid and invalid meshing for structural interpolation. a) Invalid mesh: nodes 3 and 7 will not be constrained because there is no edge connecting them to an internal node. b) Valid mesh: all the elements have at least one internal node.

Figure 56 :

 56 Figure56: Implicit stratigraphic modeling of a section from the Ribaute model ofCaumon et al. (2009). a) Result obtained on a triangle mesh using Laplacian regularization (eq. 112). b) Result obtained on a triangle mesh using the regularization operator (eq. 114). c) Result obtained on a triangle mesh using the proposed regularization operator (eq. 117). d) Result obtained on a quadrilateral mesh using the proposed regularization operator (eq. 117).

Figure 58 :

 58 Figure 58: Handling large thickness variations by explicitly taking into account the orientation of horizons. This benchmark model is from Renaudeau et al. (2017b). a) Interpolation result obtained using the isotropic regularization operator proposed in chapter II.1. b) Estimation of the stratigraphic vector field from input data. c) Interpolation result obtained using the isotropic regularization operator proposed in chapter II.1 in addition to constraining the stratigraphic field to follow the estimated vector field. d) Interpolation result obtained using an anisotropic version of the regularization operator proposed in chapter II.1; anisotropic smoothing weights are determined from the estimated vector field.

Figure 59 :

 59 Figure 59: Illustration of structural interpretation uncertainties. The seismic image above supports N geologically consistent structural interpretations. The problem of appraising structural interpretations using seismic data consists in determining which among the N structural models are more likely than others using seismic data.

Figure 60 :

 60 Figure 60: Seismic uncertainties due to complex geology. a) Thinlayered velocity model. b) Macrolayered velocity model used to build the thinlayered model. c) Seismic image (reverse time migration) computed from data generated in the thinlayered model. d) Seismic image (reverse time migration) computed from data generated in the macrolayered model. The problem of structural interpretation is that of identifying the structural surfaces observed in the macrolayered model b) from the seismic image in c).

Figure 61 :

 61 Figure 61: Seismic image generated from the sandbox model ofColletta et al. (1991)

Figure 63 :

 63 Figure 63: a) Structural interpretation 2. b) Structural model 2.

Figure 64 :

 64 Figure 64: a) Structural interpretation 3. b) Structural model 3.

Figure 65 :

 65 Figure 65: a) Structural interpretation 4. b) Structural model 4.

Figure 66 :

 66 Figure 66: a) Structural interpretation 5. b) Structural model 5.

Figure 67 :

 67 Figure 67: a) Structural interpretation 6. b) Structural model 6.

Figure 68 :

 68 Figure 68: a) Structural interpretation 7. b) Structural model 7.

Figure 69 :

 69 Figure 69: a) Structural interpretation 8. b) Structural model 8.

Figure 70 :

 70 Figure 70: The reference velocity model. It is derived from the sandbox model of Colletta et al. (1991).

Figure 71 :

 71 Figure 71: a) Migration velocity model obtained by smoothing the slowness of the reference model in Figure 70. b) Depth migrated image (Kirchhoff depth migration) of data computed from the reference model. The black box is the region of interest.

Figure 72 :

 72 Figure 72: a, b) Two possible structural interpretations/models of Figure 71b from different interpreters. c, d) Macro-layered velocity models built from the migration velocity model and the structural models.

Figure 73 :

 73 Figure 73: a) Reference model (M ref ). b) Migration model (M mig ) obtained by smoothing the slowness of the reference model.

Figure 74 :

 74 Figure 74: a) A candidate interpretation represented as a block macro-layered velocity model (M b ) built using the migration model in Figure 73b, as described by equation 122. The red arrow identifies the part of the horizon that was mispositioned during interpretation.b) Model residual, difference between a) and the reference model in Figure 73a. The red arrow identifies the interpretation error while the green arrows identify picking errors.

Figure 75 :

 75 Figure 75: The same candidate interpretation from Figure 74a represented using a wire macro-layered model (M w ), as described by equation 123.

  Mb is a smooth velocity model obtained by smoothing the slowness of M b . The idea behind the reflectivity macro-layered model comes from the Born approximation which approximates the true velocity model by the sum of a smooth background model plus a reflectivity function (Bleistein et al., 2001, chapter 2); the method proposed here attempts to replace the interpretation-dependent smooth background model Mb with the migration velocity model M mig .

Figure 76 :

 76 Figure 76: The same candidate interpretation from Figure 74a represented using a reflectivity macrolayered model (M r ), as described by equation 124.

Figure 77 :

 77 Figure 77: Geometrical interpretation of the VSP misfit weight function. a) Receivers are positioned only in one well. The green rays correspond to source-receiver pairs that will actually contribute to the data-misfit 125. According to equation 129, data-misfit values are projected along entire rays since every point x along the green rays has ω(s, r, x) = 1. b) Receivers are positioned in two wells, thus illuminating the region of interest from different directions. Now points x along the green rays have ω(s, r, x) = 2 if x ∈ Ω, otherwise ω(s, r, x) = 1.

Figure 79 :

 79 Figure 79: Data residuals of two traces from Figure 78. a) Comparison between observed (red) against synthetic (blue) traces. b) Difference between observed and synthetic traces. The top green arrow identifies picking errors (δp), the bottom green arrow identifies picking errors in addition to velocity errors (δv) in the macro-layered model in Figure 74a, and the red arrow identifies horizon mispositioning error (δh). c) Same as b) but using synthetic data from the more kinematically accurate macro-layered model in Figure 76. Velocity errors are removed as highlighted by the green circle. d) Dynamic-time-warping phase residual using synthetic data from the model in Figure 76.

Figure 80 :

 80 Figure 80: a) Projection of L 1 data residuals (equation 131, with p=1) into the image-space using L * (defined in equation 134). Synthetic data were computed in the macro-layered model in Figure 74a. Green arrows identify picking errors while the red arrow identifies velocity errors. b) Projection of L 1 data residuals into the image-space using L * . Synthetic data were computed in the more accurate macro-layered model in Figure 76. The effect of velocity errors is removed.

Figure 81 :

 81 Figure 81: a) Projection of Dynamic-time-warping data residuals (equation 132) into the image-space using KDM (defined in equation 162). Synthetic data were computed in the macro-layered model in Figure 76. The black line highlights the bipolarity of the residuals by separating positive residuals above from negative residuals below. b) Projection of Dynamic-time-warping data residuals into the image-space using L * (defined in equation 134). Synthetic data were computed in the macro-layered model in Figure 76.

Figure 82 :

 82 Figure 82: Eight possible structural interpretations/models of the seismic image in Figure 71b from different interpreters.

Figure 83 :

 83 Figure 83: Appraising the structural models in Figure 82 using VSP data. Each point represents a structural model from Figure 82. The "observed" data were computed in the reference model in Figure 70, with sources positioned at the free surface and receivers positioned in a vertical well at x = 7500 m. a) The data-space ranking obtained using the L 1 -norm of the entire data set is not consistent with the ranking expected in the model-space. b) The data-space ranking obtained using the localized data-misfit 125, along with equation 127 and equation 129, is consistent with the ranking expected in the model-space. The region of interest Ω is the black box in Figure 71b.

Figure 84 :

 84 Figure 84: a) Reference model. b) Candidate model from interpretation. c) Difference between reference (a) and candidate model (b).

Figure 85 :

 85 Figure 85: a) "Observed data" computed in the reference model in Figure 84a. b) Synthetic data computed in the candidate model in Figure 84b. c) Difference between observed (a) and synthetic data (b).

Figure 86 :

 86 Figure 86: a) Migration velocity model obtained by smoothing the reference model (Figure 84)a. b) Dynamic image warping (DIW) residuals computed from the observed and synthetic data in Figure 85a and Figure 85b. c) Reverse time migration of the DIW data residuals (b). d) Kirchhoff depth migration of the DIW data residuals (b).

Figure 87 :

 87 Figure 87: a) Reference model. b) Candidate model from interpretation. c) Difference between reference (a) and candidate model (b).

Figure 88 :

 88 Figure 88: a) "Observed data" computed in the reference model in Figure 87a. b) Synthetic data computed in the candidate model in Figure 87b. c) Difference between observed (a) and synthetic data (b).

Figure 89 :

 89 Figure 89: a) Dynamic image warping (DIW) residuals computed from the observed and synthetic data in Figure 88a and Figure 88b. b) Migration velocity model obtained by smoothing the reference model (Figure 87)a. c) Kirchhoff depth migration of the DIW data residuals (a).

Figure 90 :

 90 Figure 90: a) Reference model. b) Candidate model from interpretation. c) Difference between reference (a) and candidate model (b).

Figure 91 :

 91 Figure 91: a) "Observed data" computed in the reference model in Figure 90a. b) Synthetic data computed in the candidate model in Figure 90b. c) Difference between observed (a) and synthetic data (b).

Figure 92 :

 92 Figure 92: a) Dynamic image warping (DIW) residuals computed from the observed and synthetic data in Figure 91a and Figure 91b. b) Migration velocity model obtained by smoothing the reference model (Figure 90)a. c) Kirchhoff depth migration of the DIW data residuals (a).

Figure 93 :

 93 Figure 93: a)Adjoint-source (right-hand side of equation 152b) computed from the traces in Figure and the phase-shifts in Figure 86b. b)Reverse time migration of the adjoint-source in a). c)Reverse time migration of the data-difference in Figure 86c.

Figure 94 :

 94 Figure 94: a)Adjoint-source (right-hand side of equation 152b) computed from the traces in Figure and the phase-shifts in Figure 89b. b)Reverse time migration of the adjoint-source in a). c)Reverse time migration of the data-difference in Figure 89c.

Figure 95 :

 95 Figure 95: a)Adjoint-source (right-hand side of equation 152b) computed from the traces in Figure 91 and the phase-shifts in Figure 92b. b)Reverse time migration of the adjoint-source in a). c)Reverse time migration of the data-difference in Figure 92c.

Figure 96 :

 96 Figure 96: a) Velocity model obtained by adding the migration velocity model in Figure 86a to the error map in Figure 84c. b) Synthetic data generated from the model in a). c) Velocity error obtained by waveform inversion of the synthetic data set. d) Absolute value of the velocity error map in c.

Figure 97 :

 97 Figure 97: Estimating the error map in Figure 84c by waveform inversion. a) Velocity error map obtained using synthetic data generated in a block macro-layered velocity model. b) Absolute value of the error map in a.

Figure 98 :

 98 Figure 98: Estimating the error map in Figure 84c by waveform inversion. a) Velocity error map obtained using synthetic data generated in a reflectivity macro-layered velocity model. b) Absolute value of the error map in a.

Figure 99 :

 99 Figure 99: Estimating the error map in Figure 84c by waveform inversion with total variation regularization. a) Velocity error map obtained using synthetic data generated in a reflectivity macro-layered velocity model. b) Absolute value of the error map in a.

Figure 100 :

 100 Figure 100: Estimating the error map in Figure 87c by waveform inversion. The figures in the first column are velocity error maps; the second column is the absolute value of the first column. a) Target error map obtained by means of the experiment in Figure 96. b) Error map obtained using synthetic data generated in a block macro-layered velocity model. c) Error map obtained using synthetic data generated in a reflectivity macro-layered velocity model. d) Error map obtained using synthetic data generated in a reflectivity macro-layered velocity model plus total variation regularization.

Figure 101 :

 101 Figure 101: Long term global workflow for reducing structural interpretation uncertainties using seismic data. 1: generating synthetic seismic data from structural models. 2: computing data misfit values using synthetic and observed seismic data. 3: updating the structural model in a manner that reduces misfit values.

Figure 102 :

 102 Figure 102: Getting started with SIGMA for the user. We launch the viewer, and then open a Linux terminal below it. The terminal is used not only to run stand-alone commands but also to interact with the viewer.

Figure 104 :

 104 Figure 104: a) Example of a complete header file. The red boxes highlight the only required fields, in addition to those describing the geometry. b) Loading the RSF file. The velocity property values shown in the model are read from the data file pointed to by the "in" field of the header file.

Figure 107 :

 107 Figure 107: Bash script used to produce the example in Figure 106.

Figure 108 :

 108 Figure 108: a) Geometry file describing an array of source coordinates. b) Viewing the source coordinates (highlighted by the red box) described the given geometry file.

Figure 109 :

 109 Figure 109: a) Geometry file describing an array of receiver coordinates. b) Viewing the receiver coordinates (highlighted by the red box) described the given geometry file.

Figure 110 :

 110 Figure 110: Shot gathers recorded at receivers. Each section in the rec-time plane is a shot gather; that is, a collection of data recorded from the same shot.

Figure 111 :

 111 Figure 111: The same shot gather cube of Figure 110 after muting first arrivals (direct waves).

Figure 112 :

 112 Figure 112: Smoothing the velocity section of Figure 104 to simulate an imaging velocity model.

Figure 113 :

 113 Figure 113: a) Geometry file describing source coordinates for computing traveltime maps; it is the same as that in Figure 109. b) Computing traveltime maps at sources described by the given geometry file. Each section in the x-z plane is a traveltime map from one shot. The cube has been displayed with a stratigraphy colormap for a better visualization of wavefronts (S2GoV SetColormap cubes.cubename.propertyname colormap=stratigraphy).

  ζ(s, r, x) F (s, r, x),

Figure 114 :

 114 Figure 114: Shot common image gather. Each section in the x-z plane is an image obtained by migrating data from the same shot. The parameters oe, de, ne are the origin, sampling rate, and size of the "extension" axis; i.e. the shot axis in this figure.

Figure 115 :

 115 Figure 115: Offset common image gather. Each section in the x-z plane is an image obtained by migrating data that have the same offset. The parameters oe, de, ne are the origin, sampling rate, and size of the "extension" axis; i.e. the offset axis in this figure.

Figure 118 :

 118 Figure 118: Common offset migration. This is the result of stacking (summing) the images (x-z planes) in Figure 115 along the offset axis.

Figure 119 :

 119 Figure 119: Common angle migration. This is the result of stacking (summing) the images (x-z planes) in Figure 116 along the angle axis.

Figure 122 :

 122 Figure 122: Bash script used to produce the example in Figure 121.

Figure 123 :

 123 Figure 123: a)Documentation of the command in lines 3 and 7 in Figure 122. b)Documentation of the command in lines 17 and 31 in Figure 122.

  Figure124bshows results obtained using equation 26 in chapter II.1 [equation 5 inIrakarama et al. (2018)]. Figure124cshows results obtained using equation 30 in chapter II.1 [equation 6 inIrakarama et al. (2018)]. The results can be reproduced using the bash script in Figure125. Lines 1 to 11 of the bash script set the correct path of the RSF binary data files in the corresponding RSF headers; the input data are shipped with the tutorial. The only new command in Figure125is the command in line 21, which creates "2D curves file" from a collection of RSF points files. A "2D curve file" is a text file in which each line has 3 entries: x y curve id.

Figure 124 :

 124 Figure 124: Interpolation of the Ribaute model(Caumon et al., 2009). This experiment can be reproduced using the bash script in Figure125. a) Input geometry file (region of interest). a) Input data consisting of three horizons and three faults. b) Results obtained using equation 26 in chapter II.1 [equation 5 in Irakarama et al. (2018)]. c) Results obtained using equation 30 in chapter II.1 [equation 6 in Irakarama et al. (2018)].

Figure 126 :

 126 Figure 126: Interpolation of concentric spheres. This experiment can be reproduced using the bash script in Figure 127. a) Input data consisting of an inner and an outer sphere. b) Results obtained using equation 34a in chapter II.1 [equation 5 in Irakarama et al. (2018)]. c) Results obtained using equation 34b in chapter II.1 [equation 6 in Irakarama et al. (2018)].

Figure 127 :

 127 Figure 127: Bash script used to produce the example in Figure 126.

Figure 128 :

 128 Figure 128: Interpolation of the Balzes fold model(Ramón et al., 2015). This experiment can be reproduced using the bash script in Figure129. a) Input data consisting of two horizons. b) Results obtained using equation 34a in chapter II.1 [equation 5 in Irakarama et al. (2018)]. c) Results obtained using equation 34b in chapter II.1 [equation 6 in Irakarama et al. (2018)].

  Figure 128: Interpolation of the Balzes fold model(Ramón et al., 2015). This experiment can be reproduced using the bash script in Figure129. a) Input data consisting of two horizons. b) Results obtained using equation 34a in chapter II.1 [equation 5 in Irakarama et al. (2018)]. c) Results obtained using equation 34b in chapter II.1 [equation 6 in Irakarama et al. (2018)].

  La troisième partie du manuscrit est consacrée à l'évaluation des modèles structuraux à l'aide des données sismiques. Je propose deux contributions dans ce domaine. Premièrement, dans le chapitre III.2, je développe une théorie quantitative qui formalise le problème d'évaluation des modèles structuraux par données sismiques. La stratégie globale proposée pour résoudre ce problème est illustrée dans la Figure6: l'idée de base est de générer des données synthétiques pour chaque interprétation structurale, puis comparer ces données aux données observées par calcul d'une fonction-objectif adéquate. Je propose un ensemble de conditions qui doivent être satisfaites par la fonction-objectif pour une évaluation réussie, puis j'argumente qu'une fonction-objectif susceptible de satisfaire ces conditions est de forme

  où E(x) désigne la carte d'erreurs d'interprétation, R un opérateur de résidus qui calcule la différence entre les données observées d o et les données synthétiques d i calculées dans le i eme modèle structural, P un opérateur de projection qui projette les résidus des données depuis

l'espace des données vers l'espace physique (le sous-sol); Ω est la region que nous aimerions évaluer dans le sous-sol. Deuxièmement, dans le chapitre III.3, j'étudie les aspects pratiques du problème et montre que cela reste difficile à proposer une implémentation de l'équation 3 qui satisfait les conditions imposées pour une évaluation réussie des modèles structuraux. Malgré mon échec à démontrer une implémentation pratique qui satisfasse les conditoins requises, je conclus que l'implémentation la plus prometteuse de la l'équation 3 est de forme

  Vers la réduction des incertitudes d'interprétation structurale à l'aide de données sismiques Résumé : Les modèles géologiques sont couramment utilisés pour estimer les ressources souterraines, pour faire des simulations numériques, et pour évaluer les risques naturels; il est donc important que les modèles géologiques représentent la géométrie des objets géologiques de façon précise. La première étape pour construire un modèle géologique consiste souvent à interpréter des surfaces structurales, telles que les failles et horizons, à partir d'une image sismique; les objets géologiques identifiés sont ensuite utilisés pour construire le modèle géologique par des méthodes d'interpolation. Les modèles géologiques construits de cette façon héritent donc les incertitudes d'interprétation car une image sismique peut souvent supporter plusieurs interprétations structurales. Dans ce manuscrit, j'étudie le problème de réduire les incertitudes d'interprtation à l'aide des données sismiques. Particulièrement, j'étudie le problème de déterminer, à l'aide des données sismiques, quels modèles sont plus probables que d'autres dans un ensemble des modèles géologiques cohérents. Ce problème sera connu par la suite comme le problème d'évaluation des modèles géologiques par données sismiques. J'introduis et formalise ce problème. Je propose de le résoudre par génération des données sismiques synthétiques pour chaque interprétation structurale dans un premier temps, ensuite d'utiliser ces données synthétiques pour calculer la fonction-objectif pour chaque interprétation; cela permet de classer les différentes interprétations structurales. La difficulté majeure d'évaluer les modèles structuraux à l'aide des données sismiques consiste à proposer des fonctions-objectifs adéquates.

The discussion in section I.3.5 was implemented and tested by Thomas Beraud as part of his master's project. If we had more time, implementing and testing the discussion in section I.3.6 would have made another interesting master's project.

Condition

implies that the data misfit function has no local minima; therefore, the relation between data and model cannot be highly nonlinear. However, because the ratio in condition 2 need not be constant, the relation between data and model need not be linear either.

This chapter resulted from an ongoing technical report co-authored with Paul Cupillard, GuillaumeCaumon, and Paul Sava.

CPU code is available for some commands as well but is considerably slower.

Acknowledgements

Acknowledgements I would like to thank Christelle Butault and Paul Angrand for the seismic interpretations presented in this general introduction, and Thomas Beraud for carrying out the velocity modeling experiment.

II.1.10 Acknowledgements

The authors would like to thank Total for the data used in Figure 44, Hao Huang (ExxonMobil) for the data used in Figure 52, and Benjamin Chauvin (previously RING, currently Harvard), IFPEN, C&C Reservoirs for the data used in Figure 53.

Modeste Irakarama would also like to thank Paul Cupillard (RING) for his encouragement to publish this work.

This work was done in the frame of the RING project at Université de Lorraine. We would therefore like to thank the sponsors of the RING-GOCAD Consortium managed by ASGA for their support. Software corresponding to this chapter is available to consortium members.

III.2.9 Acknowledgements

The authors would like to thank Ludovic Métivier and Pierre Thore for fruitful discussions about the subject. This work also greatly benefited from constructive discussions with Thomas Bodin and Yann Capdeville in the frame of the HIWAI Project. We are very thankful to Elizabeth L'Heureux, Isabelle Lecomte, Stan Davis, and Xinming Wu for the time they dedicated to review this chapter. We would also like to thank Christelle Butault, Francois Bonneau, and Paul Angrand, for some seismic interpretations used in this work. The traveltime maps needed to evaluate Kirchhoff integral operators were computed using the Madagascar software freely available from www.ahay.org. This work was done in the frame of the RING project at Université de Lorraine. The sponsors of the RING-GOCAD Consortium managed by ASGA are hereby acknowledged for their support. We would also like to thank Paradigm for providing the SKUA-GOCAD software used for structural modeling.

III.3.6 Acknowledgements

The authors would like to thank Ludovic Métivier and Pierre Thore for fruitful discussions about the subject. This work also greatly benefited from constructive discussions with Thomas Bodin and Yann Capdeville within the framework of the HIWAI Project.

This work was done in the frame of the RING project at Université de Lorraine. The sponsors of the RING-GOCAD Consortium managed by ASGA are hereby acknowledged for their support. We would also like to thank Paradigm for providing the SKUA-GOCAD software used for structural modeling.

IV.3.3 Acknowledgements

This work was done within the framework of the RING project at Université de Lorraine. I would like to thank the industrial and academic sponsors of the RING-Gocad Consortium managed by ASGA for their support.

data: if a macro-layered velocity model introduces a velocity error, then all reflectors below this error will be mispositioned in the data-space even if they were positioned at the right place in the image-space by the interpreter. This eventually leads to a violation of condition 6 after back-projecting data errors into the image-space. An illustration of these problems is presented in Figure 79 where we compare individual traces. Figure 79a compares a trace from Figure 78a against a trace from Figure 78b; the difference between these two traces is shown in Figure 79b, where it can be noted that picking errors (δp) are as strong as errors due to mispositioning of horizons (δh). Figure 79c shows a similar experiment but using synthetic data computed in the more kinematically accurate reflectivity macro-layered model in Figure 76; in this case, the additional velocity errors are removed as highlighted by the green circle in Figure 79c.

III.3.1 Summary

Structural models built from interpretation of seismic images inherit uncertainties propagated from interpretation uncertainties. As a result, it is often possible to propose multiple structural models from a single seismic image. We propose a method for appraising different geologically valid structural models built from a single seismic image. The main idea of the proposed method is to simulate synthetic data for each structural model, and then to compare the resulting synthetic data to observed data. In this chapter, we discuss the practical aspects of appraising structural models using surface seismic data.

Part IV

Seismic Imaging for GeoModel Analysis -SIGMA the geometry file, (2) the header file, and (3) the data file.

1) The geometry file has a ".hd" extension. "hd" is short for header. The geometry file describes the axes of the grid. An axis is described by the following parameters:

• oi : origin of the i-th axis.

• di : sampling step of the i-th axis.

• ni : number of samples on the i-th axis.

• labeli : label (name) of the i-th axis, optional.

• uniti : units of the i-th axis, optional.

Where "i" is an integer starting from 1. An example of a simple geometry file called "geom example.hd" is shown in Figure 103.a. This geometry file describes a cube of size 100x100x100, with d1=dx=10, d2=dy=5, d3=dz=5, and an origin at the point (0,0,0); Figure 103.b shows one command that can be used to view this cube. 2) The header file also has a ".hd" extension. "hd" is short for header. The header has all the fields of the geometry file and additional ones including:

• in : path to the actual (binary) data file.

• esize : size, in bytes, of each element in the (binary) data file.

An example of a header file named "vel section.hd", and how to view the corresponding data, are shown in Figure 104.a and Figure 104.b respectively. The velocity model in Figure 104.b is derived from the sandbox of model of Colletta et al. (1991).

3) The data file has a ".dt" extension. "dt" is short for data. The data file is in binary format and contains property values of every grid point. It is found from the "in" field of the associated header file (Figure 104).

In this tutorial

In this tutorial, we will test SIGMA's 2D seismic modeling and imaging programs. In particular, we are going to perform the following activities:

• Modeling of acoustic seismic data by finite difference wave simulation.

• Building extended images, i.e. common image gathers (CIGs), from seismic data.

• Stacking CIGs to produce depth migrated images.

IV.2.2 Seismic modeling

IV.2.2.1 Preparing for wave simulation

Source and receiver coordinates First, we need to create an array of source and receiver coordinates before running the simulation. We can create an array of points by defining its geometry file and then invoking the command "S2Go GeomToPointSet", as shown in Figure 105. It is good practice to always view the resulting points to confirm that sources and receivers are positioned where we want them before running simulations.

Source wavelet

Now, we need to create a source wavelet. This can be achieved by invoking the command: "S2Go GetRickerWavelet fmax=15 dt=0.004 tsim=5 tdelay=0.4 > wlt.hd", which creates a Ricker wavelet signal that is 5 seconds long, sampled at 0.004 seconds, centered at 0.4 seconds, with a dominant frequency of 15 Hertz, and saves it in a file called wlt.hd. The acoustic wave simulation commands are implemented using staggered finite differences (Virieux, 1984), perfectly matched layers boundary conditions (Collino and Tsogka, 2001), and run on GPU [START_REF] Micikevicius | 3D finite-difference computation on GPUs using CUDA[END_REF]; they implement 8 th order difference operators in space and

Running the simulation

Angle common image gathers

The angle common image gather is obtained from Eq.164 by defining ζ(s, r, x) := acos ∇T (x, s).∇T (x, r) ∇T (x, s).∇T (x, s) , T is the traveltime cube in Figure 113.

(167)

The angle common image gather of our data set (Figure 111) is shown in Figure 116 and was obtained by the command "S2Go 2dFAKDMACIG gpu" (short for 2d First Arrival Kirchhoff Depth Migration, Angle CIG, runs on GPU).

IV.2.3.2 Stacking common image gathers

IV.2.4 Acknowledgements

This work was done within the framework of the RING project at Université de Lorraine. I would like to thank the industrial and academic sponsors of the RING-Gocad Consortium managed by ASGA for their support.

Chapter IV.3

Tutorial II: Structural Interpolation

IV.3.1 In this tutorial

In this tutorial, we will test SIGMA's finite difference structural interpolation method presented in chapter II.1 and in Irakarama et al. (2018). We start with an illustrative example of interpolating two concentric spheres as shown in Figure 126; the experiment can be reproduced using the bash script in Figure 127. The second example is the interpolation of geological data from the the Balzes fold model (Ramón et al., 2015) as shown in Figure 128; the experiment can be reproduced using the bash script in Figure 129. Note that the data set in Figure 128a is loaded from the "horizons vals.txt" file created in line 17 of Figure 129; such a "3D value file" can be loaded in the viewer using the command "S2GoV LoadFDIVals3d horizons vals.txt".

IV.3.2 Structural interpolation IV.3.2.1 Interpolation in 2D