Depth model building by constrained magnetotelluric inversion

. DELL’AVERSANA, P., S. MORANDI, L. IMPROTA and A. ZOLLO (1999): Tomographic images from «Global Offset» seismic data inversion in Southern Apennines, in 61st EAGE Conference, Helsinki (extended abstracts). DELL’AVERSANA, P., E. CERAGIOLI, S. MORANDI and A. ZOLLO (2000): A simultaneous acquisition test of high density «Global Offset» seismic in complex geological settings, First


Introduction
An apparent paradox can arise when many different data have to be used to build a model of a complex geological system: the more data are available the less reliable and stable the model. Apparently it seems that collecting more «objective» information has the effect of increasing the «subjective» character of the solution for a given problem. This situation is common for many categories of complex problems, where there is a predominance of non-linear relationships and non-linear processes between data and model spaces and when the data set is itself something like an advanced interpretation. In exploration in a thrust belt environment, that situation is frequent, especially when the quality of the near vertical reflection seismic data is poor or very poor: in that case many efforts are addressed to apply alternative geophysical tools, such as wide angle reflection seismic, magnetotelluric, gravity and so on. Also, during the last ten years, special efforts have been oriented to improve the efficiency of the seismic acquisition layouts in order to improve the quality of the data in difficult areas. In particular, the technological improvements of the recording systems (especially, but not only, in terms of number of available recording channels) yielded seismic data characterised by a fold much higher than in the past. For the same reasons, also the possibility to collect large offset seismic data has been increased significantly, without losing the benefits offered by multi-fold data sets.
Based on that point, during the last five years Enterprise Oil developed an acquisition technique, called «Global Offset», that is aimed to optimise simultaneously the benefits of both near-vertical reflection seismic, wide angle seismic and transmission-reflection tomography (Dell'Aversana et al., 1999. That goal can be obtained by using multi-channel recording systems typical of standard reflection seismic and conventional dynamite shots. The particularity is that the geophone layout is kept active throughout the seismic profile, so that each shot can be recorded, simultaneously in the whole range of available offsets, by densely spaced stations. The final data sets recorded by that kind of layout contain a huge amount of information (in terms of large offset turning rays and wide angle reflections) that is commonly lost when standard seismic, short offset layouts are used. At same time, the receiver coverage is maintained very high. A similar data set can be considered optimal for joint inversion of first breaks and reflections or, in other words, for transmission/reflection tomography. Based on that idea, in the last few years, Enterprise Oil has been acquiring Global Offset data in S. Apennines (S. Italy), also using additional standalone stations in case of long seismic profiles and/or 3D surveys.
On the other hand, strong improvements have recently been obtained for non seismic methods too, such as in the case of magnetotelluric, in terms of both recording equipment, acquisition scheme, processing tools and inversion algorithms. It is clear that the effects of these improving technologies reflect also on the exploration approach. Especially in difficult areas where standard seismic methods fail to produce interpretable data volume, a possible solution can be to try to integrate the benefits derived by complementary methodologies and data sets, such as near vertical reflection seismic, wide-angle seismic, tomography, magnetotelluric and so on. In that viewpoint, in the period 1998-2000, Eni Agip and Enterprise Oil, together with the European Union, founded a research project aimed to optimise the integration of Global Offset seismic data with magnetotelluric and gravity data in a difficult thrust belt area of S. Italy .
In this paper we discuss a similar approach based on the integration of a tomographic seismic 2D model, obtained by Global Offset data, and a magnetotelluric section obtained by MT data collected along a profile parallel to the seismic line. The two parametric models (velocity and resistivity sections) have been linked together using empirical relationships between the two physical parameters, obtained by sonic and resistivity logs. The integration approach was based on an iterative and interactive procedure: before running any 2D MT inversion, a reliable reference model was obtained by inverting the whole data set recorded by a Global Offset acquisition layout. After obtaining a reliable tomographic velocity model, a resistivity section was derived from it using empirical relationships. Also 1D MT inversions and borehole data were used to constrain the resistivity initial model. The successive step was to use that section as a reference model for 2D MT inversion based on a perturbative approach. Using that procedure, the final resistivity model can be considered the results of an integrated process based on borehole, MT and Global Offset data. In that sense, we can say that the 2D MT inverse problem is solved constraining the solution in the direction of the Global Offset tomographic model. Obviously that approach does not guarantee a unique solution for such a complex inverse problem, but at least it produces a consistent model honouring independent data set. It can be considered itself an interesting objective to realise in a difficult exploration area.
In the present paper, we discuss mainly the magnetotelluric aspects of our integration approach. Additional details on Global Offset data, transmission/reflection tomographic section, resolution synthetic tests and gravity data integration can be found in other specific papers (Dell'Aversana and Morandi, 2000; Dell'Aversana, 2001). Figure 1 shows the final stack section produced applying a conventional processing flow on a seismic line about 15 km long recorded, in S. Apennines (S. Italy), by using a standard roll along the layout. Low reflectivity and poor quality of the stacked data are quite evident, although the theoretical fold was 60%. It is clear that an interpretation based only on that section could lead to a wrong model in time and in depth domains. The complexity of the geophysical/geological problem requires a more appropriate approach integrating the standard seismic data with other independent information. In our specific case, the acquisition layout consisted of a sym-metrical spread (maximum offset = 3.6 km, move-up = 60 m, intertrace = 30 m, shot distance = 60 m). That layout was integrated by a series of stand-alone stations in order to define a «Global Offset» acquisition layout. Vertical 4.5 Hz geophones every 90 m and horizontal inline geophones every 450 m were kept active along the whole seismic line, in order to record the seismic energy produced by each single dynamite shot, in the whole range of available offsets (Dell'Aversana et al., 1999Improta et al., 1998Improta et al., , 2000aImprota et al., ,b, 2001.

Magnetotelluric data
Assuming the presence of a conductive flysch above the resistive carbonate target (on the basis of a regional understanding of the geological setting), it was expected that magnetotelluric would have been useful to highlight the top of the platform. To obtain that goal an MT survey was performed in the same area. In particular, about 20 MT sites were acquired in a profile parallel to the seismic line shown in fig. 1. The first interpretation step of the MT data was 1D inversion for each site along the profile. Assuming that TE mode is less affected or distorted than TM by 3D lateral heterogeneity, 1D inversions were performed on Transverse Electric mode. In a second step, inversion tests of TM and 1D invariant mode were also performed.
The one-dimensional models show the main features of the explored medium ( fig. 2), with resistive thrusts alternating with more conductive layers. In particular, a thick (2-3 km) conductive body seems to be present above the resistive basement in almost all MT sites. The basement shows a lower resistivity (both apparent and modelled by 1D inversion) than expected; in fact the resistive layer at about 6.5-7.5 km depth could correspond to the limestone platform, which usually presents a much higher resistivity. A possible explanation is that the thick conductive layer overlaying the carbonates induces an effect of under-estimation in the modelled resistivity associated with the platform. Also, the apparent resistivities are pushed down at middle-low frequencies. It is clear that the 1D inversion process alone cannot solve the problem of a correct definition of a so complex structure, although it discloses some of its important general features. Lateral and vertical complexity is quite evident from the examples that are shown: MT curves vary significantly with location and frequency. The presence of resistive thrusts and conductive layers is clear, down to the resistive basement, although 1D modelling is not sufficient for a quantitative estimation of true resistivity and thickness: 2D and 3D effects can still affect the model parameter definition. In any case the presence of a thick conductive layer covering the resistive basement is evident, especially in the first two examples. It produces an under-estimation effect on the apparent resistivity of the basement itself.

MT rotational analysis
In order to improve the model building process, an accurate 2D magnetotelluric inversion flow was performed.
In that approach, a 2D structure was assumed, although a 3D hypothesis would be more appropriate for the investigated geological setting. However, our hypothesis seems to be partially confirmed by rotational analysis. In fact we analysed several impedence parameters for varying rotational angles. Figure 3 shows impedence skew and ellipticity, for site BG06BR (at km 12 on the profile). The Tipper, resistivity and phase curves are also shown. The skew can be considered a measure of three-dimensionality. It is an invariant for rotation of co-ordinates and should be zero for noise-free data, in 1D and 2D cases. On the other hand, impedence ellipticity varies with setup direction: it is zero only for 1D data and in 2D case when x or y axis is along the strike direction. In the example of fig. 3, the skew is about constant, for different rotational angles and frequency values, showing low values of about 0.1, due to some noise in the data. Ellipticity varies with angles: the lower values appear for rotation of 15°, indicating the direction of the strike (or the orthogonal direction) in the area. Tipper values vary between 0.1 and 0.4, remaining within an acceptable range for the main range of frequency.
In any case, the 2D hypothesis is not perfectly valid for all MT sites along the profile and for the whole range of frequency. For example, in site AA02R, the skew values range between 0.1 and 0.5 and increase at low frequency, indicating    3D effects. They are also confirmed by looking at polar diagrams ( fig. 4), where it is clear that, especially at low frequencies, Z xx impedence is far from zero, as it should be in a perfect 2D case. It means that our 2D-inversion approach can solve only partially the problem of reconstructing a very complex three-dimensional structure, and that 3D inversion would be more appropriate to that exploration case. At the moment we are testing the possibility to apply 3D inversion algorithms to our data, but it will be the subject of a future study.

2D magnetotelluric inversion
Returning to our profile we see that MT data distribution is quite irregular. In some cases, the MT stations projected on the line appear very close (less than 1 km) whereas, in other parts of the profile, the spacing is 6 km or more (AA16R and AA22R). That is a consequence of the variations in data quality (strongly dependent on local noise effects) that induced a necessary selection of the data to be inverted: in fact some bad data were rejected producing a gap of information. In that condition, it is highly probable that the MT 2D inverse problem is, in the best case, mixed determined (Menke, 1989) and, in the worst case, underdetermined; as a consequence non-uniqueness problems have to be expected. This is a quite common problem in MT surveys, apart from those cases where a redundant acquisition layout is used. It is clear that using good constraints and a reliable starting model could be a good strategy to reduce the space of acceptable geophysical models.
From that standpoint, we performed a series of 2D MT inversions using a robust reference model built on the basis of all available Global Offset seismic data. In particular, the large amount of seismic data collected within a wide range of offsets (using stand alone stations too), proved very useful to image shallow and deep structures by transmission/reflection tomography. All the first breaks and some reflection events were inverted to produce a parametric In other words we performed a transmission/reflection tomography using all first breaks available, near vertical and wide-angle reflections. After obtaining that tomographic model we transformed it into a resistivity model using an empirical relationship between velocity and resistivity. For that purpose, we used sonic and resistivity log information available from a well drilled along the profile. The analytic form of that relationship is v = a [ln (ln (rho))] + b where a and b are two constants, v is the velocity by sonic logs, rho is the resistivity by logs.
The above resistivity section was completed by using the indications derived by 1D MT inversions. Finally it was used as a reliable reference model for 2D magnetotelluric inversion.
The used MT inversion method was a non-linear Conjugate Gradient Algorithm (Rodi and Mackie, 1998). It is based on searching a «regularised solution» minimising an objective function, , defined by for given , V and L. F is a forward modelling function, d is the data vector. is called «regularisation parameter» and is a positive number. V is a variance-covariance matrix of the error vector e. The second term of (m) is a «stabilising functional», where Lm approximates the Laplacian of log . A satisfactory fit has been obtained for almost all MT stations (figs. 6a and 6b) for both TE and TM models. Figure 5 shows the final MT model. Colours represent the «true» resistivity distribution. High resistivity thrusts (in yellow) are very clear. Some low resistivity (in blue) basins appear in the shallowest part of the section. Also a quite continuous low resistivity layer is present below the thrusts; at a depth of about 6 km b.s.l., the high resistivity layer of the carbonate platform appears.
On the basis of synthetic tests, the accuracy of the final resistivity model was estimated. The final model was used to simulate the real MT experiment using the true acquisition geometry: the synthetic curves were inverted using the same reference model based on seismic data. The final model was reproduced with an accuracy of about 5%. Although this result does not represent a true estimate of the resolution, it confirms the stability of the 2D inversion when a reliable starting model is used. That is also true in the case of sparse data, as in our real experiment. Finally, another local confirmation of the good model accuracy was obtained by comparing the final MT section with the available borehole information.

Conclusions
In complex geological settings, such as thrust belt environment, the standard seismic approach often produces low quality images in both time and in depth domains. In those cases, a MT survey can significantly improve the process of depth model building. On the other hand, many additional problems can be associated with MT data utilisation, such as non-uniqueness and low resolution of the final models. 1D inversions present all the limits associated with a strong simplification of the real 3D nature of the problem. On the other hand, especially when the distance between stations is large and irregular, the reliability of 2D inversion can be very low and the final solution of the MT inverse problem can be extremely unstable. If the inversion is based on a perturbative approach, the reference model plays a fundamental role in the process.
In this sense, the choice of a robust starting model is very important in order to obtain reliable final 2D models. In this paper, we discussed an approach aimed to produce consistent and reliable magnetotelluric models by 2D MT inversions based on reliable reference models. The basic idea to build realistic reference sections was to use the whole available seismic information contained in a range of offsets much wider than the standard range of the conventional layout. In fact, increasing the maximum offset is a way to produce supercritical reflections and turning rays exploring at considerable depth. The tomographic inversion of the first breaks and of the wide-angle reflections produced a reliable parametric model in depth that guided the process of building a reliable starting model for 2D MT inversion. This approach was applied to a real case of exploration in a thrust belt environment (S. Italy). An integrated model honouring the consistency between Global Offset, MT, borehole and geological data sets was obtained. It reflects the real complexity of the explored geological system, although additional 3D analysis could be appropriate in our case. Strong lateral and vertical geophysical variations appear, consistently with the presence of high resistive thrusts alternating with low resistive layers. The top of the resistive carbonate platform has been clearly highlighted. The resolution and the level of accuracy have just been estimated: more accurate and appropriate tests are required to establish the level of non-uniqueness of the final solution. The process of depth model building described in this paper is clearly iterative and based on a non-linear approach. It means that a full comprehension of the real accuracy and resolution is not easy. On the other hand, it was possible to test the stability of the solution of the MT inverse problem by synthetic tests. They confirmed that the model obtained using a reliable starting model is very stable and consistent with the borehole information available along the line.