A realistic inversion algorithm for magnetic anomaly data : the Mt . Amiata volcano test

The aim of this work is the formulation of a 3D model of the Mt. Amiata volcanic complex (Southern Tuscany) by means of geomagnetic data. This work is shown not only as a real test to check the validity of the inversion algorithm, but also to add information about the structure of the volcanic complex. First, we outline briefly the theory of geomagnetic data inversion and we introduce the approach adopted. Then we show the 3D model of the Amiata volcano built from the inversion, and we compare it with the available geological information. The most important consideration regards the surface distribution of the magnetization that is in good agreement with rock samples from this area. Moreover, the recovered model orientation recall the extension of the lava flows, and as a last proof of validity, the source appears to be contained inside of the topographic contour level. The credibility of the inversion procedure drives the interpretation even for the deepest part of the volcano. The geomagnetic signal appears suppressed at a depth of about 2 km, but the most striking consequence is that sub-vertical structures are found even in different positions from the conduits shown in the geologic sections. The results are thus in good agreement with the information obtained from other data, but showing features that had not been identified, stressing the informative power of the geomagnetic signal when a meaningful inversion algorithm is used. Mailing address: Dr. Osvaldo Faggioni, Istituto di Geofisica Marina, Università di Pisa, Consorzio Universitario della Spezia, Via Vailunga 1, 19125 La Spezia, Italy; e-mail: of@castagna.it


Introduction
Geophysical data, and above all geomagnetic data, are particularly useful in volcanology to define the characteristics of a volcanic complex.A geomagnetic data set is very informative, in fact, if connected to magmatic rocks or lytotypes that show a particular intense value of the magnetic susceptibility.
Using adequate inversion software, one can find solutions that seem to correspond to the tectonic-structural model of the subsurface, but additional information can be obtained over areas that cannot be investigated directly.Practically, the geomagnetic data are used to find solution about an equivalent source that should be responsible for the observed data, through the inversion process.The transfer from the observed data to a meaningful set of information on the generating source is not obvious, and has a fundamental relevance.This is particularly true for the geomagnetic field, that because of its di-polar nature is far from being interpretable «by eye».
It is well known moreover that any inversion process suffers from non-uniqueness problems.We know from Gauss' theorem that when the magnetic field is known only on a bounding surface, there exists an infinite number of source distributions inside the volume that can generate this field.We call this a theoretical ambiguity.However, another kind of ambiguity has an algebraic origin: the number of available data is often under-sampled so that the description is not complete to fully describe the field.Many different source distributions can reproduce a particular discrete data set with the same degree of accuracy.In practice, it can easily be shown that any magnetic anomaly measured on the surface of Earth, might be reproduced by considering a superposition of (even negative) magnetic dipoles in a thin layer just beneath the surface.Today there are two principal strategies in order to develop an algorithm able to overcome these ambiguities.One possibility is given by parametric inversion.Knowing the anomalies generated by a simple causative body, it is possible to approximate the source finding the «best one» (Bhattacharyya, 1980;Wang and Hansen, 1990;Zeyen and Pous, 1991).Non-uniqueness is avoided by restricting the number of degrees of freedom to be fitted.In the second methodology, the Earth is divided into a large number of prismatic cells with constant magnetization inside and an algorithm minimizes some functional representing a suitable objective function of the model.For example, it is possible to avoid dispersion of the source by minimizing its total volume (Last and Kubik, 1983) or its momentum of inertia (Guillen and Menichetti, 1984).Realistic models are obtained contrasting the depth decay of the kernel function with an appropriate weighting function (Li and Oldenburg, 1996) or introducing a multi-layer depth resolution in field measurements (Fedi and Rapolla, 1999).In this work, the result is achieved using an inversion scheme that has been developed by Caratori Tontini et al. (2003) and in the following section an outline of this inversion scheme will be shown, to provide the formal context of the performed analysis.The data used for the inversion procedure come from the ENI-AGIP Division, and are used according to the scientific collaboration between ENI and the Istituto di Geofi sica Marina of La Spezia (Eni Explor. & Prod. Div. and IGMar, 2002).The geologic map and sections used for the interpretation come from the work by Enel «Direzione Studi e Ricerche, Centro di Ricerca Geotermica» and CNR «Istituto Internazionale per le Ricerche Geotermiche» (Calamai et al., 1970;Gianelli et al., 1988).
Figure 1 shows the topographic contour level of the Mt.Amiata complex projected onto a Gauss-Boaga reference system.The Mt. Amiata shows two major peaks, La Vetta (1735 m), and La Montagnola (1571 m) aligned along a WSW-ENE direction.It is a moderate size silicic volcanic complex, whose main activity had two principal phases at about 300 000 and 200 000 years ago (Ferrari et al., 1996).The first period is characterized by the emplacement of a Trachi-Dacitic Basal complex (BTC), made of two separate sub-units (fig.3).The whole BTC complex lies over a 90 km 2 area with an average thickness of about 150-200 m, and probably it is the result of a light explosive activity.The second phase happened after 100 000 years and its main products were domes and lavic flows with composition variable from Trachi-Dacitic, Trachitic and Latitic (DLC).These domes and flows, where visible, are aligned along a WSW-ENE direction as the BTC complex.At the end of the eruptive epoch, two smaller Olivine-Latitic Lava flows (OLL) were emitted.
The detection of fault systems and their origin in this area is difficult because of the vast erosional cover surrounding the domes and flows and for the dense vegetation.However the Mt.Amiata is oriented WSW-ENE and two faults, with a length of about 2.5 km, seems to be the origin of the BTC collocation.A certain number of other smaller faults were activated during the second period of activity (Marinelli et al., 1993).

Outline of the inversion algorithm
In this section, we show briefly the method used to perform the inversion, its constraints and the way it works.A full treatment of the topic can be found in Caratori Tontini et al. (2003), together with some tests, based on both synthetic and real data that gives reliability to the method.
Here we report the formalism and the basic information about the initialization and the utilization of the algorithm, to fully understand the inversion of the Mt.Amiata geomagnetic data.Let us suppose that we have partitioned the data set into a 2D grid that we call d(i, j) (observations).We want to interpret this data set in terms of the magnetization of the generating source.We adopt a mapping of the region presumably containing the source, dividing it into a certain number (n x , n y , n z ) of prismatic cells with unitary magnetization.Obviously, the magnetic anomaly field generated by a particular subsurface configuration is obtained multiplying the unitary cell contribution for the magnetization of the cell, and summing over the whole number of cells (2.1) where DF(i, j) is the grid of the predicted data and M(l, m, n) is the cell magnetization.The ) ∑ analytic expression of the contribution at one cell ∆ ˜, , , , )is well-known and here is not reported for conciseness, but can be found in Bhattacharyya (1964); Rao and Babu (1991).
The inversion procedure works minimizing the difference between the observed data set d(i, j) and the prediction grid DF(i, j) in the least squares sense, provided that some constraints over the magnetization are imposed, otherwise it is well-known that the problem has multiple solutions.
The algorithm adopts a smoothness constraint and the bound of positivity imposed by using a Gaussian approximation.Practically the source is expanded using a predetermined number of 3D Gaussian functions that are positively constrained.This particular choice introduces a non-linear relationship among the parameters defining the source, but has the fundamental relevance of drastically reducing the degrees of freedom of the problem, permitting however the reconstruction of the most articulated sources.The smoothness constraint moreover is imposed naturally in this way.Within this choice, the magnetization is parameterized as where the set of 7 parameters (A k , x k , y k , z k , sx k , sy k , sz k ) define completely the k-th Gaussian function, while the number n g of functions used as the basis of the envelope is chosen by the user according to some empirical considerations.If the algorithm starts with some unnecessary function, i.e. n g is redundant, the amplitude of these superfluous functions is compelled to a null value through a process of annihilation that is explained in detail by Caratori Tontini et al. (2003).
It is a valid geological observation that often the spatial variations of physical properties are smooth.This is particularly true for some geological scenarios: for example smooth magnetization variations can be seen near magmatic intrusions, and smooth density variations may apply to sedimentary formations under different level of compaction.The positivity constraint, that is of fundamental relevance to reduce the ambiguity domain, in the case of geomagnetic data is obvious because the paramagnetic effects of matter are usually stronger than the diamagnetic ones.This further limitation is imposed through the use of a positivity mapping function for the amplitude parameters A k (Oldenburg and Li, 1994;Li and Oldenburg, 1996).The algorithm proceeds by means of a modified Levenberg-Marquardt iteration scheme for non-linear least squares minimization to determine the best-fit source, together with the statistical confidence level for the parameters.Depending on the number n g of functions involved, even the most composite models can be reconstructed by the inversion algorithm to achieve a fine resolution.

The Mt. Amiata aeromagnetic data inversion
The data used for the inversion were collected during the survey performed by AGIP in 1978-1979 in a large area, including Tuscany and the Tyrrhenian sea, using a CFS Cesium magnetometer.The upper part of fig. 2 shows the magnetic anomaly using Gauss-Boaga kilometric coordinates reference system.The total field data (time corrected) were levelled and filtered through standard procedures, and then a polynomial regional field (order 2) was residuated to obtain the magnetic anomaly.The survey in this area was made through lines (direction SW-NE) with an average spacing of 2 km, orthogonally intersected by tie lines (direction NW-SE) with an average spacing of 5 km.The sampling pass along each profile was about 500 m, at an altitude of 1600 m.
The data were gridded through a minimum curvature method with a cell size of 750 m, for a spectral resolution of 1.5 km.The magnetic anomaly is very articulated, showing different positive local maxima with the associated negative lobes, suggesting the idea of a very articulated source for geometry and depth.For these reasons, these 22 ¥ 16 data are inverted using a prismatic region defined by the east interval (1708; 1716) km and north (4747; 4755) km with 4 km depth, divided into 6912 cubic  volcanic edifice.We suppose that the highest magnetized recovered bodies belong to the DLC and OLL volcanic complexes.Two interpretative geological sections (Calamai et al., 1970;Giannelli et al., 1988) show two main conduits located in proximity of «La Vetta» and «La Montagnola» that represent the locations of the shallower portion of magma ascent during the emplacement of domes structures.Probably those vertical structures reach a depth of 5-7 km, where a magmatic intrusive granitic body (Giannelli et al.,1997) that supplied the dome pile up is supposed to be located.A magnetic signature is given only by the upper portion of those vertical structures, because of the very shallow Curie isotherm of the Amiata geothermal region (Baldi et al., 1995).
In fig. 3 the recovered model has been spread over the topographic surface for comparison with the rock magnetization data sampled during field measurements that the Istituto di Geofisica Marina performed in the area of the Mt.Amiata, using a GMS-2 susceptibility meter.The magnetic susceptibility of different lytotypes coming from rock samples of BTC, DLC and OLL was measured, and an average magnetization was calculated from the different samples to avoid statistical errors coming from geometrical irregularities of the sample surface.The first three horizontal cross-sections of the recovered model have been fused into a unique map, in order to overlap the data found during the inversion procedure on the topographic surface of Mt.Amiata.The high magnetic behaviour of the rocks coming from DLC and OLL, with respect to the lowly magnetized complex of BTC, is confirmed by the recovered model.The model shows the presence of the generating source in the central part of the area as arguable from the geologic map.The source geometry appears aligned along the WSW-ENE direction in good agreement with the orientation of the principal fractures that caused the emplacement of the volcanic complex; moreover the magnetization peaks seems to confirm the presence of localized sources that should be referred to the volcanic conduits system mentioned before.The inversion has highlighted that the magnetization of subvertical sources located inside Amiata volcanic edifice is stronger than that measured on the surface; the deepest part of the domes that have a surface topography expression like «La Vetta» and «La Montagnola», show a higher degree of magnetization before it reaches the shallow Curie isotherm: this is a new fundamental constraint to be considered to understand the petrographic evolution of the Amiata volcanic system, imposed only by the recovered inversion model from magnetic anomaly data.

Conclusions
From the data obtained by ENI-AGIP Division we have evaluated the magnetic anomaly of the Mt.Amiata volcanic complex.From the inversion of the extracted data set we have shown a 3D model of the magnetic source inside the apparatus.The model shows the geometry of the magnetic sources of the volcanic edifice from the altitude of 1.5 km to the depth of 2.2 km, where the geomagnetic signal loses its informative capability.The recovered model is in good agreement with the information coming from the geologic knowledge and the susceptibility measurements of the area of interest, showing however characteristics and features that were not still found.In particular, we were able to suppose the presence of sub vertical structures that should be connected to intrusions or localized magma ascent conduits characterized by a magnetization greater than the surface portion one.The surface distribution of the magnetization is in good agreement with the geologic map and the rock samples measurements.This study is thus fundamental in the definition of the geologic submerged structures of the Mt.Amiata apparatus, otherwise not accessible by direct methods.

Fig. 2 .
Fig. 2. Magnetic anomaly and recovered model compared with the topographic contour level (Gauss-Boaga projection).The broken lines represent the geologic sections of fig. 4.

Fig. 3 .
Fig. 3. Surface distribution of the recovered source compared with the geologic map of the Mt.Amiata area and the geologic sections (Gauss-Boaga projection).

Fig. 4 .
Fig. 4. Geologic sections of the area of the Mt.Amiata.