Study on the limitations of travel-time inversion applied to sub-basalt imaging

Abstract. The difficulties of seismic imaging beneath high velocity structures are widely recognised. In this setting, theoretical analysis of synthetic wide-angle seismic reflection data indicates that velocity models are not well constrained. A two-dimensional velocity model was built to simulate a simplified structural geometry given by a basaltic wedge placed within a sedimentary sequence. This model reproduces the geological setting in areas of special interest for the oil industry as the Faroe-Shetland Basin. A wide-angle synthetic dataset was calculated on this model using an elastic finite difference scheme. This dataset provided travel times for tomographic inversions. Results show that the original model can not be completely resolved without considering additional information. The resolution of nonlinear inversions lacks a functional mathematical relationship, therefore, statistical approaches are required. Stochastic tests based on Metropolis techniques support the need of additional information to properly resolve sub-basalt structures.


Introduction
Sub-salt and sub-basalt imaging has been a key objective during the last 2 decades for the oil exploration (Rousseau et al., 2003;Williamson, 2003;Sava and Biondi, 2004).Oil exploration has revealed the imaging difficulties in the presence of high velocity features (such as salt and/or basalts).Low velocity structures under relatively high velocity features are poorly constrained by conventional processing and/or inversion schemes (Flecha et al., 2004).Velocity provides the link between seismic images and rock types.Ray trac-ing theory, based on Fermat's principle, states that regions surrounded by higher velocities are under-sampled by rays.Seismic images of the subsurface strongly benefit from well resolved estimation of seismic velocities.These seismic velocities are currently determined by velocity analysis and, in the best case, by travel time tomography (or by the inversion of travel time of first arrivals) of wide-angle seismic reflection/refraction shot-gathers (Zelt and Smith, 1992).The determination of the velocity models requires the interpretation/identification of the seismic arrivals within a shotgather.Furthermore, the mathematical inversion schemes require digitized travel times, offset pairs, to calculate velocities.Usually, a standard crustal velocity model features an increasing velocity with depth, however in the presence of salt and/or basaltic intrusions this assumption fails.Intrusions often represent the emplacement of a high velocity body in the crust, therefore zones beneath these structures may feature low velocities.This is the case of basalt covered areas as erupted basalt buries previous structures that may feature low velocities such as in the Faroe Shelf.In the Faroe Shelf, covered areas represent potential hydrocarbon reservoirs, therefore this topic is of special interest for the industry.
This manuscript develops a theoretical study investigating the reliability and effectiveness of travel time inversion of wide-angle seismic reflection/refraction data.We investigated when and to what extent the low velocity structure beneath a high velocity basalt layer can be resolved by using travel times alone.This manuscript is not aimed to present an specific inversion scheme, the validity and calibration of the inversion scheme used here is presented elsewhere (Trinks et al., 2005).Also relatively recent new developments which are more computationally expensive such as full waveform inversion could be employed to constraint sub-basalt features, however, these approaches are beyond the scope of this contribution.

Geological setting and imaging problems
In the Faroe-Shetland Basin, Mesozoic and Tertiary sedimentary sequences fill the basin but, close to the Faroe Shelf, these sedimentary sequences are covered by Paleocene-Eocene basaltic lavas of which the Faroe Islands are composed.The previous topography of the basin was dominated by normal faults as a consequence of the extension and subsidence during Cretaceous and Paleocene (Richardson et al., 1999).As huge amounts of molten rock were extruded, after filling the lows between fault blocks, lava flows extended over long distances in the basin.Basalt flows were erupted in several episodes and three major units have been identified: Lower, Middle and Upper Series.The composition and thickness differs from one unit to another.Moreover, in periods without igneous activity, lacustrine shales and coals were accumulated and sediments were emplaced filling the basin floor deeps (White et al., 2003).The resulting structure in the Faroe-Shetland Basin may be considered as a relatively thin wedge/finger of basaltic rocks emplaced within a sedimentary sequence This structure developed within the tectonic framework of the evolution of the North Atlantic Igneous Province and the opening of the NE Atlantic rift (Jolley and Bell, 2002).
In the Faroe Shelf, geologic and geophysical data suggest that a layer of basalt is placed within two low velocity sedimentary sequences (Hughes et al., 1998;Richardson et al., 1998Richardson et al., , 1999;;Fliedner and White, 2003;Smallwood et al., 2001;Sørensen, 2003;White et al., 2003;Raum et al., 2005).The velocity structure for sediments above the basalt can be resolved by conventional techniques.The top of the basalt layer can be determined very effectively due to the high contrast in seismic physical properties between the basalt and the overlying sediments.However, the high velocity basalt layer represents a complex scenario for seismic imaging methodologies, acting as a barrier so that the underlying structures can not be imaged.The high velocity that characterises the basalt contrasts with relatively low velocity of the surrounding materials.This causes that most of energy is reflected and/or travels along this layer.Furthermore, the heterogeneous structure of the basalt layers scatter 100 the higher seismic frequencies of the source signal (Pujol and Smithson, 1991;Hobbs, 2002).The lack of penetration and the multiple scattering within the basalt layer obscures the potential seismic events generated bellow the basalt.This could represent potentially prospective sedimentary structures.
Although basalt flows tend to be sub-horizontal on large scale, at small scale, rugged interfaces cause scattering and disperse the elastic energy destroying any lateral coherency of possible sub-basalt events.In addition to the differences between the three major Series, within every unit, basaltic bodies are highly heterogeneous in composition and physical properties.These heterogeneities strongly disperse the seismic energy and destroy the signal coherence in the seismic wave-field (Pujol and Smithson, 1991).The outer parts of the basalt flows are affected by weathering causing a decrease in velocity, this contrasts with the internal parts which cooled slowly and without any external influence preserving a high velocity feature.Interfaces between individual flows in a basalt block produce internal multiples and wave conversions.Also some intrusive basalt flows were emplaced as sills within previous structures providing an additional cause for scattering at and beneath the base of the basalt.
In addition to these major imaging issues, the usual problems of marine seismic reflection data acquisition must be also considered (tidal noise, multiples, peg-leg, reverberation, converted-waves, etc.).In the Faroe Shelf, conventional seismic reflection techniques are insufficient to study subbasalt structures.Sub-basalt imaging is very sensitive to acquisition and processing parameters.In acquisition, long offset 2-D and 3-D seismic data can contribute to an improvement in the seismic image below top basalt.Increasing the source energy at low frequencies by: towing the source and receiver cables at deeper levels and/or using bubble-tuned rather than conventional peak-tuned source arrays.Further improvement can be provided by High frequencies (dominantly noise) are filtered out of the data early in the processing to concentrate on the low frequency data.Careful multiple removal is important with several passes of de-multiple being applied to the data using both Surface-Related Multiple Elimination (SRME) and Radon techniques.Velocity analysis is performed as an iterative process taking into account the geological model.In summary sub-basalt imaging has undergone remarkable advances in last years, these improvements consist in designing new geometry acquisition patterns (White et al., 2003), designing new sources (Staples et al., 1999;White et al., 2002;Ziolkowski et al., 2003), understanding the scattering caused by the basalt (Martini et al., 2001;Martini and Bean, 2002) or combining several geophysical methodologies (Jegen-Kulcsar and Hobbs, 2005).Nevertheless, studying sub-basalt structures requires a detailed velocity model to obtain valuable information and to apply more sophisticated approaches such as prestack depth migration.

Theoretical geologic model and synthetic seismic data
A synthetic dataset was acquired using an idealised velocity model.The model consists of a wedge of basalt layer within a sedimentary column.The basalt layer features an irregular upper surface, and very high velocities and densities that contrast with the velocities and densities of the surrounding sediments.The choice of a thin wedge for the basalt Table 1.Velocities and densities used in the synthetic model.Physical properties were taken from Carmichael (1982).The Poisson ratio was 0.25 and the density was calculated using the Christensen relation Christensen and Mooney (1995) is justified because it represents the best geological setting for exploration and exploitation for the oil industry.In order to simulate a highly variable structure, we consider the small-scale top basalt topography to be a random field which was generated using Von Karman functions (Goff and Jordan, 1988).As there is a high velocity contrast between the sedimentary cover and the basalt, can this topography be recovered?In addition, in the case of Shetland-Faroe Basin, the area where the basalt thins is close to the center of the basin where geology is well known and can be extrapolated to suggest the existence of sub-basalt sedimentary structures.The P-wave velocities were taken from laboratory measurements (Carmichael, 1982).All the physical properties used in the simulation are summarised in Table 1.
A second order finite difference solver of the elastic wave equation using Sochacki's interface scheme (Sochacki et al., 1987(Sochacki et al., , 1991) ) was used to generate a seismic dataset acquired over the velocity model.The 2-D model was of 100 km wide and 8 km deep and a 10 × 10 m grid was used (Fig. 1).After intense calculations, 80 shot-gathers were simulated with more than 3000 traces per shot (one trace every 30 m) and 30 s of recording time.The sampling for this synthetic dataset was 1 ms.The parameters for these simulations are summarised in Table 2.In this data several phases were identified (Fig. 2).The source wavelet is a minimum phase wavelet with a frequency content of 5 to 25 Hz.Travel time data was picked using the peak amplitude as the reference and an estimated uncertainty of 0.010 s was used in the inversion scheme.
Water multiple and peg-leg signal were generated by the elastic finite difference algorithm, no additional noise was included in the data.Although in nature, basalts appear as highly heterogeneous layered structures, in order to simplify the problem, the basaltic wedge was considered as an homogeneous feature in its internal velocity distribution.Under these conditions we obtained a quite ideal dataset.Thus at this points, the thickness of the basaltic wedge was known.This information was used in the inversion (see text for more explanation).

Tomographic inversions
The main aim of synthetic simulations was checking the possibility of recovering the original model using tomographic techniques.As a first attempt, first arrivals travel time seismic tomography was applied using the TTT software package (Trinks et al., 2005).This code is based on initial value ray tracing in velocity models constructed of Delaunay triangulated grids and interfaces.The tomography is implemented as a joint interface and velocity inversion using the biconjugated gradient method.The TTT algorithm presented in Trinks (2003) and Trinks et al. (2005) is based on the method by McCaughey and Singh (1997); Hobro et al. (2003) with the advantage that allows for irregular parametrized-velocity grids using Delaunay triangulation.The velocities are defined at triangle vertices, then a linear interpolation is used to calculate the slowness squared within the triangle.Note that the velocity field is defined by the slowness squared (sloth models (Muir and Dellinger, 1985) this parametrization results in fast calculation of the analytic solution to the ray tracing equations.This approach avoids the discontinuities in ray paths that would appear when using constant velocity cells ( Ćervený, 2001).The appendices in Trinks et al. (2005) go to further detail in mathematics used to fill in the velocity grid.Depending on the ray density the algorithm allows for a denser model parametrization of the densely sampled regions, increasing the resolution and efficiency of the algorithm.The inversion follows Vesnaver (1994); Böhm (1996) scheme that reduces the non-uniqueness of the travel-time inversion result by adapting the grid.
The forward problem is solved by using an analytical ray tracing in a medium with a linear gradient of slowness squared (Farra, 1990;Ćervený, 2001), an initial-value raytracing scheme with traveltime interpolation is used.The model parametrization and the ray-tracing approach used by this algorithm allow for an efficient analytical computation of the Frechet derivatives of travel-time with respect to model parameters.Trinks et al. (2005) demonstrates that the calculation of the Frechet derivatives are only needed at the points along the ray paths that correspond to the travel-times, not at every point of the model grid.Therefore, the inversion approach is very efficient.Finally, Trinks et al. (2005) demonstrates this points using synthetic and real data tests.Further details on the TTT tomographic algorithm can be found in Trinks et al. (2005).
The final model 3 by this adaptive tomographic inversion scheme was reached when it was able to reproduce the the travel time data within the picking uncertainty and χ was approaching 1.The average velocity structure of the zone over the basalt layer was recovered as well as the top of the basalt layer where a sharp velocity contrast is displayed.However, no low velocity can be reproduced under the basalt, toward the right end of the model, where no basalt exists, some realistic information about velocities can be obtained for the deepest part of the model.These results suggest that there are physical limitations in constraining sub-basalt structures by seismic travel time tomography.TTT code can also invert additional phases.In order to include all the information from phases identified in synthetic shots, a layer by layer striping inversion was performed.As the problem is a specific one the starting model was chosen relatively close to the expected solution to determine if the travel-time data alone was able to resolve the sub-basalt structures, taking into account the high velocity contrasts between the sediments and the basalt.In the seismic exploration area, where this study is applicable, the 1 dimensional average structure is relatively well constrained.Therefore, a relatively simple 1-D model was chosen as the starting model.This consists in a five layer cake model these 5 layers correspond to: a water layer, a sediment layer, a basalt layer characterised by high velocity, and second sedimentary sequence and the basement.So velocities typical of this layers were considered.A series of starting models characterised by a maximum variations of 20 % on the velocity values and 15 % in the thickness where tested with no significant variations in the resolved final models.Top basalt interface was inferred from the model obtained using only first arrivals (Fig. 3).Even though in the theoretical model some sub-layers were included in sedimentary sequences, the contrast in velocity between these sub-layers is quite smooth which makes it difficult to identify events from these interfaces, therefore this minor discontinuities were not considered in the inverted model.
The additional information that TTT can use corresponds to the identification of travel-time branches related to the  different features within the model.Therefore, a subjective interpretation of the travel time branch is required, and then the picks corresponding to the branch are associated to a particular structure.

Inverting phases over the basalt layer
Firstly, only the travel time branches interpreted to correspond to the sedimentary cover were included in the inversion.The results show a good recovery of the original model (Fig. 4).In the first part of the model (thin water layer, 0-30 km marked as (A) this phase appears as first break and the results are similar to the first arrivals inversion (Fig. 3) while in the last part of the model (thick water layer, 70-100 km marked as (B) the picks used were not considered in first arrivals inversion, hence, the additional data provides further constraints on the sedimentary cover over the basalt layer.
The TTT code can include also reflected arrivals in the inversion scheme.As reflections from the top of the basalt are displayed as a very high amplitude events, the travel times of the reflected phases can be identified and picked at normal incidence.
Considering the final model of the previous case as starting model, and, without modifying the velocity values for this model, reflections from the top of the basalt layer were inverted in order to obtain the topography for this interface.A detailed structure was achieved which reproduces in some degree the rugged topography featured by the original model (Fig. 4).

Inverting refractions inside the basalt
In this case, the main aim is to constrain the base of the basalt using the refracted waves inside this layer.Raypaths are very sensitive to high velocity anomalies, therefore some constraints on the base of the basalt should be gained by introducing these refractions.As no noise is present in this synthetic dataset, refraction from basalt can be followed up to far offsets (Fig. 2).The maximum offset to stop picking is arbitrary because there is no way to separate basalt refraction from base basalt reflection.In a first picking stage, refractions from the basalt were picked as far as possible and inverted.The results do not fit the theoretical model, overestimating the basalt thickness (Fig. 5).
To simulate a more realistic situation, some arbitrary noise was added to the data (Fig. 6).As a result, the noise limited the maximum offset for picking which was considerably reduced compared with the noise-free case.The inverted model correlates better with the synthetic model (Fig. 7) which is closer to the real model.However, the range between 30 and 75 km features again an overestimation of the basalt thickness.This effect is caused by the difficulty in differentiating basalt refractions from base basalt reflections which interfere at short offsets.In any case, in the travel time branch, the limit between the basalt refraction (head wave) and the subbasalt reflection is completely arbitrary (subjective, depends on the interpreter).The travel time picking of this arrival is complicated even farther by the existence of noise, and thus, the maximum picking offset depends critically on the quality of the data.The consequence of this subjectivity is that the thickness of the basalt layer is proportional to the offset picked, and this effect is stronger where the basalt layer thins.
This contrasts with the conventional idea of using large apertures for sub-basalt imaging.This strategy may be useful in zones with thick basalt layers but, in the case of thin basalt layers, the basalt head wave and reflection from the base of the basalt interfere (see below) and there is no benefit by using large apertures to infer basalt seismic properties.

Inverting refractions and reflections from the basement
As shown above, the basalt layer cannot be resolved properly only considering refractions within this layer.In the case of a thicker basalt layer, reflection from the base of the basalt could be differentiated from the basalt refraction which may contribute to better constrain the base of this layer.However, for thin layers there are no possibilities of deducing the basaltic structure using refraction data.At this point, additional information is required to constrain the basalt layer thickness, in a real case, this additional information could be provided by drilling through the basalt layer, fixing in this way the velocity of the basalt, the thickness of the basalt layer and the velocity of the sediments beneath the basalt.We introduced additional information in the inversion scheme, we assumed that two wells were drilled located at x = 10 km and x = 50 km (Fig. 1).In the last part of the model there is no basalt layer, hence the signal coherence is preserved making it possible to identify and pick normal incidence reflections from the basement.Inverting this phase, a reliable estimation of the top of the basement was obtained for the last 25 km of the model which, jointly with the velocity obtained in first arrival inversions, constrained the model in this part.This results were extrapolated under the basalt layer and used as starting model to invert reflections and refractions from the basement which yielded to our final model where the theoretical model is reasonably well recovered (Fig. 8).Note that additional information is required by the travel time inversion methods to obtain reliable models.

Metropolis simulations
Statistical methods maybe used in order to assess the reliability of inverted velocity models.Among them, probably Monte-Carlo based simulations are the most used.These consist of generating a relatively large number of random velocity models from the same starting model and performing an inversion for each starting model.Finally, the results are compared to asses which model fits the data the best.This method is also used to test the reliability of the starting model in a inversion scheme and the stability of the results (Korenaga et al., 2000;Sallarès et al., 2003;Martí et al., 2006).In this method every iteration is independent from the previous one and no information is inherited for every new case.
Metropolis algorithms (Metropolis et al., 1953) take advantage of the a priori knowledge of the previous model in the iterative scheme.In the first stages of the simulation the influence of the starting model is clear but, after some iterations, this influence decreases considerably, this is known as "burn-in".Using this technique the whole region of allowed parameters is visited which yields a random walk within the region of possible parameters.A scheme of the algorithm is shown in Fig. 9 and a detailed description of the methodology can be found in Pearse.In this case, Rayinvr (Zelt and Smith, 1992) was used to solve the forward problem and to calculate χ 2 which provides the likelihood for every model.The iterative process starts using the starting model (in this case the real model used to build synthetic data) as the current model, then this was randomly modified to generate a new model.Then, the reliability of the model was tested based on the ratio obtained from dividing the likelihood of the current model versus the likelihood of the new model: -If the ratio was equal or larger than one, the new model was accepted and the new model became the current model.
-If ratio was smaller than one, a random number between 1 and 0 was calculated and another test was performed: -If the ratio was larger than the random number, the new model was accepted and the new model became the current model.-If the ratio was smaller than the random number, the new model was rejected and the current model remained unchanged.
The process is repeated for a large number of iterations giving as a result a set of different models that reasonably fit the picked travel times when a picking uncertainty is given.In this study, 40 000 iterations were performed for every case in order to obtain a set large enough to have a statistical value.The allowed variations in velocity, thickness and velocity gradient were: 0 % for water, 2 % for layer over basalt, 5 % for basalt, 10 % for sub-basalt layer and 10 % for basement (Table 3).Variations were chosen increasing in depth to account for the lost of accuracy in deeper layers.Some synthetic shots have been generated using different 1-D models and two different frequencies 10 Hz and 20 Hz (Figs. 10 and 11): Fig. 9. Scheme used in metropolis calculation as described in Pearse (2002).The random modification is subject to prior defined degrees o freedom.In the present case we divide the likelihood of the current versus likelihood of the new model, if ratio is greater than 1 then mode is accepted, if ratio is between 0 and 1 then a random decision is taken based on a prior probability function (linear in this case) that wil preferentially accept models that are close to the accepted boundary.The likelihood ratio is compared with a normalised random number, i ratio is larger than this number then model is accepted.In any other case, the model is rejected.Fig. 9. Scheme used in metropolis calculation as described in Pearse (2002).The random modification is subject to prior defined degrees of freedom.In the present case we divide the likelihood of the current versus likelihood of the new model, if ratio is greater than 1 then model is accepted, if ratio is between 0 and 1 then a random decision is taken based on a prior probability function (linear in this case) that will preferentially accept models that are close to the accepted boundary.The likelihood ratio is compared with a normalised random number, if ratio is larger than this number then model is accepted.In any other case, the model is rejected.
-Model 1: model with sub-basalt low velocity layer.
-Model 2: the same as model 1 but with a thicker basalt layer.

Model 1: results
The first analysis was done on the data from model 1 and 10 Hz.In this case, refraction from basalt layer seems to be very clear as shown in Fig. 10.However, this yields a phase identification which is not correct because the phase identified as a refraction from basalt layer is the interference between two phases: refraction from basalt and reflection from the base of the basalt.To emphasize this effect, a shotgather was calculated using model 1 but replacing the sub-basalt velocity by a layer of 5 km s −1 , in this way no reflection in the base of the basalt layer was generated obtaining a pure refraction in the basalt layer.This pure refraction was picked and  compared with the identified in the previous case (Fig. 12).Thus, in the first case, we have identified an interference between basalt refraction and base basalt reflection as a pure refraction which is erroneous.
Using the picks from the data obtained with model 1 in the Metropolis approach with a high picking error (see Table 4) and considering 40 000 different cases, we obtained  an overestimation in both, velocity and thickness of the subbasalt layer (Fig. 13).By reducing the picking error (case with low uncertainty) we would expect a better correlation between the more probable model and the theoretical one.In practice, under the same conditions but reducing the picking uncertainty, a worse result was obtained where there was no need of a sub-basalt low velocity layer in order to fit the data (Fig. 14).This effect can be explained because considering a bigger error in the picks, the range of times can include both, refraction within the basalt and reflection from the base of the basalt.Therefore, what was labelled as basalt refraction is within the allowed range of times for this phase.On the other hand, by reducing the picking error, the range of times do not include the real refraction and then our erroneous phase identification yields the unexpected result of Fig. 14.
The same analysis was repeated using model 1 and data generated with 20 Hz.In this case, results are better and fit the right model (Fig. 15).

Model 2: results
In model 2 the thickness of basalt layer was increased to test if it was possible to use a well identified base basalt reflection to constrain better the thickness and velocity of this layer.Two different simulations have been performed: one using a very conservative picking and avoiding picks in the "interference zone" (basalt refraction/base basalt reflection) displayed in Fig. 12 and another including more picks.For the first case (Fig. 16) the velocity gradient for basalt layer was well reproduced while for the second case (Fig. 17) this gradient did not fit the real model.For the basement velocity gradient the result was the opposite, obtaining a better result when considering a larger number of picks.In both  cases, sub-basalt velocity layer is not reliably recovered.As in model 1, better fit is obtained for 20 Hz data (Fig. 18), where the velocity gradient for basalt and basement are well reproduced.Again, in both cases, the sub-basalt layer is not well recovered.The Metropolis study reveals that the phase identification is a critical step.Despite objectivity provided by mathematics used in the inversion, phase identification turns travel time tomography in a subjective procedure.Additionally, uncertainty is also a critical parameter in the inversion which can influence the inversion algorithm.Moreover, considering data with different frequency content also has an effect on the selection of the most probable model.Not all models required a low velocity layer and there was an unresolvable trade-off between thickness and velocity, even for models where the base basalt reflection could be identified.Synthetic shots were created using a full-waveform code while likelihood was calculated using a ray tracing code.Full-waveform techniques are more accurate than ray tracing methods because they take into account information about the amplitudes, which are ignored in ray tracing simulations.Due to the high computational cost of full-waveform methodologies, ray tracing methods are still conventionally used to obtain velocity models (Pratt et al., 1996).These results suggest that conventional travel time inversion (tomography) schemes without additional information are not sufficient to constrain the base of basalt or sub-basalt geological structures.
In the synthetic tests performed in this study no noise has been considered.In this sense, some features that commonly are present in real data as tidal noise, electrical noise and some other effects affecting the quality of the data, have not been included.Results derived from these tests must be    interpreted as results obtained in ideal conditions and in consequence, the best ones expected for a real case using this methodology.

Conclusions
This study, which involves synthetic data suggests that there are some physical limitations to obtain a reliable velocity model for sub-basalt zones in areas covered by high velocity rocks (like basalts and salts).In the case of thin basalt layers, the base basalt reflection is totally masked within the waterwave cone and it cannot be separated from the basalt refraction.There are several subjective factors that can affect and condition the results from the inversion as maximum picking offset or picking uncertainty.Another important point is the frequency content of the signal, our Metropolis simulations suggest that the original model is best recovered using high frequencies and thicker basalt.This result is relevant because the actual tendency is using and designing airguns that produce low frequency data as single bubble source.The most critical point in the travel time inversion is the phase identification/interpretation in the shot record.Differentiating in the travel time branch, between the head wave travelling within the basalt (refraction) and, the base basalt reflection is a key element in determining the correct thickness of the basalt.The uncertainty associated to the travel time picks is also a relevant issue, as it can not distinguish between high and low velocity sub-basalt structures.Moreover, a wrong determination of the basalt thickness and velocity has a direct influence on the resulting model for layers under the basalt.Reliable sub-basalt imaging with wide-angle reflection/refraction datasets requires additional information as the knowledge on the thickness of the basalt at some point and its internal velocity distribution.This could be achieved by using other methodologies to infer basalt properties.

Fig. 1 .
Fig. 1.Synthetic velocity model resampled using 100 × 100 m squared cells.The original model used to run simulations was sampled by 10 × 10 squared cells.Vertical lines at 10 and 50 km show the location of hypothetical wells drilled through the basalt layer.Thus at this points, the thickness of the basaltic wedge was known.This information was used in the inversion (see text for more explanation).

Fig. 2 .
Fig. 2. Synthetic shotgather.Different phases can be identified: water wave (blue), refraction from sediments over the basalt layer (green), refraction from basalt (red), reflection from the top of the basalt (yellow), refractions from basement (purple) and reflections from the top of the basement (orange).

Fig. 3 .
Fig.3.Velocity model obtained using TTT packageTrinks et al. (2005) considering only first arrivals.Interfaces between layers are the ones that were used in the theoretical model.Velocities over the basalt wedge are well constrained.A prominent velocity discontinuity can be observed for the first 70 km between 3 and 4 km in depth which provides the location for the top of basalt layer.A 1 dimensional five layer-cake model was used as the starting model (water, sediments, a wedge of basalt, sediments and basement).Note that the interface between the sediment and the underlying basalt was determined using only the first arrivals.No significant differences in the final model were observed by changing the velocity model 20% in the velocities and 15% in the layer finickinesses.The unresolved parts of the model at both ends of the transect are determined from forward modelling, ray tracing.

Fig. 4 .
Fig. 4. Velocity model obtained using TTT package Trinks et al. (2005) after inverting refractions from the sediments over the basalt and reflections from the top of the basalt layer.Dashed lines represent layers from inversion and continous lines the theoretical layer interfaces Note the coincidence between dashed line and continous line in the top of the basalt layer.The unresolved parts of the model at both ends of the transect are determined from forward modelling, ray tracing.

Fig. 3 .
Fig. 3. Velocity model obtained using TTT package Trinks et al. (2005) considering only first arrivals.Interfaces between layers are the ones that were used in the theoretical model.Velocities over the basalt wedge are well constrained.A prominent velocity discontinuity can be observed for the first 70 km between 3 and 4 km in depth which provides the location for the top of basalt layer.A 1 dimensional five layer-cake model was used as the starting model (water, sediments, a wedge of basalt, sediments and basement).Note that the interface between the sediment and the underlying basalt was determined using only the first arrivals.No significant differences in the final model were observed by changing the velocity model 20 % in the velocities and 15 % in the layer finickinesses.The unresolved parts of the model at both ends of the transect are determined from forward modelling, ray tracing.

g. 4 .
Velocity model obtained using TTT package Trinks et al. (2005) after inverting refractions from the sediments over the basalt and ections from the top of the basalt layer.Dashed lines represent layers from inversion and continous lines the theoretical layer interfaces.te the coincidence between dashed line and continous line in the top of the basalt layer.The unresolved parts of the model at both ends of e transect are determined from forward modelling, ray tracing.

Fig. 4 .
Fig. 4. Velocity model obtained using TTT package Trinks et al. (2005) after inverting refractions from the sediments over the basalt and reflections from the top of the basalt layer.Dashed lines represent layers from inversion and continuous lines the theoretical layer interfaces.Note the coincidence between dashed line and continuous line in the top of the basalt layer.The unresolved parts of the model at both ends of the transect are determined from forward modelling, ray tracing.

Fig. 5 .
Fig. 5. Results from the basalt refraction inversion.White dashed lines represent layers from inversion and continuous lines the theoretical layer interfaces.The base of the basalt layer, which is overestimated, should be delineated by raypaths.The black area represents the part of the model sampled by rays.

Fig. 6 .
Fig. 6.Synthetic shotgather with noise.Different phases can be identified: water wave (blue), refraction from sediments over the basalt layer (green), refraction from basalt (red), reflection from the top of the basalt (yellow), refractions from basement (purple) and reflections from the top of the basement (orange).Note the difference with Fig. 2 in the refractions from basalt (red).

Fig. 7 .
Fig. 7. Results from the basalt refraction inversion using picks from data with noise.White dashed lines represent layers from inversion and continuous lines the theoretical layer interfaces.The black area represents the part of the model sampled by rays.The base of the basalt layer should be delineated by raypaths.The constrain on the thickness of the basalt layer is failing where the layer is thinner, probably due to overpicking refractions.

Fig. 8 .
Fig.8.Final result obtained using all the phases after fixing the base of the basalt layer considering that two wells were drilled through this layer at 10 and 50 km.Dashed lines represent layers from inversion and continuous lines the theoretical layer interfaces.Introducing additional information, the theoretical model is recovered quite accurately.Energy dispersion caused by the topography of the basalt wedge masks the reflected energy from the basement.Therefore, the recovered basement has an irregular top which is not real but the influence of the overlying velocity heterogeneities.

Fig. 10 .
Fig. 10.1-D model used to generate synthetic data (top).Shots generated using the model and different frequencies: 10 Hz (left) and 20 Hz (right).Main phases were identified: sea bottom reflection (blue), top basalt reflection (green), basalt refraction (red), top basement reflection (yellow) and basement refraction (orange).Under the top of the basalt no phases were picked within the water-wave cone because in real data this phases are difficult to identify.

Fig. 11 .
Fig. 11.1-D model used to generate synthetic data (top).Shots generated using the model and different frequencies: 10 Hz (left) and 20 Hz (right).Main phases were identified: sea bottom reflection (blue), top basalt reflection (green), basalt refraction (red), base basalt reflection (purple), top basement reflection (yellow) and basement refraction (orange).Under the top of the basalt no phases were picked within the water-wave cone because in real data this phases are difficult to identify.

Fig. 12 .
Fig.12.Basalt refraction picks for model with sub-basalt low velocity layer (red) and for a model without sub-basalt low velocity la (cyan).In the case without sub-basalt low velocity layer the picks represent a pure refraction while in the other case, the phase tha identified as a refraction is made by the basalt refraction interfering with the base basalt reflection.

Fig. 12 .
Fig. 12.Basalt refraction picks for model with sub-basalt low velocity layer (red) and for a model without sub-basalt low velocity layer (cyan).In the case without sub-basalt low velocity layer the picks represent a pure refraction while in the other case, the phase that is identified as a refraction is made by the basalt refraction interfering with the base basalt reflection.

Fig. 13 .
Fig. 13.Results obtained after using the Metropolis algorithm on data from model 1 and 10 Hz for 40 000 cases.Red line represents the real model and every black line a modified model.The colour scale stands for the number of times that a model (or part of it) is visited.The preferred model (blue colours) overestimates the subbasalt layer thickness as well as the velocity for this layer.

Fig. 14 .
Fig. 14. Results obtained after using the Metropolis algorithm on data from model 1 and 10 Hz for 40 000 cases.Uncertainties were reduced in comparison with the previous case (Fig. 13).Red line represents the real model and every black line a modified model.The colour scale stands for the number of times that a model (or part of it) is visited.The preferred model consists in a velocity gradient which includes basalt, and basement layers, avoiding the need of a low velocity layer under the basalt.

Fig. 15 .
Fig. 15.Results obtained after using the Metropolis algorithm on data from model 1 and 20 Hz for 40 000 cases.Red line represents the real model and every black line a modified model.The colour scale stands for the number of times that a model (or part of it) is visited.The preferred model consists in a velocity gradient which includes basalt, sub-basalt and basement layers.The sub-basalt low velocity layer is reasonably well recovered in velocity and thickness.

Fig. 16 .
Fig. 16.Results obtained after using the Metropolis algorithm on data from model 2 and 10 Hz for 40 000 cases considering a "conservative" picking avoiding picks in the "interference zone".Red line represents the real model and every black line a modified model.The colour scale stands for the number of times that a model (or part of it) is visited.

Fig. 17 .
Fig. 17. Results obtained after using the Metropolis algorithm on data from model 2 and 10 Hz for 40 000 cases considering more picks than in the previous case (Fig. 16).Red line represents the real model and every black line a modified model.The colour scale stands for the number of times that a model (or part of it) is visited.

Fig. 18 .
Fig. 18. Results obtained after using the Metropolis algorithm on data from model 2 and 20 Hz for 40 000 cases.Red line represents the real model and every black line a modified model.The colour scale stands for the number of times that a model (or part of it) is visited.The preferred model consists in a velocity gradient which includes basalt, sub-basalt and basement layers, avoiding the need of a low velocity layer under the basalt.

Table 2 .
Parameters used to generate synthetic data.

Table 3 .
Allowed variation to generate modified models for velocity, velocity gradient and thickness for every layer in models 1 and 2.

Table 4 .
Picking uncertainty in ms for every layer considered in metropolis algorithm.