Nicolas Clausolles

Pauline Collon

Modeste Irakarama

Guillaume Caumon

Stochastic velocity modeling for assessment of imaging uncertainty during seismic migration: application to salt bodies

Keywords:

Variations in the migration velocity model directly affect the position of the imaged reflectors in the subsurface, leading to structural imaging uncertainties. These uncertainties are not explicitly addressed when trying to deterministically build an adequate velocity model. This paper presents a new stochastic geology-controlled velocity modeling method handling the possible presence of a salt weld. This permits to generate a large set of geological scenarios and associated velocity models. Each model is used to remigrate the seismic data. Then, a statistical analysis of the resulting seismic images is performed to quantify the local variability of the seismic responses. The approach is applied to the imaging of salt diapirs, in an iterative scheme (migrate, pick and update). The results show that, similarly to stacking common mid-point gathers, the statistical analysis preferentially preserves recurrent features from an image to another. In particular, this analysis permits to distinguish between connected and detached diapirs without prior knowledge about their connectivity, highlighting the potential of the method to resolve important aspects about basin and reservoir architecture. More generally, it provides quantitative information on the parts of the seismic image most sensitive to migration velocity variations, which opens interesting perspective to quantitative interpretation uncertainty assessment. Finally, the presented application also suggests that it is possible to significantly improve the quality of the generated seismic images by sampling many possible geological scenarios.

I

Obtaining a faithful seismic image of the structures at depth requires the definition of a set of accurate subsurface parameter fields, which are used during the seismic migration. Building a reliable velocity model is, therefore, a seismic processing step in its own right [e.g., [START_REF] Jones | Estimating subsurface parameter fields for seismic migration: Velocity model building[END_REF]. The presence of inaccuracies in the velocity model, and the intrinsic limits of seismic acquisition and imaging result in potentially uncertain seismic images (see [START_REF] Thore | [END_REF] for an overview of seismic uncertainty sources). All authors agree that propagating these uncertainties along the E&P workflow is primordial, as it should be taken into account when assessing risks and making decisions [e.g., [START_REF] Osypov | Model-uncertainty quantification in seismic tomography: Method and applications[END_REF].

Numerous authors have documented the sensitivity of seismic imaging to variations in the migration velocity model [e.g., [START_REF] Grubb | Estimating and interpreting velocity uncertainty in migrated images and AVO attributes[END_REF][START_REF] Paffenholz | SIGSBEE_2a Synthetic Subsalt Dataset -Image Quality as Function of Migration Algorithm and Velocity Model Error[END_REF], Bube et al., 2004a,b]. Some criteria have been proposed to (in)validate a velocity model (such as image gather flatness), and other authors alternatively proposed to use geostatistical approaches to check the spatial consistency (and correct inconsistencies) of the migration velocity fields [e.g., [START_REF] Leron | Spatial Quality Control of Seismic Stacking Velocities Using Geostatistics[END_REF], Sandjivy et al., 2003]. Still, assessing the impact of inaccuracies in the velocity model on the resulting migrated image remains a complex task [START_REF] Guo | Becoming effective velocity-model builders and depth imagers, Part 2 -The basics of velocity-model building, examples and discussions[END_REF]Fagin, 2002, Jones, 2015].

Early works often focused on deriving formulas in simple model configurations to describe the reflection positioning in-accuracies [e.g., [START_REF] Landa | Interpretation of velocity estimates from coherency inversion[END_REF], Bube et al., 2004a,b, Pon and Lines, 2005]. They are, however, not suited for use with real data sets. More recent approaches often aim at determining an envelope around specific structural elements (the faults and horizons of interest) to describe their expected variability. A variety of approaches have been proposed to perform this task. Based on a thorough analysis of the seismic processing chain, [START_REF] Thore | [END_REF] parameterize the uncertainty about the position of structural elements in terms of magnitude, orientation, and spatial correlation. Other authors use physical considerations to define these envelopes, such as the concept of Fresnel zone around faults [START_REF] Weinzierl | Volumetric impact of fault perturbation in the first Fresnel zone: Interpretation[END_REF]. A more computational and common family of approaches consists in introducing small perturbations into the velocity model which is obtained after tomographic inversion, and to remigrate the structural elements of interest using map migration [e.g., [START_REF] Osypov | Model-uncertainty quantification in seismic tomography: Method and applications[END_REF].

As noted by [START_REF] Thore | [END_REF], most of these geophysical approaches only focus, however, on the uncertainty about specific model elements and do not capture the overall model uncertainty. An alternative to capture this global uncertainty resides in the use of geostatistical approaches, which have been largely developed in the geomodeling community (see [START_REF] Wellmann | 3-D Structural geological models: Concepts, methods, and uncertainties[END_REF] for a review). These approaches also have, however, their own limits: they often require an exhaustive description of each individual source of uncertainty, and generally fail to produce models in accordance with seismic data.

In this paper, we focus particularly on the structural uncertainties related to the presence of salt bodies in seismic images, which are known to be particularly difficult [START_REF] Jones | Seismic imaging in and around salt bodies: Interpretation[END_REF]. The presence of salt often introduces large and abrupt velocity contrasts that impair automated tomographic inversion. This leads to an iterative processing during which the salt boundaries are manually picked one after another, the domain beneath a boundary is flooded with either salt or sediments (depending on the iteration), and the seismic data are remigrated [START_REF] Jones | Estimating subsurface parameter fields for seismic migration: Velocity model building[END_REF]. In this context, accounting for interpretation uncertainties is difficult because building the velocity model requires manual interaction [START_REF] Jones | Seismic imaging in and around salt bodies: Interpretation[END_REF]. Full waveform inversion (FWI) is an appealing alternative and several recent methods have been proposed to automate FWI of salt structures [e.g., ???]. Even so, the combined effects of cycle skipping and physical limitations do not guarantee that the obtained solution corresponds to the global optimum, and the Hessian-based uncertainty measures around that solution are likely to be under-estimated.

To tackle this problem, a relatively uncommon approach is adopted in this paper. Rather than considering the velocity model in a geophysical data-driven way, we consider the velocity models as directly controlled by the subsurface structures. This permits to use stochastic structural modeling to define a set of geologically reasonable velocity models, and to quickly explore the impact of a large range of salt interpretations on seismic imaging. In this context, velocity variations are directly related to the sampled structural uncertainties. Two questions form the primary focus of this paper: Can a stochastic geological velocity modeling framework help to reduce imaging structural uncertainties in a reproducible way and improve a seismic image? How to extract reliable information from multiple, knowingly inaccurate, seismic images?

To answer the first question, we compare the difference in response of two different seismic data sets to a single set of migration velocity models. These two data sets correspond respectively to a "connected" salt diapir and a "detached" salt diapir that is connected to its primary salt layer by a weld. Performing this comparison requires to answer the second question. A common way to assess the quality of a seismic image consists in computing the focusing of the reflections in the image, often based on the measurement of the residual moveout in common image gathers [Biondi, 2006]. Though such approaches are well suited to assess the local quality of an individual seismic image, they do not permit to capture the overall behavior of a seismic image set. [START_REF] Protasov | Subsalt imaging in the presence of salt-body uncertainty[END_REF] propose to define a "mean" seismic image as a weighted function of multiple perturbed seismic images obtained by stochastically sampling salt boundary geometries (the weights are determined using a local focusing measure). Though this permits to handle multiple inaccurate seismic images, it is not applicable to the assessment of structural imaging uncertainties as it focuses on providing a single deterministic image. In order to extract both qualitative and quantitative information from a set of perturbed seismic images, we investigate the use of a pixel-based statistical description of the image set. The idea is to compute, for each pixel, the distribution of its values over the image set and to summarize this information into a comprehensive form. In particular, for a given image set, two new images are generated by respectively taking for each pixel the mean and the standard deviation of its distribution. This is illustrated on two examples. In the first one, the "mean" and "standard deviation" images are computed from the seismic amplitude images. It is very similar, in principle, to stacking common mid-point gathers. Directly working with the seismic amplitude is, however, extremely sensitive to phase changes implied by local velocity changes. To solve this problem, a second example illustrates the use of seismic attributes as a substitute to amplitude.

The remainder of this article is organized as follows. The first section details the geological modeling, the velocity modeling and the seismic migration scheme. The second section illustrates the statistical description of a set of perturbed seismic images. The third section presents a combined application of the two previous sections to the iterative refinement of a salt diapir interpretation. Finally, the last section discusses some elements toward the integration of the proposed methodology into seismic imaging workflows.

C

Overview

The work is carried out on a synthetic case study that is 4 km long and 6 km deep. It is composed of subhorizontal sediment strata, a basal salt layer at the bottom, and a diapir originating from this salt layer. At this stage of the imaging workflow, it is assumed that the quality of the seismic image is relatively poor, and that only the top of the diapir and the top of the salt layer away from the diapir are properly imaged. The main task is, therefore, to elucidate what is happening along the diapir flanks and in particular, to determine whether the diapir is connected to the salt layer or has been detached and is now overlying a salt weld [START_REF] Hudec | Terra infirma: Understanding salt tectonics[END_REF]].

To answer this question, two acquisition data sets (corresponding to a connected diapir and a detached diapir) are generated using seismic modeling. The difference in response of these two data sets to seismic migration is investigated by generating an independent set of 100 migration velocity models which are used to remigrate both data sets. An important facet of the velocity modeling phase is that the velocity is assumed to be directly controlled by the geological structures. Velocity modeling thus includes an initial geological modeling step, to simulate the structures that guide the velocity modeling. In the end, two collections of 100 seismic images (obtained from inaccurate velocity models) have been generated, one for each diapir type. The overall workflow used to generate the data sets is summarized in Figure 1. The remaining of this section details each step of this workflow.

Structural modeling

Modeling strategy

The modeling of salt bodies is performed on a 2D Cartesian grid of 800 by 1200 cells (cells have a dimension of 5 meters by 5 meters) with the method proposed by [START_REF] Clausolles | Generating variable shapes of salt geobodies from seismic images and prior geological knowledge: Interpretation[END_REF] (Figure 2). It is based on an a priori definition of the structural uncertainties, represented as a buffer zone encompassing the salt boundary. This buffer zone, referred to as the Uncertain region in the following, separates two other regions (respectively called Sediments and Salt) where the interpretation of the seismic image is considered unambiguous. The limits of the Uncertain region are manually picked and define a general diapiric shape (Figure 2.a). Within this region, the salt boundary corresponds to the 0 level set of the field D per t,i defined Figure 1 Overview of the creation of the two sets of seismic images. Starting from a manually picked definition of the uncertain salt boundary position, a set of 102 structural interpretations of the boundary and the associated velocity models is stochastically generated. Two of these models (respectively a connected salt diapir (left) and a detached salt diapir (right)) are extracted to create two synthetic seismic data sets by simulating a wave propagation. The remaining velocity models are used to migrate each seismic data set using Reverse Time Migration (RTM). This provides two sets of 100 seismic images obtained with inappropriate seismic velocities that are later compared to the two "true" solutions.

as:

D per t,i = D -ϕ i . ( 1
)
where D is a reference scalar field, and ϕ i is a spatially correlated random field used as a perturbation.

The reference field D (Figure 2.b) is computed from the manual picks defining the Uncertain region: its value is set to 0 at the contact with the Salt region, and to 1 at the contact with the Sediments region. A few additional points are set to 0.5 in the stem of the Uncertain region to ensure an approximately equivalent probability of simulating connected and detached diapirs. The remaining values of D are then interpolated using the finite difference-based structural modeling method proposed by [START_REF] Irakarama | Finite Difference Implicit Modeling of Geological Structures[END_REF].

The random field ϕ i (Figure 2.c) is generated using a sequen- The top of salt boundary is usually well imaged and its interpretation is not prone to much uncertainties [e.g., [START_REF] Mosher | Iterative imaging for subsalt interpretation and model building[END_REF]. As a consequence, the top of salt is considered to have already been interpreted at an earlier imaging stage. Imposing the value of ϕ i to take the value of D at these interpretation points [START_REF] Clausolles | Generating variable shapes of salt geobodies from seismic images and prior geological knowledge: Interpretation[END_REF] freezes the position of the top diapir and of the basal salt layer (Figure 2.a, black dots).

In principle, this perturbation may induce some bubbles in the perturbed geometry, which may be useful to seed allochtonous diapirs in the uncertain zone. The use of a relatively smooth perturbation field ϕ i prevents such occurrences, as does the velocity extension method used by citetyang2019. The produced realizations include models both with connected and detached salt diapirs, as described in the next section.

Assessment of the structural interpretation set

The set of 102 salt structural interpretations contains 54 models of detached salt diapirs and 48 models of connected salt diapirs. Figure 3 illustrates the two models which have been selected for seismic modeling (Figure 3.a, white dotted and dashed lines), and the overall variability of the remaining 100 models which are used for seismic migration.

Two indicators are used to assess the variability of the structural interpretation sets: the posterior probability P sed for each pixel to belong to the sediments (called indicator probability function in [START_REF] Wellmann | Towards incorporating uncertainty of structural data in 3d geological inversion[END_REF] or E-type in Journel [1989]), and the information entropy H [START_REF] Wellmann | Uncertainties have a meaning: Information entropy as a quality measure for 3-D geological models[END_REF]. The P sed indicator (Figure 3.b) directly describes the distribution of the simulated boundaries in the Uncertain region. For any pixel x in the seismic grid, it is defined as [adapted from Wellmann et al., 2010, equation 11]

P sed (x) = 1 100 100 i=1 1 if D per t,i (x) > 0 0 otherwise , (2)
where i identifies a given realization. The information entropy H (Figure 3.c) can be interpreted, in this application, as the predictability of the output category of a pixel (i.e., salt or sediments). It is defined as H(x) = -P sed (x)×log(P sed (x))-(1-P sed (x))×l o g(1-P sed (x)).

(3) The higher is the entropy value, the lower is the predictability of the category. High information entropy values around the central part of the Uncertain region therefore indicate that the structural modeling method samples this model part without favoring specific model geometries. It is interesting to note, however, that locally imposing the position of the salt boundary has a significant impact on the possible salt geometries, due to the relatively large ranges of the variogram model used for simulating the perturbation fields ϕ i .

Velocity modeling

Many ways exist to generate a velocity model, depending on the model purpose and the domain of application [e.g., [START_REF] Dubrule | Geostatistics for seismic data integration in earth models[END_REF]. In this paper, velocity is considered to be directly controlled by the geological structures (another common view being to consider the velocity model as explaining the seismic data): the final velocity model is obtained by overlaying a background sediment velocity model by the salt velocity, where salt is produced by the stochastic structural modeling framework. This approach allows for the fast generation of multiple velocity models from the structural realizations. Following the acoustic approximation, the models only consider the P-wave velocity, which is assumed isotropic. In order to investigate the impact of salt boundary uncertainty, the exact background sediment velocity model V sed is fixed. This choice has two important implications. First, the sediment velocity model is never known in practice and sediment velocity uncertainty should also be investigated [e.g., [START_REF] Grana | Probabilistic inversion of seismic data for reservoir petrophysical characterization: Review and examples[END_REF]. Second, a constant lithology is assumed within salt bodies whereas the effect of unpredicted lithological variations (such as cap rock formation) is known to drastically alter seismic images (e.g., [START_REF] Jones | Seismic imaging in and around salt bodies: Interpretation[END_REF], [START_REF] Jackson | Origin of an anhydrite sheath encircling a salt diapir and implications for the seismic imaging of steep-sided salt structures, Egersund Basin, Northern North Sea[END_REF]). These two points are further discussed in section "Discussion and perspectives".

The background sediment velocity model is generated from a finely layered stratigraphic model S (Figure 4, first column) used to populate a constant per layer velocity model V sed l (Figure 4, second column). The sediment velocity ranges from roughly 1900 to 4100 m.s -1 and globally increases with depth. To reproduce natural variability, V sed l is then perturbed by a field φ (Figure 4, third column) obtained with a Sequential Gaussian Simulation (SGS). The SGS uses a Gaussian distribution model ( (1, 0.01)) and a Gaussian variogram model having relatively large ranges along the stratigraphic horizons (2000 meters), and a short range orthogonal to the layers (100 meters). The resulting finely layered velocity model V sed is defined as (Figure 4, last column)

V sed (x) = V sed l (x) × φ(x). ( 4
)
The presence of welds in the structural models (i.e., when simulating detached diapirs) brings an additional difficulty by introducing a discontinuity in the velocity models, especially since the position of this discontinuity varies from one realization to another. To handle varying weld locations, three different background sediment velocity models (V sed cer t , V sed + and V sed -) are defined and then combined depending on the output of the structural modeling workflow. The V sed cer t model manages the sediment velocity in the Sediments region of the seismic image (where the interpretation is considered certain, cf. Figure 2.a). Each pixel belonging to this region will therefore always have the same velocity from one realization to another (Figure 4, top row). The two other models (the V sed + and V sed models) are complementary and manage the velocity values within the Uncertain region, on each side of the simulated weld or diapir (Figure 4, middle and bottom rows). Details about their generation are provided in Appendix.

The final velocity model is defined as a combination of the background sediment velocity models and the simulated salt boundary(ies), and is obtained by pointwise value transfer (Figure 5). The combination depends on the field D per t,i . Any pixel x of the seismic grid has four possibilities. First, if D per t,i (x) ≤ 0, x belongs to the salt body(ies) and its velocity V (x) is equal to V sal t , which is assumed constant at 4480 m.s -1 . Then, if x belongs to the Sediments region, V (x) is given by V sed cer t (x). If x does not fall into one of these two categories, it means that it corresponds to sediments located in the Uncertain region. The side of the diapir (or weld) to which it belongs is determined by operating a binary partition P + /P -of the field D per t,i . This partition is defined as the subvertical minimum value ridge located in the stem of the Uncertain region. It is determined by applying a watershed transform to the field D per t,i (algorithm proposed by [START_REF] Cousty | Watershed cuts: Minimum spanning forests and the drop of water principle[END_REF], coupled with a local maximum filtering strategy (algorithm proposed by [START_REF] Najman | Building the component tree in quasilinear time[END_REF]). To extract the minimum ridge, the transform and the filtering are applied on the opposite of the field D per t,i . The part of this ridge connecting the salt base and the bulb determines the exact location of the simulated weld (or simply crosses through the diapir in the case of a connected stock) [Clausolles, 2020, chapter 3]. Depending on the part P + or P -it belongs to, x respectively takes the value of V sed + or V sed -. The final velocity model V can thus be written as

V (x) =          V sal t if D per t,i (x) ≤ 0 V sed cer t (x) else if D(x) ≥ 1 V sed + (x) else if x ∈ P + V sed -(x) otherwise . ( 5
)

Seismic modeling

Seismic modeling is performed with a finite difference acoustic wave simulation assuming constant density, using the velocitystress formulation described by [START_REF] Virieux | SH-wave propagation in heterogeneous media: Velocitystress finite-difference method[END_REF]. The implemented code uses 8 th order difference operators in space and 2 nd order difference operators in time, uses perfectly matched layers boundary conditions [START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]] and runs on GPU [START_REF] Micikevicius | 3D finite difference computation on GPUs using CUDA[END_REF].

In order to ensure a proper illumination along the diapir flanks, the velocity model is laterally extended by horizontally extrapolating the first and last traces. The final model geometry is 1200 by 1200 cells, with a spatial grid discretization h of 5 meters (both vertically and horizontally).

The seismic modeling is performed considering 13 seismic sources and 31 receivers, evenly distributed every 500 and 200 meters respectively. The source wavelet is a derivative of Gaussian with a dominant frequency of 25 Hz [START_REF] Virieux | SH-wave propagation in heterogeneous media: Velocitystress finite-difference method[END_REF]. Each shot is slightly delayed (so that the wavelet is centered at 0.2 second) and is recorded at the receivers location for 5 seconds. The time discretization d t used for the wave simulation is determined by d t = 0.2 × h V max ≈ 0.22 millisecond. The result of the simulation is a set of 13 shot gathers. These gathers are resampled with a time step of 2 milliseconds and the first arrivals (corresponding to direct waves) are muted. This is done by defining a cone starting at t 0 = 0.4 second with a slope equal to the offset divided by the minimum speed V 0 = 1900 m.s -1 , and muting any signal outside this cone.

Seismic migration scheme

Due to the complex shape of salt bodies, a pre-stack depth migration method is applied: the Reverse Time Migration [RTM, [START_REF] Baysal | Reverse time migration[END_REF]. The RTM implementation follows the same rules as the wave simulation, but it is performed on a 600 by 600 cell grid to limit the computational cost. The velocity models used for migration are obtained by subsampling the stochastically generated velocity models (cf. section "Velocity modeling") by a factor 2. The time discretization used for the RTM is thus about 0.44 millisecond. A post-processing is applied to the RTM images to normalize the amplitude, which includes a high-pass filtering step and automatic gain control based on the illumination map.

S

The idea underlying the use of multiple stochastically generated seismic images is relatively simple. A seismic image is the result of two contributions: one from the seismic data, and one from the migration velocity model. If one can stack the different images into a single one, then it should be possible to increase the contribution of the seismic data as compared to the individual velocity model contributions. In this section, we investigate two avenues to perform this stacking. Both of them consist in a statistical analysis performed at the pixel scale.

Statistical description of the image set

A way to locally describe the seismic response variability is to characterize -for each pixel -the distribution of amplitude values over all the seismic images. This information is summarized using two quantities: the mean and the standard deviation of this distribution. They are used to generate two images: the mean amplitude image I amp and the amplitude standard deviation image σ amp .

Mean amplitude image

The mean amplitude image I amp of a set of N seismic amplitude images I amp,k , k = 1, . . . , N is defined as

I amp (x) = 1 N N k=1 I amp,k (x) (6)
for any pixel x in the image. Computing I amp from a set of stochastically generated seismic images is very similar, in principle, to stacking common mid-point gathers. In the image parts where all the velocity models have approximately the same response (i.e., where the signal is locally in phase from an image to another), the reflections are preserved. On the contrary, in the image parts having disparate responses (i.e., where the signal is out of phase from an image to another), summing the different contributions tends to cancel the signal.

In the I amp image, the structures which are consistent from one realization to another are therefore emphasized, while the inconsistent ones fade away. Figure 6 illustrates the I amp images computed for the two (connected and detached) diapir data sets. Reflections within the Sediments region are well preserved, which is expected because the velocity models do not vary within this region. Reflections are thus always relocated at the same (exact) position. For the same reason, the top of salt reflections are coherent in the parts where their position was imposed. The main interest is, therefore, to characterize what happens in the Uncertain region and especially in the part between the interpreted top of the diapir and the base salt layer (the stem region, Figure 6, black frame). In the upper part of the diapir (where the modeling method always simulates salt, cf. Figure 3), both images exhibit a very smooth, low amplitude texture, which is consistent with the presence of salt. In the stem region, reflections are attenuated (due to the relatively low illumination and the presence of salt in some of the velocity models). Several differences between both scenarios can, however, be noticed. First, the detached diapir image displays a strong inverse polarity reflection located beneath the upper Salt region. This reflection is typical of a base of salt reflection, and is located at the exact salt boundary position. The correct relocation of the reflection can be explained by the fact that the migration velocity models are correct above the base of the salt bulb. Second, another relatively strong (but partly blurred) reflection appears at the bottom of the detached diapir image. Although it is slightly shifted down, it can be associated to the top of the basal salt layer. The shifting might be explained by the inaccuracy of the velocity models in the stem region: contrarily to the base of salt reflection, this reflection is not always relocated at the same position from a seismic image to another (which also explains the partial blurring of the reflection). Finally, subhorizontal attenuated reflections corresponding to sediments are visible in both images, but they exhibit different patterns. In the connected diapir image, numerous inconsistent cross-cutting reflections can be observed, whereas in the detached diapir image they approximately stop at the weld position. This last observation should, however, be mitigated, as a few cross-cutting reflections also appear in the detached diapir image.

Amplitude standard deviation image

The amplitude standard deviation image σ amp of a set of N seismic images I amp,k , k = 1, . . . , N is defined as

σ amp (x) = 1 N N k=1 I 2 amp,k (x) -I amp 2 (x) (7)
for any pixel x in the image. The analysis of I amp only yields a qualitative information about the average seismic response.

The computation of σ amp offers an additional quantitative information to help interpret I amp . As it relates to the average dispersion of the amplitude values around the mean, this measure can be used to characterize the sensitivity of the data sets to variations in the velocity models. Two (complementary) interpretations can be derived from this sensitivity. The first one, which is the one preferred in this paper, is that this measure can be used as a proxy to characterize the imaging-related structural uncertainties. As it is sensitive to amplitude changes from one image to another, it is well suited to track the vertical reflection shifts induced by changes in the migration velocity models. However, a larger dispersion of the amplitude values also means that seismic data are more sensitive to inaccuracies in the velocity model, and thus that the range of models that correctly focus the seismic energy at depth is narrower. Accordingly, large σ amp values can also be related to lower imaging uncertainties. To summarize, large σ amp convey both a larger structural uncertainty and a lower refocusing uncertainty. Figure 7 illustrates the σ amp images computed for each data set. The first observation is that several spots exhibiting large σ amp values are present in the Uncertain region of both data sets (close to the Sediments region). More precisely, these spots are located at the terminations of the reflections corresponding to the largest model velocity contrasts (thus the largest amplitude values), and around the uncertain salt boundary. An explanation is that for each velocity model, the potential reflections tend to be large in the sedimentary regions, and are significantly lower in salt regions. This results in a different response from one model to another, hence the large σ amp values. In the stem region (Figure 7, black frame), the most interesting element is the presence of two regions of very large σ amp values in the detached diapir image. They correspond to the two (bottom and top) salt reflections observed in the I amp image. These two regions quantitatively confirm that the seismic image content strongly varies depending on the input velocity model, which is consistent with the presence of (un)predicted salt boundaries. Another interesting observation is that, except the several spots described above, the σ amp values are globally very low in the Uncertain region of the connected diapir image. This shows that the low amplitude values observed in the I amp image of the connected diapir are due to globally low amplitude values over all the realizations rather than to the canceling of opposite phase signals, which is consistent with the presence of salt. On the contrary, in the stem region of the detached diapir image, two relatively high value spots can be observed (on the right of the weld). These spots indicate that, even though the reflections are attenuated in the mean seismic image, moderate reflections occur in this part of the model, which would be a sign of structured stratigraphic reflections in this area.

It can also be interesting to jointly visualize the mean seismic image and its variability at the same time, as this is known to affect human perception [START_REF] Viard | Adjacent versus coincident representations of geospatial uncertainty: Which promote better decisions?[END_REF]. Figure 8 illustrates a composite display where I amp is overlaid with the variance image (the squared σ amp ). Using the variance permits to emphasize the image parts exhibiting the highest variations.

Assessment on the use of statistics on raw seismic images

To summarize, the mean amplitude image I amp obtained from a set of stochastic geology-controlled velocity models (1) yields a reliable image in the zones without uncertainty (mainly where the velocity model is imposed), and (2) only preserves recurrent structures in the uncertain image parts. This permits to isolate, for example, the position of the salt boundaries. The information provided by I amp can be completed with the computation of the amplitude standard deviation σ amp , which allows for the distinction between unresolved image parts (e.g., due to low illumination) and zones of erratic reflections. In the presented example, this permits to distinguish between the stem of the connected diapir and the sediments located beneath the bulb of the detached diapir.

These conclusions, however, depend on the hypotheses used to create the data set, and in particular the use of the exact background sediment velocity model. In practice, the background sediment velocity model would also be affected by uncertainties. Using a partially incorrect velocity model would result in a systematic bias in the interpretations. To contemplate the use of stochastic approaches in real situations, it is therefore necessary to also account for sediment velocity uncertainty, which would introduce local shifts of the reflections from one realization to another. As the computation of statistics directly on seismic images is sensitive to such variations, the next section considers the use of seismic attributes to analyze image variability.

Statistical description based on seismic attributes

Seismic attributes are derived from seismic data to emphasize some specific content in the seismic data or images [START_REF] Sheriff | Encyclopedic Dictionary of Applied Geophysics[END_REF]. In this study, we consider attributes that describe the changes in the seismic image texture, which are commonly used to distinguish salt from sediments. Indeed, salt usually appears as relatively smooth, low amplitude bodies without specific internal organization due to its massive and crystalline nature, whereas layered sediments often appear as packages of planar reflections [e.g., [START_REF] Brown | Interpretation of Three-Dimensional Seismic Data[END_REF]. Among all the textural seismic attributes that have been proposed, we selected the gray-level co-occurrence matrix (GLCM) proposed by [START_REF] Haralick | Textural features for image classification[END_REF], a generic descriptor widely used for salt interpretation [e.g. [START_REF] Gao | Volume texture extraction for 3D seismic visualization and interpretation[END_REF][START_REF] Eichkitz | Calculation of grey level co-occurrence matrix-based seismic attributes in three dimensions[END_REF][START_REF] Berthelot | Texture attributes for detection of salt[END_REF].

Gray-level co-occurrence matrix and derived attributes

Computing the GLCM requires first to quantize the seismic amplitude into a set of n g discrete gray levels. In this application, 32 gray levels were used, regularly distributed between the minimum and maximum amplitudes of each individual seismic image. In its general formulation, the GLCM is a two-point categorical statistics method describing for two neighboring pixels separated by a fixed lag vector the probability of each possible pair of gray levels. As commonly done in seismic interpretation [START_REF] Eichkitz | Calculation of grey level co-occurrence matrix-based seismic attributes in three dimensions[END_REF], we use a simplified formulation which only considers adjacent pixels and ignores the direction. For a given pixel p, a single GLCM M is computed by scanning all the pairs of contiguous pixels p 2 to each pixel p 1 in the neighborhood N global of p (defined as a window of fixed dimension centered at p). Let M i j be the (i, j) entry of M , where i and j are two gray levels:

M i j (p) = p 1 ∈N global (p) p 2 ∈N local (p 1 )
1 if p 1 = i and p 2 = j 0 otherwise , (8) where N local (p 1 ) contains all the pixels contiguous to p 1 (i.e., sharing an edge or corner). The normalized GLCM coefficients m i j are defined as

m i j (p) = M i j (p) n g k=1 n g l=1 M kl (p) , ( 9
)
where n g is the total number of gray levels. These normalized GLCM coefficients describe the joint distribution of the gray level associations in N global (p). The neighborhood is here set as a window of 21 by 21 cells (i.e., 10 cells on each side of the pixel). Numerous attributes can be derived from the GLCM. They can be organized into three categories depending on the characteristics they emphasize [START_REF] Eichkitz | Calculation of grey level co-occurrence matrix-based seismic attributes in three dimensions[END_REF]. The contrast category highlights the strong gray level variations between neighboring pixels. The orderliness category highlights the uniformity of the gray level associations (the opposite being the presence of predominant associations). The statistics category characterizes the gray level distribution in the neighborhood (independently of the pair associations). We selected one attribute of each category to cover the different types of obser-vation the GLCM provides: the contrast, the energy and the standard deviation.

For each seismic image, the GLCM and its three derived attributes are computed at each pixel location. This yields three sets of attribute images, that are processed exactly as the amplitude images in the previous section. For each set, we compute for each pixel the attribute value distribution and summarize it using its mean and standard deviation.

GLCM standard deviation

The GLCM standard deviation M st d is defined as

M st d =   n g i=1 n g j=1 (i -µ) 2 m i j   1 2 , ( 10
)
where µ = n g i=1 n g j=1 i × m i j is the mean gray level in N global . This statistical measure relates to the gray level distribution in the pixel neighborhood (without considering the gray level associations). It describes the dispersion of the gray level values, corresponding here to the local dispersion of the seismic amplitude in each pixel's neighborhood.

Figure 9 illustrates the GLCM standard deviation mean and its SD images computed for each data set. It shows that the very high values are localized along strong amplitude reflections in both data sets (top and base of salt and strong sediment reflection packages). As previously, the base of the salt bulb and the top of the base salt layer are well visible in the detached diapir data set. On the contrary, image regions corresponding to low amplitude variations (due to the presence of salt) appear with low M st d values. Another major difference is the connected diapir shape of the low M st d values in the first data set, whereas only a bulb and a base layer of low values are visible in the second data set.

The standard deviation images complete the information provided by the mean images. First, as the strong sediment package reflections are recurrent from one image to another, they only present a low dispersion, and thus low standard deviation values. This permits to distinguish between the strong sediment package reflections and the subhorizontal salt boundaries. Second, as already observed on the σ amp images, the presence of two moderate values spots within the stem region of the detached diapir data set indicates the existence of reflections beneath the diapir.

GLCM contrast

The GLCM contrast M cont r ast is defined as

M cont r ast = n g i=1 n g j=1 (i -j) 2 m i j . ( 11
)
It is designed to emphasize the presence of large gray level variations in the neighborhood of a pixel. As compared to the GLCM standard deviation, it emphasizes the largest gray level variations while minimizing the smaller ones, so it yields sharper images. The GLCM contrast is often used for salt interpretation as it captures well the salt boundaries marked by strong amplitude reflections [e.g., [START_REF] Gao | Volume texture extraction for 3D seismic visualization and interpretation[END_REF][START_REF] Eichkitz | Calculation of grey level co-occurrence matrix-based seismic attributes in three dimensions[END_REF][START_REF] Berthelot | Texture attributes for detection of salt[END_REF].

Figure 10 illustrates the GLCM contrast mean and SD images computed for each data set. The low amplitude regions (which are due to salt) are even more visible in the contrast mean images than in the GLCM standard deviation mean images. The top of the diapir is also clearly visible, together with the base of the bulb in the detached diapir data set. The top of the base salt layer does not appear, however, as clearly as previously (only the highest GLCM standard deviation value spots are preserved). The highest reflection sediment packages are also well visible.

The standard deviation image of the GLCM contrast permits to distinguish between these packages and the regions where the amplitude significantly varies depending on the migration velocity model. Indeed, recurrent high amplitude sediment packages have low contrast variations from an image to another, while the contrast variations in variable amplitude regions (e.g., around the uncertain salt boundaries) are much more significant. In particular, even though the top of the base salt layer does not clearly appears on the mean contrast image of the detached diapir data set, the standard deviation image records significant variations at its location.

GLCM energy

The GLCM energy M ener g y is defined as:

M ener g y =   n g i=1 n g j=1 m 2 i j   1 2 . ( 12
)
It emphasizes the presence of preferential gray level associations in the neighborhood of a pixel. It is regularly used (or some variant) for interpreting salt [e.g., [START_REF] Eichkitz | Calculation of grey level co-occurrence matrix-based seismic attributes in three dimensions[END_REF].

As amplitude is generally low within salt, salt bodies only contain few different gray levels. Therefore, only a few different gray levels associations are found in salt, which results in large energy values. On the contrary, regions with larger amplitude ranges (such as sediments) have more possible gray level associations and thus exhibit lower energy values. We chose the energy rather than other common representatives of its GLCM attribute category (namely the angular second moment and the entropy) because it formulation varies "linearly" with the m i j coefficients, contrarily to the angular second moment (which is the squared energy) and the entropy (which is a function of m i j and its logarithm).

Figure 11 illustrates the GLCM energy mean and SD images computed for each data set. As for the GLCM contrast, salt bodies can be clearly identified (i.e., a connected diapir, a salt bulb, and a base salt layer). Also, the GLCM energy values are significantly higher in the stem region of the connected diapir data set than in the detached diapir data set. This second observation shows that the reflections are more attenuated in this image part in the connected diapir data set (which is consistent with the presence of salt).

The standard deviation images are more difficult to interpret. The higher variations that are present in the stem region of the detached diapir data set seem consistent with the last observation. Higher variations indicate that the attenuation varies more in this region depending on the migration velocity model, which is more likely to happen if there is no salt (as in this case, the attenuation is mainly controlled by the prediction of the migration velocity model).

Assessment on the use of statistics on GLCM attributes

The statistical analysis of the GLCM attributes provides results similar to those obtained from the direct computation of statistics on the amplitude images. It allows to distinguish between the two types of diapirs, based on low values of GLCM contrast and standard deviation and high values of GLCM energy. Furthermore, for the detached diapir, the presence of high value spots in the contrast and standard deviation mean images and of low value spots in the energy mean image underlines the presence of two salt boundaries (the base of the salt bulb and the top of the base salt layer). The strong packages of sediment reflections are also emphasized, in a very similar way as the salt boundaries. The computation of the standard deviation images associated with the different GLCM attributes permits, however, to distinguish between these sediment and salt reflections (still under the assumption of a fixed background sediment velocity model). The only information which is not recovered as compared to the computation of statistics directly from the amplitude images is about the polarity of the salt reflection (which permits to distinguish between top and base of salt reflections). As a conclusion, the use of GLCM attributes (and more generally of seismic attributes) to characterize the variability of a stochastically generated seismic image set is a good alternative to the computation of mean and standard deviation images, and has the advantage of being more robust to local reflection shifts.

The information extracted from the different selected GLCM attributes is somehow redundant. Indeed, the GLCM attributes mainly measure two types of information (the difference between gray levels, i.e., the i and j pairs, and the distribution of the gray level associations, i.e., the m i j coefficients) with different manners of balancing them (cf. equations 10, 11 and 12). Nevertheless, we think that it is important to consider each attribute (and in general many attributes) to extract as much information as possible from the seismic images (cf. section "Toward automatic velocity model updating").

A

In this section, we present a quick application combining the generation of velocity models from stochastic geological scenarios and the statistical analysis of seismic images. Two main questions are investigated: Is it possible to discriminate between the two types of diapir end-members (connected and detached) without a priori information? Up to which point can an interpretation be refined using a stochastic structural modeling approach?

Description of the application

The application is performed using the two seismic data sets from Figure 1 (that correspond respectively to a connected and a detached diapir), and with the same a priori definition of the uncertain buffer zone (illustrated in Figure 2). Three scenarios (i.e., sets of structural modeling constraints) are successively considered to generate velocity models, and the resulting seismic images are then analyzed. Except for the last scenario, the same sets of velocity models are used to migrate both seismic data sets.

Case 1: No constraints

First, no conditions are imposed on the position of the salt boundary. A single set of one hundred velocity models is generated by simulating salt envelopes in the buffer zone, and each of the velocity models is used to migrate both seismic data sets. Finally, the I amp and σ amp images are computed for each seismic data set (as described above).

These images are illustrated in Figure 12 (first row). In both I amp images, the main observation is the presence of a strong reflection, corresponding to the top of the salt boundary, that is correctly relocated. These images do not provide much more information, as they only exhibit inconsistent, rather blurred, cross-cutting reflections beneath the top of salt reflection. The σ amp images provide more insights. First, at the position of the top of salt reflection, we observe in both images a large spot of high standard deviation values that is much thicker than the top of salt reflection. These spots show that, from one realization to another, the reflection is relocated at a noticeably different depth even though it is properly captured on I amp images. Second, even though no clear differences can be made in the lower part of the I amp images, the σ amp images exhibit a different behavior, with the presence of two additional large spots of high values in the detached diapir image. The next case investigates this difference in behavior.

Case 2: Constraining the top of the salt

As the case 1 demonstrated that the top of salt boundary can be imaged, the case 2 imposes its position when simulating the velocity models. Again, the same set of velocity models is used to image both seismic data sets. The resulting images are illustrated in Figure 12 (second row).

This case was used to illustrate the section about the statistical analysis of seismic images, so we only summarize here the main conclusions. When fixing the top of the salt boundary, the observed behavior of the I amp images differs. For the connected diapir data set, no more reflections are observed in the diapir stem region. For the detached diapir, however, the base of salt reflection of the detached diapir bulb is correctly relocated, and the top of salt reflection of the diapir base starts to appear. The σ amp images confirm the observations of the first case, indicating two different behaviors between the two data sets.

Case 3: Constraining the top of the salt and the topology of the diapir

As the case 2 permits to distinguish between the two types of diapir end-members, a constraint about the topology of the diapir (i.e., connected or detached) is introduced when simulating the velocity models. For the connected diapir data set, only velocity models corresponding to connected diapirs are generated. For the detached diapir data set, a single constraining point is added at the location of the base of salt reflection that was observed in the second case. Using this point in the stochastic structural and velocity modeling is sufficient to enforce the generation of detached diapir models. The resulting average and standard deviation images are shown in Figure 12 (third row).

In the connected diapir I amp image, the center of the diapir stem clearly appears and the only uncertainty remaining is about the exact location of the stem boundary, as it is slightly blurred. In the detached diapir I amp image, the bulb can be identified and the reflection corresponding to the top of the diapir base layer is recovered. Even the weld is pretty well recovered.

The last row of Figure 12 illustrates the reference seismic images for both seismic data sets (i.e. the image obtained using the same velocity model for seismic modeling and migration). By comparing these reference images to the I amp images obtained from stochastically generated structural velocity models, we observe that these mean images captured correctly the diapir boundaries even with a relatively low number of interpretive constraints. The only feature which is not recovered correctly is the sharpness of the reflection terminations against the salt bodies.

D

We have presented a method for assessing structural uncertainties during seismic imaging through the use of structurallycontrolled stochastic velocity modeling, and applied it to the imaging of salt diapirs. To go further, two broad types of applications can be considered for this method. The first one is a classic structural uncertainty assessment framework after geophysical imaging [START_REF] Wellmann | 3-D Structural geological models: Concepts, methods, and uncertainties[END_REF]]. The second one, which is not developed yet and would be worth investigating, consists in re-injecting the information obtained about the structural sensitivity of seismic migration to velocity variations into the seismic imaging workflow. Indeed, one of the major challenges faced by common velocity model building methods is that, for a variety of reasons, salt cannot usually be directly handled by tomographic inversion and therefore requires to be manually picked [START_REF] Jones | Estimating subsurface parameter fields for seismic migration: Velocity model building[END_REF]Davison, 2014, Jones, 2015]. This section discusses different prospects that would be necessary steps toward real case applications.

Integration of all the sources of structural uncertainty into velocity modeling

The presented application aims at testing the concept of using structure-guided velocities, so only the salt boundary was variable, as this is in general the main controlling factor for imaging. In practice, the other subsurface structures have uncertainties that should be considered as they affect the velocity model accuracy.

Sources of uncertainty to consider

Among the different sources of structural uncertainties, the main ones are the stratigraphic layer geometry, the stratigraphic layer topology (i.e., faults, unconformities, etc.) and, here, the potential salt internal heterogeneity. Furthermore, the velocity distribution within the sediment layers should also be investigated, especially in the absence of well velocity measurements.

Many approaches to perform such investigations have been proposed for sampling multiple geometries and topologies of horizons and fault networks [e.g., [START_REF] Lecour | Modelling of stochastic faults and fault networks in a structural uncertainty study[END_REF], Holden et al., 2003[START_REF] Godefroy | Multi-scenario interpretations from sparse fault evidence using graph theory and geological rules[END_REF]. Horizon and fault uncertainty could be integrated by perturbing the stratigraphic grids which are generated from the horizon picks and hold the background sediment velocity models [e.g. [START_REF] Mallet | Solid earth modeling and geometric uncertainties: Presented at the SPE annual technical conference and exhibition[END_REF]. It would require, however, to parametrize the uncertainty associated to the different structural model elements, which is a recurrent concern with stochastic modeling methods [START_REF] Thore | [END_REF]. Furthermore, it might introduce some technical difficulty to ensure that the different stratigraphic grids (S cer t , S + and S -) are consistent at the contact between the Sediments and Uncertain regions.

We are not aware of methodologies to introduce heterogeneities within the simulated salt bodies in the geomodeling literature. However, several first order solutions can be considered depending on the type of heterogeneity to introduce. For example, the generation of a cap rock [such as an anhydrite sheath, see e.g. [START_REF] Jackson | Origin of an anhydrite sheath encircling a salt diapir and implications for the seismic imaging of steep-sided salt structures, Egersund Basin, Northern North Sea[END_REF] at the top of the diapir can be done by simulating a spatially correlated random field along the salt boundary to define the thickness of this cap rock formation, and then to assign a specific velocity to this cap rock region. Similarly, the presence of sediment inclusions within salt could be simulated using a scalar field and a threshold, to create blobs within the salt bodies.

Impact on the statistical analysis results

Producing a set of stochastic velocity models with varying sediment velocity could potentially affect the results of the statistical analysis. As mentioned in the conclusion of section "Statistical description of the image set", the local velocity changes would introduce reflection shifts which prevent from directly performing the statistical analysis of the seismic image set. This has motivated the use of seismic attributes, presented in section "Statistical description based on seismic attributes", to work with measures that are less sensitive to velocity variations.

Nevertheless, it is not conceivable to obtain a reliable result in the entire model with a single iteration of the proposed method. Indeed, a given reflection in the subsurface can generally be relocated correctly in the subsurface only if the migration velocity model is correct above it (unless inaccuracies in the velocity model cancel out). If this is not the case, its position will depend on how the velocity inaccuracies affect the wavefield. Given that stochastic models were used, it is probable that the lower parts of the seismic images will suffer from the inaccuracies of the upper parts, and that the amplitude variations will be too erratic to extract any reliable content. One such example is illustrated in Figure 12 (first row): in the detached diapir mean amplitude image, the top of the base salt layer is not recovered due to the large inaccuracies that affect the top of the salt bulb.

This leads us to suggest a top-down velocity model building approach (at least for salt modeling, as salt introduces large velocity contrasts). This is very similar to current practices in the industry [e.g., [START_REF] Mosher | Iterative imaging for subsalt interpretation and model building[END_REF][START_REF] Jones | Estimating subsurface parameter fields for seismic migration: Velocity model building[END_REF]. A first iteration could be used to fix the top of salt, a second one to fix the bottom of salt, etc. After the main salt boundaries have been identified (and the topology of the salt body(ies) is known), further iterations can focus on refining the feebly illuminated salt boundaries and the sediment velocity. One could also consider to turn to more classic methods, such as migration velocity analysis, to locally improve the sediment velocity model [e.g., Biondi, 2006].

Toward automatic velocity model updating

One of the current bottlenecks for automatically building (and especially updating) a velocity model in salt tectonics contexts is that salt bodies often introduce large velocity contrasts over too short scales for being handled by tomographic inversion, and therefore require to be manually picked [START_REF] Jones | Estimating subsurface parameter fields for seismic migration: Velocity model building[END_REF]. The proposed method permits to automatically generate multiple velocity models, and can therefore be considered as a potential step toward automatic velocity model updating. The main question to consider now is: How can it been integrated into the imaging loop?

A first solution might consist in using the salt modeling method only as a surface reconstruction tool, to quickly build a salt boundary surface. This is, however, not really interesting as manual picking from a seismic image would probably do better than one stochastic realization. The strength of the method is that it permits to explore a large range of scenarios. The question of the updating can thus be rephrased as: How can we use the information obtained from multiple seismic images to "update" the velocity model? And what does "update" mean when dealing with a collection of velocity models instead of a single deterministic one?

A possible way to address these questions could be to use the stochastic velocity models as starting models for full waveform inversion methods [e.g., ???]. The proposed level-set representation of the salt body could also be integrated within the FWI iterations to extend the approach of ? to the presence of salt welds.

Updating of the Uncertain region

Applying the methodology in an iterative seismic imaging workflow would call for shrinking the Uncertain region at each step. The statistical analysis of the image set can provide insights to perform this update. The geological information that it provides (such as about the nature of the imaged diapir) can also be used to constrain the topology of the simulated salt boundary interpretations.

One of the main limitations of the method is, however, that performing multiple seismic migrations implies intensive computations. Therefore, the number of images should remain relatively small for the method to be practical. The use of multiple seismic attributes can help with this task. As already stated, a large amount of attributes have been proposed in the literature and have proven to help interpret salt. Each attribute being designed to extract a specific content out of seismic images, the cross analysis of multiple seismic attributes can help to improve the accuracy of the observations (the same idea motivates the use of automatic multi-attribute classification). As the computational cost of attribute computation is negligible as compared to that of RTM, this can be a good way to make the best possible use of a limited number of seismic images.

Another (complementary) solution could also consist in using automatic salt interpretation methods (see e.g. [START_REF] Wang | Successful leveraging of image processing and machine learning in seismic structural interpretation: A review[END_REF] and [START_REF] Clausolles | Stochastic seismic interpretation of salt bodies: detection, sampling and impact on seismic imaging[END_REF], chapter 1, for partial reviews) to separately interpret each individual perturbed seismic image, and to compute facies probabilities from these interpretations to isolate the most uncertain image parts.

Updating strategies

As already mentioned, the use of stochastic velocity modeling requires to use a top-down iterative approach for building the velocity model. Optimizing the sampling of velocity models is important to limit the number of iterations. At each iteration, one should therefore focus on specific objectives. An example of strategy for the first iteration(s) can be to try to determine the bulk salt volume. This can be done by simulating salt bodies with a very variable volume (see [START_REF] Clausolles | Generating variable shapes of salt geobodies from seismic images and prior geological knowledge: Interpretation[END_REF] for details about salt modeling parameter selection). A large range of salt volumes can permit to quickly determine the "breaking point" after which the salt volume is too large and systematically deteriorates the seismic images. On the contrary, for the next iterations, it is mandatory to locally constrain the salt boundaries to ensure a proper imaging below them.

Links between stochastic velocity modeling and geophysics

The proposed method can be seen as a geological approach to a geophysical problem (i.e., stochastic structural modeling to build velocity models). In order to concretely move toward automatic velocity model updating, this approach should integrate usual geophysical criteria of velocity model validation [Biondi, 2006]. Furthermore, geophysical knowledge can help to improve the definition of the different input parameters.

Building the background velocity models A first question relates to the way to populate the structural models with velocities. A solution is to use classic velocity estimation methods such as using picked time horizons, stacking velocities, interval velocities, etc. [START_REF] Jones | Estimating subsurface parameter fields for seismic migration: Velocity model building[END_REF]. Another one is to integrate velocity uncertainty in addition to structural variability using, for example, value ranges rather than fixed values during velocity modeling [e.g., [START_REF] Fomel | Structural uncertainty of time-migrated seismic images[END_REF].

Links between velocity models and seismic data From a geophysical point of view, the major limitation to the use of stochastic velocity modeling is that the simulated models may not always be consistent with the seismic data. However, exploring multiple scenarios is potentially essential to address the ill-posedness of the full-waveform inversion problem. Saltrelated uncertainty cannot be investigated when they are manually picked as in current practice. One shall consider, however, ways to make stochastically generated velocity (more) consistent with seismic data. Another perspective to alleviate the computational cost of this approach could be to directly compare the recorded seismic data with synthetic seismic data obtained on each candidate velocity model, see [START_REF] Irakarama | Appraising structural interpretations using seismic data -Theoretical elements[END_REF]. The idea of dehomogenization recently proposed by [START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-D numerical illustrations[END_REF] opens interesting avenues to achieve this goal.

C

We propose to use geology-controlled velocity models to assess structural uncertainties during seismic imaging. The use of stochastic structural modeling as a base for velocity modeling permits to quickly generate a large range of significantly different velocity models to sample the structural imaging uncertainty space. We couple it with a pixel-based statistical analysis of the seismic images obtained after remigration in order to characterize imaging structural uncertainties.

The application of the method on a 2D synthetic example of salt diapir shows that, at least in almost ideal conditions, it helps to extract reliable information about the structures at depth. Particularly, it permits to distinguish between different diapir topologies (i.e., connected or detached) and to identify the seismic image parts which are most sensitive to migration velocity variations, and thus locating most of the structural uncertainties.

This method is a promising step to move toward automatic velocity uncertainty in salt tectonics contexts. Although seismic attributes likely address uncertainties in the sediment velocities, these results should be confirmed by additional studies considering this possible source of inaccuracy. Another important avenue relates to further consider geophysics-based methods (velocity analysis and full-waveform inversion) to constrain geologically-controlled velocity models.

A

This work was performed in the frame of the RING project (https://ring.georessources.univ-lorraine.fr/) at Université de Lorraine. We would like to thank for their support the industrial and academic sponsors of the RING Consortium managed by ASGA. The software corresponding to this paper is available as the Goscope plugin of the SKUA-GOCAD geomodeling software suite and the SIGMA software suite. We also acknowledge AspenTech for the SKUA-GOCAD software and development kit. Finally, we thank two anonymous reviewers and the Associate Editor Marcilio Castro de Matos for their constructive comments.

tial

  Gaussian simulation (SGS). It is controlled by two parameters: a distribution model (which describes how much the salt boundary can move away from its reference position) and a variogram model (which controls the spatial correlation, and thus the local shape of the boundary). Here the SGS uses a triangular distribution model (min = 0, mode = 0.5, max = 1), and an isotropic Gaussian variogram model with principal ranges of 2000 meters.

Figure 2

 2 Figure 2 Generation of the salt structural interpretations. (a) Input manual picks used to define the Uncertain region (dots colored with interpolation constraint values) and to locally impose the position of the salt boundary (black dots). (b) Reference field D interpolated from the Uncertain region picks. (c) Perturbation field ϕ i obtained by conditional sequential Gaussian simulation (SGS) using a triangular distribution model and an isotropic Gaussian variogram model. (d) Corresponding field D per t,i (equation 1) and simulated salt boundary (dashed curve).

Figure 3

 3 Figure 3 Variability of the simulated salt boundary interpretations. (a) Reference field D and the two interpretations which have been selected for seismic modeling (white dotted and dashed lines). (b) Posterior sediment probability P sed (equation 2). (c) Information entropy H (equation 3). The white parts of the Uncertain region in (b) and (c) indicate that the facies at the pixel location is constant over the 100 realizations.

Figure 4

 4 Figure 4 Building of the background sediment velocity models. Top: models used in the Sediments region. Middle and bottom: models used in the Uncertain region, on each side of the weld. Left: stratigraphic models. Middle left: constant per layer velocity models. Middle right: velocity perturbation. Right: finely layered velocity models (equation 4).

Figure 5

 5 Figure 5 Stochastic velocity modeling. The velocity model associated to a given scalar field D per t,i is obtained by overlaying an adaptive background velocity model representing the sediments with the salt velocity according to the simulated salt body described by D per t,i . The adaptive background model is the combination of three sediment velocity models, governing the different parts of the seismic grid (equation 4, Figure 4): V sed cer t is used in the Sediments region, V sed + and V sed are used in the Uncertain region on each side of the simulated weld or stock (equation 5). The dashed line separating the V sed + and V sed models indicates the limits of the P + and P -partitions defined by the watershed transform of the field D per t,i . The two illustrated velocity models are those selected for seismic modeling (respectively corresponding to the connected diapir data set (left) and the detached diapir data set (right)).

Figure 6

 6 Figure 6 Mean amplitude images I amp (equation 6) corresponding to the connected diapir (left) and the detached diapir (right) data sets. The black dashed lines are the Uncertain region boundaries. The white dashed lines are the reference salt boundaries used for seismic modeling. See text for details.

Figure 7

 7 Figure 7 Amplitude standard deviation images σ amp (equation 7) corresponding to the connected diapir (left) and the detached diapir (right) data sets. The black dashed lines are the Uncertain region boundaries. The white dashed lines are the reference salt boundaries used for seismic modeling. See text for details.

Figure 8

 8 Figure 8 Variance images overlaying mean amplitude images I amp for the connected diapir (left) and the detached diapir (right) data sets. This display permits to visualize the image parts which are the most sensitive to variations in the migration velocity model. The dashed black lines are the uncertainty envelope boundaries. The dashed white lines are the reference salt boundaries used for seismic modeling.

Figure 9

 9 Figure 9 GLCM standard deviation mean and SD images for the connected diapir (left) and the detached diapir (right) data sets. See text for details.

Figure 10

 10 Figure 10 GLCM contrast mean and SD images for the connected diapir (left) and the detached diapir (right) data sets. See text for details.

Figure 11

 11 Figure 11 GLCM energy mean and SD images for the connected diapir (left) and the detached diapir (right) data sets. See text for details.

Figure 12 I

 12 Figure12I amp and σ amp images obtained using different sets of structural constraints when simulating the structural velocity models, respectively: no constraints (first row), top of salt constraints (second row), top of salt and topological constraints (third row). The "reference" seismic images (obtained by using the same model for seimic modeling and migration) are illustrated for comparison (last row).

APPENDIX A: P

To handle varying weld locations, three different background sediment velocity models (V sed cer t , V sed + and V sed -) are generated and combined depending on the output of the structural modeling workflow. The following details the practical aspects of their generation.

The manually picked horizon data (Figure 4, S cer t , black lines) are locally pulled up in the central part of the image to reproduce the upturned geometry of sedimentary strata associated to diapir downbuilding. The resulting stratigraphic grid S cer t is used to build the constant per layer velocity model V sed l,cer t and its perturbation φ cer t , that are combined using equation 4 to obtain V sed cer t . It is crucial for the velocity values to be continuous at the contact between the Sediments and Uncertain regions. To ensure this continuity, elements of the S cer t , V sed l,cer t and φ cer t models are used to build V sed + and V sed -(Figure 4, S + and S -, black lines -for horizon data-and black dotted lines -for new picked data). The geometry of the horizons is here again upturned. To handle the contact between these horizons and the S cer t data constraining the upper model part, an unconformity in the S + and S -models is introduced (Figure 4, S + and S -, thick red line). The position of the unconformity is chosen such that it is not visible in the final velocity model: the strata truncations only occur within the simulated salt bodies and in the V sed + and V sed model parts that are not used (i.e., on the opposite side of the diapir or weld, as described in the next section). As stratigraphic continuity is ensured between the different models, using the layer velocities of V sed l,cer t for V sed l,+ and V sed l,-guarantees the global consistency. The last condition to fulfill is that the φ + and φ -perturbations also have to be continuous with φ cer t at the contact between the Sediments and Uncertain regions. This is ensured by using some of the φ cer t values as constraining data when performing the SGS (Figure 4, φ + and φ -, whitened areas). The V sed + and V sed are finally obtained using Equation 4, reminded here: