November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

INVERSE MODELING IN GEOPHYSICAL APPLICATIONS

G. CURRENTI ,R. NAPOLI, D. CARBONE, C. DEL NEGRO, G. GANCI Istituto Nazionale di Geosica e Vulcanologia, Sezione di Catania, Catania, Italy E-mail: [email protected] maglab.ct.ingv.it The interpretation of the potential eld data is an useful tool that allows for both investigating the subsurface structures and providing a quantitative evaluation of the geophysical process preceding and accompanying period of volcanic unrest. Potential eld inversion problem are required to combine forward models with appropriate optimization algorithms and automatically nd the best set of parameters that well matches the available observations. Indeed, investigations on the mathematical equations to be inverted, have revealed that these models are ill-posed and highly non-linear. Numerical methods for modeling potential eld observations are proposed and applied on real dataset. Keywords: Geophysical inversion; Potential Field; Mount Etna.

1. Introducton The inversion problem in potential eld modeling suers from the ambiguity and instability of its solutions. The ambiguity arises from the inherent property of potential elds for which dierent combinations of parameters may lead to similar anomalies. Moreover, potential elds inversion is notoriously unstable. Because of ambiguity and instability of solutions, the inversion problem can be secured by narrowing the set of all possible solutions to a predened solution class that allow for a unique and stable solution. The priori denition of the geometry (simplied bodies: spherical, rectangular, prismatic) of searched source and the priori recognition of the involved physical mechanism allows reducing the number of likely solutions considerably.1 Potential eld inversion method can be classied into two main categories depending on the type of the unknown parameters to be retrieved. Firstly, there are some methods looking for magnetization or density contrast values with xed geometry sources. In this case linear inversion technique can be successfully applied. Secondly, there are those

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

methods that look for geometrical determinations of the source. In this latter case the inversion problem is highly non-linear and robust nonlinear inversion techniques are needed. The rst category of interpretation models is widely used in providing physical property distributions of subsurface geological structures, while the second one is generally applied in modeling potential eld sources in volcanic areas. In the following two inversion approach are described to deal with linear and non-linear inversion problems and are applied on two case studies to appraise the goodness of the proposed methods. 2. 3D linear inversion of potential eld data Potential eld inversion are aimed to provide a model which reconstructs, as well as possible, subsurface geological structures having a density/magnetization contrast with their surroundings. For the sake of simplicity, we consider the magnetic inverse problem. However, similar formulation can be applied to the case of gravity inversion method. The computational domain V which is supposed to surround the magnetic source is discretized using a nite number of m = Nx Ny Nz rectangular prisms whose magnetization Jj is uniform inside each prism.1 Using a 3D discrete numerical approach, the total anomaly eld at i th observation point is computed by:
n

Ti =
j=1

Ji aij

(1)

where elements aij quantify the contribution to the total magnetic eld anomaly T at ith observation point due to the magnetization of the j th prism.2 Therefore, the inverse problem can be formulated as the solution of a system of n linear equations as: Ax = T (2)

where x is the m vector of unknown magnetization values of the prisms, T is the n vector of observed magnetic data, and A is a matrix with elements aij . The analytical expression of the aij term for a prismatic body was devised by Rao and Babu.3 Based on discretization, the number of prisms is usually larger than the number of observation points, thus the linear inverse problem in Eq. (2) turns out to be unavoidably underdetermined. In such a case, the linear system leads to a solution with m n degree of

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

freedom. A further diculty in solving the system in Eq. (2) is due to the inherent non-uniqueness of the potential eld: there are an innite number of inverse models that can explain the same observed magnetic anomaly within error limits. This highlights that the magnetic inverse problem is ill-posed. That calls for some regularization techniques. In such a case, it is necessary to impose further constraints taking into account a priori knowledge about the solution when available. The idea of reducing the class of possible solution to some set on which the solution is stable lies on the fundamental concept of introducing a regularizing operator. The inverse problem can be re-formulated as an optimization problem aimed at nding the unknown magnetization values x that minimize a functional composed of a data mist d and a smoothing functional m : 1 T [(Ax Tobs )(Ax Tobs ) + (x x0 )T W T W (x x0 )] (3) 2

= d + m =

where W is a weighting matrix and is the regularization parameter. Because of the lack of depth resolution, W is aimed to counterbalance the contribution of deeper prisms with respect to shallower ones. To ensure that the solution is geologically reasonable it is advisable to prescribe realistic bounds on the magnetization values on the basis of rock samples or available information about the local geology. The minimization of the quadratic functional in Eq. (3) subjected to bound constraints can be solved by using a Quadratic Programming (QP) algorithm based on active set strategy:4 1 min = min[ xT Qx dT x], L <= x <= U 2

(4)

where Q = AT A + W T W and d = AT T + W T W x0 , and L and U are the vectors of lower and upper bounds. The quadratic formulation of the problem is solved iteratively by generating a sequence of feasible solutions that converge toward the optimal solution. The iteration is stopped when no relative improvements in the functional are achieved. The inversion procedure described above was applied to analyze the anomalies detected by a ground magnetic survey of the Ustica island covering an area of about 9 km2 . The total-intensity anomaly eld, obtained after data reduction process, shows the presence of a W-E striking magnetic anomaly in the middle of the island and other two intense anomalies, which seem to continue oshore, in the south-western and the north-eastern sides (Fig. 1).

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

Fig. 1.

Map of the total intensity magnetic eld after reduction process.

In order to allow the maximum exibility for the model to represent geologically realistic structures the island was represented as a crustal block, 4x3 km2 in area and about 1.2 km in thickness, and was discretized into a set of rectangular prisms (0.125 x 0.125 x 0.15 km3 in size). The anomaly eld was inverted assuming that the average direction of total magnetization is close to Earths present-day eld direction (approximate inclination of 54o N, declination of 1o W). Moreover taking advantage of magnetizations measured on volcanic rocks of the island, the inversion of the data set was performed constraining the range of the magnetization value for the assumed sources from 0 to 10 A/m. The iterations were continued until the functional does not show signicant improvements. The residual eld, the dierence between observed eld and calcualte eld is shown in Fig. 2. Only in the N-E and S-W area of the island two small residual anomalies remain probably due to a low resolution of the volume discretization in these zones. The 3D model resulting from the constrained inversion of data set reveals the presence of three main magnetic trends, N-S, E-W and NE-SW, which are in good agreement with the surface geology, and are coincident with the main regional structural lineaments (Fig. 3). In particular, the N-S and E-W trends, which are of more recent origin,

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

N[m]

E[m]
Fig. 2. Residual, total-eld magnetic anomaly map produced by subtracting the computed eld from observed eld.

prevail in the shallow part of the model, while the NE-SW is relevant below 0.4 km b.s.l. At this depth the model reveals a magnetized body in the central area, which may be interpreted as the preferential area for magma storage and ascent and which supplied the feeding systems of the main subaerial volcanic centres of the island, Mt. Guardia dei Turchi and Mt. Costa del Fallo. Two other magnetized volumes were identied and ascribed to the small submarine/subaerial eruptive centres of the western island and to the younger cone of Capo Falconiera, respectively. These ndings highlight how the regional tectonics has strongly inuenced the structural and magmatic evolution of the Ustica volcanic complex producing preferential ways for magma ascents. Moreover,considering that the NE-SW trend prevails only in the deeper part of the model and is replaced at the shallow depths by the N-S and E-W orientations, it is possible to assume that a change of tectonic style and/or its orientation took place in the past5 .

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

4287000

4287000
N

4286000

4286000

4285000

4285000

4284000 339500 340500 341500

r Tyr

ia hen

Sea
4284000 Z= - 300 m 340500 341500 342500 343500

342500

343500

339500

4287000

4287000

4286000

4286000

4285000

4285000

4284000 339500 340500 341500

Z= - 500 m

4284000 343500 339500 340500 341500

Z= - 700 m

342500

342500

343500

4287000

4287000

4286000

4286000

4285000

4285000

4284000 339500 340500 341500

Z= - 900 m

4284000 343500 339500 340500 341500

Z= - 1100 m

342500

342500

343500

Am-1

Fig. 3. The 3 D magnetization model of the Ustica volcanic complex. Horizontal sections of the uppermost part (above sea level) of the model and from 300m depth until 1100m depth.

3. Non-linear inversion of gravity anomaly by genetic algorithm technique Attempts to model potential elds expected to preceed and accompany volcanic eruptions often involve a great deal of eort due to the complexity

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

of the considered problem. The inversion problem deals with the identication of the parameters of a likely volcanic source that leads observable changes in potential eld data recorded in volcanic areas. Indeed, investigations on the analytical models describing the involved geophysical processes, have revealed that the models are highly non-linear and, usually, characterized by several parameters. When non-linear models are involved, the inverse problem becomes dicult to solve through local optimization methods. Hence, elaborated inversion algorithms have to be implemented to eciently identify the source parameters. We have investigated the use of Genetic Algorithms6 (GAs) which perform a broad search over the parameter space using a random process with the aim of minimizing an objective function that quanties the mist between model values and observations. The GAs inversion strategy was applied on a gravity anomaly that grew up progressively in 5 months before the 2001 ank eruption of Mt. Etna along a East-West prole of stations on the southern slope of the volcano.7 Between January and July 2001, the amplitude of the gravity change reached 80 Gal, while the wavelength of the anomaly was of the order of 15 km (Fig. 4). Elevation changes observed through GPS measurements during a period spanning the 5-months gravity decrease, remained within 4-6 cm all over the volcano and within 2-4 cm in the zone covered by the microgravity prole. Since February 2001, an increase in the seismicity was observed,8 with many earthquake locations clustered within a volume at a depth of about 4 km bsl and focal mechanisms pointing towards a prevailing tensile component. Therefore, we review both gravity and elevation changes by the Okubo model9 which mathematically describes the eect of uniformly distributed tensile cracks on gravity and ground deformation eld in an elastic homogeneous half-space medium.10 The gravity solution consists of four contributions: (i) the free-air eect proportional to the uplift, (ii) the Bouguer change caused by the upheaved portion of the ground, (iii) the gravitational attraction of the crack-lling matter and (iv) the gravity eld due to the redistribution of mass associated with the elastic displacements. This model implies inversion for nine parameters, m = (A, X, Y, Z1 , Z2 , L, W, U, ), whose description is reported in Table 1. To make the GA converge towards a solution which (i) best ts the observed data and (ii) has a good chance to be realistic from the volcanological point of view, we suitably restrict the parameter space to be investigated. The ranges of variability for the model parameters to be found are set on the grounds of the available geophysical and geological evidence. For the inversion procedure adopted in

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

3 Km 4

NEC VOR BN SEC

3000

1.8
2000
20 00

N [m]
4185000
2.7 4.4 2.3 4.1

2550 m vent

Valle del Bove

1.8 5.8

Summit -1.5 Craters

6.5

-1.5

4180000

2.7
Legend Pyroclastic cones 2001 Lava Flows Eruptive fissure Fracture

Valle del Bo ve

4175000

B
1.6 2.2

C
4170000

A
100 0

Gravity Station GPS Station

0 5 km

4165000

Dg [mGal]

-20 -40 -60 -80 485000 490000 495000 500000 505000 510000

E [m]
Fig. 4. Sketch map of Mt. Etna. Top: the gravity stations of the microgravity network along the East-West Prole (ABC) and the GPS stations measured in July 2000 and June 2001. The elevation changes higher than 1.5 cm, observed during the same interval, are also reported. The inset on the left shows the 2001 lava ow. Bottom: gravity change observed between January and July 2001 along the East-West Prole.

the present study, we have set the objective function equal to the 2 value that accounts for the measurements error dened by the standard deviation , as:

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

2 =
k

(Mk Ck )2 2 k

(5)

where Mk are the measured data, Ck are the computed gravity variations and k is the number of available measurements. The parameters of the best model found through the GA are reported in Table 1, together with the edge values of each search range.
Table 1. The best model parameters found by the GA. Minimum 0 1000 3000 500 -45 4170000 495000 -2400 0 Maximum 2500 5000 7000 2000 0 4180000 505000 -2700 2 Best value 409 5000 7000 500 -30 4173940 503499 -2500 2

Parameter Z1 - Depth of the top (m b.s.l.) L - Length (m) H - Height (m) W - Thickness (m) A - Azimuth (from the North) X - Northing Coordinate(m) Y - Easting Coordnate(m) - Density contrast (kg/m3 ) U - Extension (m)

The t between observed and calculated gravity changes is very good (Fig. 5), with a residual of 7.43 Gal, well within the error on temporal gravity dierences along the East-West prole (10 Gal at the 95% condence interval). Some of the best values found by the GA correspond to one of the edge value of the corresponding search range. This is in principle not ideal since it states for the GA to have been unsuitably conned during its search procedure. A sensitivity analysis of the parameters whose best value coincides with an edge of the chosen range is accomplished. Sensitivity tests are computed for the L-H and W-U couples of parameters, while the other remaining parameters are held xed at their best values. The tests highlighted that while U is a highly sensitive parameter, W appears not to be so much sensitive. This observation reects the assumptions behind the model: U denes the volume increase (volume of new voids, if the density of the lling material is set to zero) predicted by the model, while W only changes the volume within which the new fractures are uniformly distributed. In principle, changes of W do not inuence the net eect of the model-body, if W remains much smaller than the depth of the mass center of the body itself. As for the L and H, they are sensitive parameters and the obtained best values appears to be a good compromise between the needs of (a) not

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

10

rising the value of these parameters towards an implausible extent and (b) assessing a satisfactory t. In conclusion, the sensitivity tests conrm that the GA accomplished its task satisfactorily. Results show that, although it is possible to explain the observed gravity changes by means of the proposed analytical formulation, calculated elevation changes are signicantly higher than observed. This nding could imply that another mechanism, allowing a signicant density decrease without signicant deformation, coupled the tensile mechanism due to the formation of new cracks, increasing the negative gravity eect while leaving the displacements at the surface unchanged. The only mechanism allowing a density decrease at depth without surface displacements is the loss of mass from a magma reservoir. To keep the maximum elevation change within 2-4 cm, a value of extension U less than 1 m should be assumed. Using this value as an upper bound for the extension of the fracture zone, only about 50% of the observed gravity decrease can be explained using the Okubo model. Under the assumption that the new-forming tensile cracks and the loss of mass come from the same inferred source, it results that an about 1011 kg mass should be lost to contribute the missing 50% of the observed gravity decrease. The estimated mass corresponds to a minimum magma volume of 35*106 m3 . It is worth to note that, during late 2000 and the rst months of 2001, an almost continuous activity, with lava emission and Strombolian explosions from the summit Southeast Crater (SE), was observed.11 The volume of the products emitted between January 2001 and the start of the main ank eruption is estimated to range between 13 and 20*106 m3 . On the basis of our calculation, the volume emitted from the SE crater is of the same order of magnitude as the volume lost from the estimated gravity source. The hypothesis that the volume emitted from the summit SE crater, during the January-June 2001 period, was supplied by the same source in which tensile cracks formed needs to be further investigated to understand (i) which phenomenon could have triggered the mass transfer from the deep source volume to the surface and (ii) how this phenomenon relates to the formation of tensile fractures. 4. Conclusion We illustrate two procedures dealing with linear and non-linear inversion of potential eld data. Although the inversion problem suers from the ambiguity and instability of its solutions, numerical methods allow for narrowing the set of all possible solutions and providing a unique and stable solution. As for large-scale linear inverse problem, a 3D inversion of magnetic eld

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

11

Gravity Stations
GPS Stations
0 5 km

1.8 2.7 4.4 2.3

N [m]
4185000

N
1.8 5.8

4.1 -1.5 -1.5

4180000

10

6.5 2.7

-5
2.2

4175000

A
5

1.6

C
0
4170000

10

4165000

Dg [mGal]

-20 -40 -60 -80 485000 490000 495000 500000 505000

Computed Observed

510000

E [m]
Fig. 5. Results based on the best model. Top: contour map (at 5 cm intervals) of computed elevation changes; surface projection of the source (dashed line); observed (July 2000 - June 2001) elevation changes higher than 1.5 cm. Bottom: observed (January - July 2001) and computed gravity anomalies.

data was performed by means of QP algorithm to produce a magnetization model that provides useful information of the subsurface geological structure of Ustica volcanic complex. As for non-linear potential eld inversion, a GA techinique was proposed to infer the parameters of a volcanic source that justies the growth of a negative gravity anomaly preceding the 2001 Etna eruption. Our ndings demonstrate that the identication and interpretation of potential eld data can be a useful instrument both for de-

November 28, 2006

21:33

WSPC - Proceedings Trim Size: 9in x 6in

Currenti

12

tecting subsurface geological structure model and improving the monitoring of active volcanoes. References
1. R. J. Blakely, Potential Theory in Gravity and Magnetic Applications (Cambridge University Press, New York, 1995). 2. M. Fedi and A. Rapolla, Geophysics 64, 452 (1999). 3. D. B. Rao and N. R. Babu, Geophysics 56, 1729 (1991). 4. P. Gill, W. Murray, D. Ponceleon and M. Saunders, Solving reduced KKT systems in barrier methods for linear and quadratic programming, Technical Report SOL 917, Stanford University (Stanford, CA, 1991). 5. R. Napoli, G. Currenti and C. Del Negro, Bull. Volc. in print (2006). 6. G. Currenti, C. Del Negro and G. Nunnari, Geophys. J. Int. 162, 1 (2005). 7. G. Carbone, D.and Budetta and F. Greco, J. Geoph. Res. 108, 2556 (2003). 8. A. Bonaccorso, S. DAmico, M. Mattia and D. Patane, Pure. Appl. Geophys. 161, 1469 (2004). 9. S. Okubo and H. Watanabe, Geophys. Res. Lett. 16, 445 (1989). 10. D. Carbone, G. Currenti and C. Del Negro, Bull. Volc. in print (2006). 11. N. C. Lautze, A. J. L. Harris, J. E. Bailey, M. Ripepe, S. Calvari, J. Dehn, S. K. Rowland and K. Evans-Jones, J. Volcan. Geotherm. Res. 137, 231 (2004).