Three-dimensional interpretation of MT data in volcanic environments ( computer simulation )

The research is aimed, first, to find components of MT-fields and their transforms, which facilitate the imaging of the internal structure of volcanoes and, second, to study the detectability of conductivity variations in a magma chamber due to alterations of other physical parameters. The resolving power of MT data with respect to the electric structure of volcanic zones is studied using software developed by the author for three-dimensional (3D) numerical modeling, analysis and imaging. A set of 3D volcano models are constructed and synthetic MT data on the relief Earth's surface are analysed. It is found that impedance phases as well as in-phase and quadrature parts of the electric field type transforms enable the best imaging of the volcanic interior. The impedance determinant is, however, the most suitable for adequate interpretation of measurements carried out for the purpose of monitoring conductivity variations in a magma chamber. The way of removing the geological noise from the MT data by means of its upward analytical continuation to the artificial reference plane is discussed. Interpretation methodologies are suggested aimed at 3D imaging and monitoring volcanic interiors by MT data.


Introduction
Electromagnetic (EM) (in particular, magnetotelluric (MT)) fields are widely used to study the geodynamic processes in geothermal and volcanic zones (Park and Torres-Verdin, 1988;Mogi and Nakama, 1990;Mauriello et al., 1997, etc.) due to their deep penetration into the Earth and ability to resolve the parameters of complex geological media.However, in most of the studies so far, 1D or 2D interpretation tools have been used.For instance, based on MT data it has ____________ Mailing address: Dr. Vjacheslav Spichak, Geoelectromagnetic Research Institute, Russian Academy of Sciences, P.O.Box 30, 142190 Troitsk, Moscow Region, Russia;e-mail: v.spichak@g23.relcom.rurecently been possible to reconstruct a 2D resistivity cross-section beneath a profile crossing Mt.Somma-Vesuvio (Di Maio et al., 1998).Meanwhile, long-term forecast of the earth activity should be evidently based on the knowledge of the deep 3D volcanic structure as well as on our ability to interpret the measured data properly.
At a first glance, the construction of a 3D geoelectric model as well as the monitoring of crucial parameters would require synchronous MT measurements carried out at sites regularly distributed over the Earth's surface.However, forward numerical modeling indicates that if we monitor only the conductivity variation within a very limited area in the Earth (for instance, geothermal reservoir or magma chamber) it may be sufficient to interpret properly data measured even at one site (Spichak, 1999b).Since very few 3D MT surveys and interpretation experience are currently available in the world, only computer simulation of both MT field behavior and monitoring may provide a proper basis for effective 3D MT imaging and monitoring of active zones in the lithosphere.

Imaging volcanic interiors using MT data
MT sounding has been used to monitor volcanic activity and to help understand the processes leading to eruptions (Fitterman et al., 1988;Mogi and Nakama, 1990;etc.).Modeling, however, has been relatively crude.For example, Newman et al. (1985) modeled a homogeneous prism in a layered earth with a 3D integral equation method to study the detectability of a magma chamber, whereas Moroz et al. (1988) built a more elaborate analogue model to study the distortion of MT fields due to a volcanic cone.Mauriello and Patella (1999) studied the internal structure of the volcanic environment using probabilistic tomography.Spichak (1999a,b) studied the MT response in volcanic environments using a 3D geoelectri-al model of a Hawaiian type volcano.It is based on a model constructed in 1989 (with Prof. George Keller) to determine whether the MT method could also detect the internal proc-esses that had previously been found with TDEM measurements in the Kilauea volcano (Jackson and Keller, 1972).

Model of the volcano
The model represents a marine shield volca-no, characterized by a low and flat summit formed by homogeneous basaltic rocks (Fig. 1).Its flanks stretch down into the ocean; the conductivity of the ocean water was taken to be 3.6 S/m.The volcano's summit, 0.5 km thick, is formed by basaltic lavas with conductivity   By a high content of brine/seawater (this zone is 1.7 km thick and has a conductivity σ = 0.17 S/m).At 3 km depth from the volcano summit, there are dense lava formations 5.5 km thick with conductivity σ = 0.01 S/m, underlain by crystalline crust with conductivity σ = 0.001 S/m.
The conductivity distribution in the model was considered to be symmetrical with two ver-tical planes of symmetry.So, only onequarter of a 3D grid was used for calculations.

Synthetic MT pseudosections
MT fields for this model were synthesized by the program package FDM3D (Spichak, 1983) for two polarizations of the primary field at four periods T = 0,1; 1, 10 and 100 s.Then, a number of MT field transformations were cal-culated and analysed at different levels in the atmosphere to find those which are most sensi-tive to the parameters of the model.
Since the geoelectrical structure is strongly screened by sea water, direct estimation of the conductivity distribution from apparent resistivity pseudosections is very difficult.In contrast, 3D isosurfaces or 2D contour maps of the transforms based on the impedance phases and on the in-phase and quadrature parts of the horizontal electric fields are the best for imaging the structure of the volcano (Figs. 2 to 6).
In particular, Figs. 2 and 3 show the vertical cross-sections of the volcano overlapped by the contour maps of the transformed impedance phases ( φ xy and φ det , correspondingly), construct- ed for the plane at a height 0.5 km above the summit of the volcano.Although the values assigned to these isolines could hardly be interpreted in terms of the rock physical properties, their spatial gradients clearly indicate the location of the magma chamber and of the conductive formation above it.It is worth mentioning that the gradient of φ det delineates not only the magma chamber but also the flanks of the volcano (Fig. 3).This result agrees with earlier   findings by Park and Torres-Verdin (1988), Zhdanov and Spichak (1992), etc.) on the interpretation of the impedance phases.
Fig. 4 shows Re Z yx isolines.They indicate a position of the lower boundary of the magma chamber and even the interface between the dense volcanic rocks and crystalline crust.
Transforms of the in-phase and quadrature parts of the horizontal electric field component parallel to the incident electric field (Figs. 5 to 7) turn out to be even more sensitive to gradients of conductivity.Fig. 5 shows the vertical crosssection of the model overlapped by the map of the isolines of the Re y E transformation.The isolines condensation correlates with gradients of the conductivity with the local extremes marking the upper and lower edges of the magma chamber.
Isolines of the transformations of ImE y are shown in Fig. 6.There is a strong maximum located at the lower boundary of the magma chamber delineated by the isolines.A 3D pseudo-structure of the medium can be fairly well seen using an ImE y based volume image (Fig. 7).
Thus, construction of 3D pseudo-structures based on joint interpretation of transforms of the in-phase and quadrature parts of the horizontal electrical field and the impedance phases appears to be a useful tool for delineating the geometric parameters of the complex volcanic environments.This modeling supports the analytical findings of Szarka and Fischer (1989), explaining the behaviour of the MT field transformations at the Earth's surface in terms of the distribution of subsurface currents.

Methodology of interpretation of the MT data measured over the relief surface
The calculation of the transformations used in the previous section is very important in practice.They could be determined at the relief earth surface from the electric and magnetic fields measured there.However, multiple effects of the geological noise caused by near-surface inhomogeneities may greatly distort correspond-Fig.7. Im E y ( ) ×10 1 iso-surfaces for volcano model shown in Fig. 1.
ing MT field transformations and, consequently, the interpretation results obtained using some imaging technique.One way to overcome this difficulty consists of foregoing upward analytical continuation of the data to the artificial reference plane located in the non-conductive atmosphere higher than the top topographic point by means of the integral transformations of the Stratton-Chu type (Spichak, 1999a): ! H a are the anomalous electric and magnetic fields, correspondingly, determined at the reference plane, -is Green's function of free space, S is the surface of measurements, r !-is the radius-vector of the points belonging to this surface and n ! is a unit vector normal to the surface S and pointing outwards.
The analytical continuation of anomalous MT fields from the earth surface to the atmosphere could be proved to be unique (Berdichevsky and Zhdanov, 1984) and, therefore, continued data could be used for the interpretation.The corresponding numerical calculations are straightforward and stable.Another important item concerns the appropriate height of the reference plane.After many simulations it was found that it's optimal value is 250-500 m above the volcano summit.At lower heights, the pseudo-sections become distorted mostly by the nearest parts of the topography and by the noise, whereas at greater heights, details in geoelectrical structure may be lost.The choice of the optimal height for a given structure requires special investigations.
The synthetic MT fields and their transforms, calculated for a 3D model of a volcano of Hawaiian type have indicated that impedance phases and in-phase and quadrature parts of the electrical-field components are the most sensitive to the structure.Hence, the following procedure of volcano inner structure visualization can be outlined: 1) Upward analytical continuation of MTdata to an artificial reference plane located 250-500 m higher than the highest topographic point; 2) Construction of continued field transformations from a reference frame (X, Y, Z) to a new one (X, Y, Z app (or logT)); 3) Spline-interpolation of constructed onedimensional vertical profiles on a regular threedimensional grid; 4) Computer tomography of geoelectric structure by construction of transformation isosurfaces in the domain of search, limited horizontally by the size of the data array and vertically by the largest skin-depth.
The model studies show, that for a visualization of the volcanogenic zones it is necessary to use not only magnetic, but also electric fields.This makes the MT-method indispensable for studying their geoelectric structures.
Finally, the results obtained indicate that three-dimensional imaging of volcanic structure is possible not only on the basis of the EM data measured on its surface, but also at different levels, which are higher than the highest point of a survey relief.This, in turn, enables the remote sounding of regions of difficult access.

Geoelectric model of a typical volcano
In recent papers (Spichak, 1995;1999 a ) a number of 3D geoelectrical models of volcanoes were constructed and the effect of melt condition in a magma chamber on MT data was studied.The probable volcanic edifice is represented by a cone with a base diameter of 30 km and height 1.3 (5.0) km (Fig. 8).It was also supposed, that the cone is made of volcanic rocks with 1000 m ⋅ Ω  at a depth of 2.5 km from the Earth's surface there is a cubic magma chamber (5 x 5 x 5 km), filled with basalt melt of ρ = 2 m ⋅ Ω and connected with the crater by a channel with a diameter of 0.6 km.

Detection of the magma chamber by MT data
At first, the estimation of the MT data resolving power in the detection of the magma chamber was conducted for the case, in which the model has the magma channel (i.e., contrary to the case considered in the previous section).
For calculation of synthetic MT fields the software package FDM3D (Spichak, 1983) was used.The calculations were made for four periods: 1, 10, 100 and 1000 s.
As previously, the pseudo-sections of various field components and their transformations were constructed.In Fig. 9 shows the pseudosection of the impedance phase constructed by the procedure described in the previous section.It can be observed that in spite of the cone influence, all basic elements of the geoelectric structure (magma chamber, channel, cone and 1D layering) are adequately contoured by the isolines of this transformation.
Effective monitoring in cases, when the observations are made on a relief surface and data contains noise, requires, at least, foregoing modeling.After the preliminary 3D model is determined, it is necessary to recognize those com- ponents of MT-field (or its function), whose measuring (or evaluation) allows a reliable monitoring of the selected parameters.

Estimation of MT data resolving power with respect to the conductivity variations in the magma chamber
To understanding how the melt state in a magma chamber can be monitored using MT-data measured on a relief surface, a simulation was carried out as follows.
For the model shown in Fig. 8, the MT-field was first calculated for the conductivity of the magma chamber ranging from 0.1 to 0.5 S/m.Then they were continued upwards to an artificial reference plane allocated 0.5 km above the highest point of the topography.On this plane, the maps of isolines of differences in field component transformations, and also of apparent resistivity, corresponding to two different values of electrical conductivity of the magma were constructed.
Let us consider how the increase of the electrical conductivity in the magma chamber from 0.1 to 0.5 S/m influences the components of electric and magnetic fields.
In Fig. 10a,b the difference maps of normalized amplitudes of magnetic field component H x at the period T = 1 s are represented in the case of cone absence (fig.10a) and cone presence ( 10b).In both cases relative changes of the field amplitudes (up to 24 % for cone absence and up to 5 % for a cone with 3 .1 = h km) have local character and are manifested in the neighborhood of the magma vent.On both maps there exist as zones, in which the sign of change of the component amplitude coincides with the sign of change of the electrical conductivity in the chamber as zones, where they are opposite (in Fig. 10a,b domains with negative values of relative change).This is caused by the distinct character of the free-space attenuation of the anomalous field for different values of electrical conductivity in the chamber.
As the period increases (in Fig. 11a,b the difference maps for T = 10 s are represented),   The results for the most expressive MT-field components in TE -mode were demonstrated.With change of the primary field polarization the configuration of isolines of differences (in particular, the location of domains, in which the sign of differences is the same) also changes.Unfortunately, interpretation in terms of the apparent resistivity components does not improve this situation, since appropriate difference maps have similar peculiarities.
From these results it is evident that for a reliable interpretation of monitoring of the electrical conductivity variation in the magma chamber it is desirable to use some scalar function of the MT-field, which takes into account field components measured at different polarizations of the external field.One (and, perhaps, unique) such function is «an apparent resistivity» (or «an apparent conductivity») based on the impedance determinant.Indeed, it naturally takes into account all horizontal components of the MT-field measured at different polarizations of the primary field, and its change has mainly the same sign as the variation of the resistivity (conductivity) in a magma chamber (as it will be shown below).
In Fig. 13a,b the maps of isolines of differences of normalized values det ρ at the period 10 = T s are represented: fig.13a for the cone absence and 13b for the cone presence.Irrespective of availability or not of the cone, the difference in values of the selected function in both cases has the same sign as the corresponding difference of resistivities in a magma chamber.Moreover, its values depend only on the horizontal distance to the center of the «crater» (at least, in a zone with diameter approximately three times larger than the horizontal sizes of the magma chamber) that can essentially simplify interpretation of monitoring results.
In order to eliminate the influence of the magma channel on the behavior of the isolines in a difference map, the same model without the channel was considered.In fig.14a,b the difference maps of ρ det for the model without the cone at the period 10 s corresponding to increase in conductivity in the magma chamber from 0.1 S/m to 0.2 S/m (fig.14a), and to 0.5 S/m (fig.14b) are represented.The analysis of the maps represented in fig.14a,b indicates: a) The maximum effect is observed in an area limited by the horizontal sizes of the mag-   Matching of fig.14b with fig.13a shows that the magma channel strengthens the effect of variation of magma melt conductivity ap-Proximately by three times.Firstly, it increases the diameter of the zone for reliable monitoring, and, secondly, reduces the period threshold sufficient for detection of even small variations of electrical conductivity in the magma chamber.

Procedure for MT-monitoring electric conductivity in a magma chamber
From the preceding sections it follows that for a reliable monitoring of a melt state variations it is worth following the next procedure: 1) To construct a geoelectric model of the region; 2) To create difference maps («templates») similar to those given in fig.14a,b; 3) To locate the sensors at a distance from the estimated centre of the magma chamber no more than three times exceeding its horizontal diameter; 4) To make measurements, whenever possible, of long-period field variations and at different primary field polarizations; 5) To continue analytically the measured MT-fields upwards to a reference plane (by the procedure described in Section 2.3); 6) To create difference maps by means of calculation of ρ det from the continued MT fields at the reference plane; 7) To find relevant variations of melt resistivity using «templates».
The use of the upward analytical continuation allows us not only to filter the data, as in the case of imaging, but also to reduce the data «to a common denominator», which is important for removing the relief effect on the transforms «templates» and, consequently, for proper interpretation of the monitoring data (at least, in the framework of the procedure suggested above).
The modeling results indicate that the apparent resistivity function derived from the impedance determinant, is most suitable for adequate interpretation of measurements carried out for the purpose of monitoring.
Lastly, the results shown in this section point to the opportunity of remote monitoring of volcanoes and other targets of difficult access.

Conclusions
Synthetic MT fields and their transforms, calculated for 3D geoelectric models of volcanoes, indicate that the impedance phases as well as in-phase and quadrature parts of the electrical-field components provide the best imaging of the volcanic interior.
On the other hand, the apparent resistivity function derived from the impedance tensor determinant seems to be the most suitable parameter for adequate interpretation of measurements carried out for the purpose of monitoring.
The use of the upward analytical continuation of MT data to the artificial reference plane allows us not only to filter the data necessary for proper imaging, but also to reduce the data «to a common denominator», which is important for removing the relief effect on the transforms' "templates" and, consequently, for robust interpretation of the monitoring data (at least, in the framework of the procedure suggested in Section 3.4).
The use of the interpretation methodologies proposed in this paper should make it possible to monitor the conductivity variations in some very restricted zones of the Earth's crust (in particular, in the magma chamber) based on appropriate measurements carried out even in one site properly located with respect to the position of the target.
It is obvious that the actual situations are much more complicated than the models considered in this paper.However, the methodologies suggested following the results of computer simulation may help, first, to avoid unnecessary interpretation errors due to noise (both geological and instrumental) in the data, irregularity and lack of observation points in regions of difficult access; and, second, to make MT survey planning more scientifically substantiated.I am grateful to the staff of the Laboratory of Electromagnetics (Dipartimento di Scienze Fisiche, Universita' di Napoli "Federico II") headed by Prof. Domenico Patella and of the Osservatorio Vesuviano (Director Prof. Lucia Civetta) who supported my research activity during my visit to Naples and Prof. Giovanni Orsi who kindly offered to me the computer facilities Finally, many thanks are addressed to Dr. Tzanis and anonymous referee for valuable comments that improved the text of the manuscript.

Fig. 5 .
Fig. 5. Re E y ( ) ×10 1 pseudo-section for volcano model shown in Fig.1.Result is based on an incident electric field, parallel to O y -axis.

Fig. 6 .
Fig. 6.Im E y ( ) ×10 1 pseudo-section for volcano model shown in Fig.1.Result is based on an incident electric field, parallel to O y -axis.
resistivity.Below there is 2.5 km thick layer of volcano-sedimentary rocks with resistivity ρ = 20 m ⋅ Ω .Crystalline rocks of the crust bottom with 1000 m ⋅ Ω resistivity and 45-50 km of thickness underlie it.Finally,

Fig. 8 .
Fig. 8. Geoelectric model of a typical volcano (vertical cross-section under the central profile).The figures show electrical conductivity values in S/m.

Fig
Fig. 10a,b.Isolines of normalized amplitude differences (in %) of the component H x at a period T = 1 s: a)in the case of the cone absence and b) -at its presence.

Fig.
Fig. 11a,b.Maps of difference isolines (in %) of normalized amplitudes of magnetic field component H x at the period T=10 s: a) in the case of the cone absence and b) at its presence.

Fig. 12 .
Fig. 12. Maps of difference isolines (in %) of normalized amplitudes of electric field component E y at the period T=10 s: a) in the case of the cone absence, b) at its presence.

Fig
Fig. 13a,b.Map of difference isolines (in %) of normalized values ρ det at the period T=10 s a) in the case of the cone absence and b) at its presence.

Fig.
Fig. 14a,b.Map of differences isolines (in %) of normalized values ρ det at the period T = 10 s for the model without the magma channel shown in a Fig. 8 at decrease of resistivity in the magma chamber: a) in 2 times; b) in 5 times.
(a) and (b), correspondingly) times exceeding the horizontal sizes of the magma chamber.