adg vol5 n02 Aver247_258.pdf ANNALS OF GEOPHYSICS, VOL. 45, N. 2, April 2002 247 Depth model building by constrained magnetotelluric inversion Paolo Dell’Aversana and Sergio Morandi Enterprise Oil Italiana S.p.A., Roma, Italy 1. 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. Mailing address: Dr. Paolo Dell’Aversana, Enterprise Oil Italiana S.p.A., Via dei Due Macelli 66, 00187 Roma, Italy; e-mail: Paolo.DellAversana@rome.entoil.com Abstract In this paper we describe an approach aimed at integrating seismic and magnetotelluric data in a complex geological setting, characterised by thrust structures, in Southern Apennine, Southern Italy. Seismic data were collected by the «Global Offset» technique that is designed to record high fold data in a wide range of offsets, without losing the benefit of near vertical reflection seismic. First arrivals picked from short to long offsets and the main reflections were inverted in order to produce a tomographic velocity-interface model. It was converted into a resistivity section applying an empirical relationship, obtained by well logs, between resistivity and velocity. That section was used as a reliable reference model for 2D inversion of magnetotelluric data collected along a parallel section very close to the seismic profile. The process was iterative and interactive and was aimed at producing consistent velocity and resistivity sections, honouring seismic and MT data set. The final MT model fits very well the observed apparent resistivity and phase, reproduces the main geological trends and is constrained by a well drilled close to the line. Key words magnetotelluric seismic inversion integration 248 Paolo Dell’Aversana and Sergio Morandi Based on that point, during the last five years Enterprise Oil developed an acquisition tech- nique, 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, 2000). That goal can be obtained by using multi-channel re- cording 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 stand- alone 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 mag- netotelluric, in terms of both recording equip- ment, 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, magneto- telluric 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 (Dell’Aversana et al., 2001). 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, ob- tained by sonic and resistivity logs. The inte- gration 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 in- version 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 t h e G l o b a l O ff s e t t o m o g r a p h i c m o d e l . 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). 2. Seismic data Figure 1 shows the final stack section produced applying a conventional processing 249 Depth model building by constrained magnetotelluric inversion flow on a seismic line about 15 km long recorded, in S. Apennines (S. Italy), by using a standard roll along the layout. Low re- flectivity 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., 1999, 2000; Improta et al., 1998, 2000a,b, 2001). Fig. 1. Conventional time-stack section obtained by a conventional processing flow in CDP domain. Many factors can negatively affect the quality of seismic imaging in thrust belts, such as low signal to noise ratio, scattering phenomena, bad static corrections and so on. The final effect can be a seismic section very difficult to interpret, like that shown in the figure. 250 Paolo Dell’Aversana and Sergio Morandi 3. 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 ma- gnetotelluric 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 con- ductive 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 plat- form. 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. 4. MT rotational analysis Figure 2 shows resistivity and phase, for both TE and TM modes on the left panel for each MT site. On the right panel, 1D models are shown: the blue curve represents the Bostic model; the purple curve is the Occam model; the green curve is the corresponding layered model (Guess layered model, fixing the number of possible layers). 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 suf- ficient 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. In order to improve the model building pro- cess, 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 251 Depth model building by constrained magnetotelluric inversion Fig. 2. Three examples of MT curves and corresponding TE 1D models (AA04R, AA22R, AB26R, starting from top). See details in the text. 1 10 10010 252 Paolo Dell’Aversana and Sergio Morandi Fig. 3. Rotational analysis for the MT station BG06BR. Several rotations were applied to the axis in order to evaluate possible 3D effects. For that purpose impedence skew and ellipticity were analysed for each rotation. A 253 Depth model building by constrained magnetotelluric inversion Fig. 4. Polar diagram for MT site AA02R, at the beginning of the profile. 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. 5. 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, underdetermi- ned; as a consequence non-uniqueness prob- lems 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 con- straints and a reliable starting model could be a good strategy to reduce the space of accep- table 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 tomo- graphy. All the first breaks and some reflection events were inverted to produce a parametric Zxy impedence Zxx impedence Impedence strike Tipper mag. Strike Ind. Arrow o 1 0 1 * 254 Paolo Dell’Aversana and Sergio Morandi depth section to be used as a guide for building a reliable reference model for 2D MT inversions. In other words we performed a transmis- sion/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 re- ference 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 (m) = (d – F(m))T V -1 (d – F(m)) + mT LT L m for given , V and L. F is a forward modelling function, d is the data vector. is called «regu- larisation 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 Fig. 5. Final magnetotelluric model. It honours MT data and, at the same time, it derives from a starting model based on Global Offset data. Note the presence of a high resistivity thrust system, low resistivity shallow basins, low resistivity layers before reaching the high resistivity platform at about 6 km b.s.l. Z (m) 1000 0 -1000 -2000 -3000 -4000 -5000 -6000 X (km) 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 4096 2048 1024 512 256 128 64 32 16 8 4 2 Rho (Ω•m) 255 Depth model building by constrained magnetotelluric inversion shallowest part of the section. Also a quite con- tinuous 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. Fig. 6a. Pseudo-sections for TE mode. TE and TM calculated pseudosections were compared with the observed pseudo-sections in terms of resistivity (upper panels) and phase (lower panels). The fit ranges from good to very good, indicating that the model of fig. 5 honours the MT data. App Rho Obs App Rho Calc Phase Obs Phase Calc 10-2 10-1 100 101 102 103 262144 131072 65536 32768 16384 8192 4096 2048 1024 512 256 128 64 32 16 8 4 2 App. Rho ( •m) 90 80 70 60 50 40 30 20 10 Phase (deg.) 10-2 10-1 100 101 102 103 10-2 10-1 100 101 102 103 10-2 10-1 100 101 102 103 256 Paolo Dell’Aversana and Sergio Morandi Fig. 6b. Pseudo-sections for TM mode. TE and TM calculated pseudosections were compared with the observed pseudo-sections in terms of resistivity (upper panels) and phase (lower panels). The fit ranges from good to very good, indicating that the model of fig. 5 honours the MT data. 6. 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. 257 Depth model building by constrained magnetotelluric inversion 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 infor- mation 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, con- sistently 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 con- sistent with the borehole information available along the line. REFERENCES DELL’AVERSANA, P. (2001): Integration of seismic, magne- totelluric and gravity data in a thrust belt interpretation, First Break, 19 (6), 335-341. DELL’AVERSANA, P. and S. MORANDI (2000): Integrated inversion of seismic and magnetotelluric data sets in thrust belt, in 62nd EAGE Conference, Glasgow (extended abstract). 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 Break, 18, 87-96. DELL’AVERSANA, P., D. COLOMBO, S. MORANDI and M. BUIA (2001): Thrust belt exploration by «Global Offset» seismic and reflection/refraction tomography, in 63rd EAGE Conference, Amsterdam (extended abstracts). IMPROTA, L., P. DELL’AVERSANA, A. HERRERO, S. MORANDI and A. ZOLLO (1998): Tomografia sismica di strutture crostali superficiali in Appennino meridionale mediante inversione di dati acquisiti in configurazione «Global Offset», in 17th GNGTS Meeting, Rome, Italy (abstract). IMPROTA, L., M. FRATTINI, A. ZOLLO, J. VIRIEUX, A. HERRERO and P. DELL’AVERSANA (2000a): Non-linear inversion of global offset reflected arrivals in highly heterogeneous shallow crustal structure, in European Geophysical Society Meeting, Nice, France (abstract). IMPROTA, L., A. ZOLLO, J. VIRIEUX, A. HERRERO and P. DELL’AVERSANA (2000b): Mapping interfaces in an overthrust region by non linear travel time inversion of reflection data, Am. Geophys. Un., Fall Meeting (abstract). IMPROTA, L., A. HERRERO, A. ZOLLO, P. DELL’AVERSANA and S. MORANDI (2001): 2D non-linear travel time inversion of first arrivals and reflection data: application to the Southern Apennines thrust-belt (Italy), in SEG, San Antonio (abstract). MENKE, W. (1989): Geophysical data analysis: discrete inverse theory, Int. Geophys. Ser., vol. 45. RODI, W. and R.L. MACKIE (2001): Non-linear conjugate gradients algorithm for 2D magnetotelluric inversion, Geophysics, 66 (1), 174-187.