Integrated inversion of ground deformation and magnetic data at Etna volcano using a genetic algorithm technique

Geodetic and magnetic investigations have been playing an increasingly important role in studies on Mt. Etna eruptive processes. During ascent, magma interacts with surrounding rocks and fluids, and inevitably crustal deformation and disturbances in the local magnetic field are produced. These effects are generally interpreted separately from each other and consistency of interpretations obtained from different methods is qualitatively checked only a posteriori. In order to make the estimation of source parameters more robust we propose an integrated inversion from deformation and magnetic data that leads to the best possible understanding of the underlying geophysical process. The inversion problem was formulated following a global optimization approach based on the use of genetic algorithms. The proposed modeling inversion technique was applied on field data sets recorded during the onset of the 2002-2003 Etna flank eruption. The deformation pattern and the magnetic anomalies were consistent with a piezomagnetic effect caused by a dyke intrusion propagating along the NE direction. Mailing address: Dr. Gilda Currenti, Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Catania, Piazza Roma 2, 95125 Catania, Italy; e-mail: currenti@ct.ingv.it


Introduction
Mt. Etna is one of the most active and investigated volcanoes in the world.Over recent decades, new modern techniques of volcano monitoring have also been implemented to improve the knowledge of eruptive processes.In particular, monitoring involves geodetic and magnetic techniques that have provided essential information on the eruption mechanism including magma storage and transport within the volcano edifice.
Continuous ground deformation monitoring at Mt. Etna is performed with permanent tilt and GPS networks (Bonaccorso et al., 2002).For the different lateral eruptions that occurred between 1981 and 2001, ground deformation data analysis and modelling show early indications of intrusions initiating months to years earlier.The modelling showed that the intrusion source, i.e. the tensile mechanisms inside the volcano edifice, always involved the same structure oriented approximately NNW-SSE, which coincides with a main tectonic trend crossing the volcano (Bonaccorso and Davis, 2004).
Since the end of 1998, a permanent magnetic network equipped with Overhauser magnetometers has been operating on Mt.Etna (Del Negro et al., 2002).During the more recent lateral eruptions of Mt.Etna, significant correlations were observed between volcanic activity and changes in the local magnetic field, up to ten nanoteslas (Del Negro and Currenti, 2003).These observations are consistent with those calculated from volcanomagnetic models, in which the magnetic changes are generated by stress redistribution due to magmatic intrusions at different depth and by the thermal demagnetization at a rather shallow depth.The magnetic data not only allowed the timing of the intrusive event to be described in greater detail but also, together with other volcanological and geophysical evidence, permitted some constraints to be set on the characteristics of propagation of shallow dikes (Del Negro et al., 2004).
Geodetic and magnetic investigations have been playing an increasingly important role in studies on Mt.Etna eruptive processes.The amount of available data collected represents a valuable database, but no major effort has been made for an effective integration of this knowledge.Ground deformation and magnetic data are usually interpreted separately from each other and the consistency of the geophysical models coming from the different methods is checked only qualitatively.When the cause of their variations can be ascribed to the same volcanic source, a joint inversion of ground deformation and magnetic data would be advisable to identify the source parameters with a greater degree of accuracy.Indeed, the integrated approach involving geophysical data of different kinds ensures a more accurate solution than when single data types are considered (Nunnari et al., 2001) and leads to the best possible understanding of the physical process adding further constraints to the interpretation.
We develop an inverse technique that is robust enough with respect to the complexity of the problem and able to model different geophysical data simultaneously.The aim of inverse modelling is to find out as much as possible about the source parameters using the observations.Elaborated inverse methods combine forward models with appropriate optimization algorithms to automatically find out the best model parameters that minimizes the difference between the model values and the observations (Tiampo et al., 2000;Jousset et al., 2003).
The forward problem involves deriving a mathematical relationship that relates observations and model.Investigations on the geophysical models have revealed that they are highly non-linear and, as a consequence, the function to be minimized exhibits inherent discontinuities and multiple local minima.In such a case, techniques based on gradient methods or linearization of the optimization problem are hardly effective.That entails the need for global optimization methods.In this regard, we investigate the ability of an inversion procedure based on Genetic Algorithms (GAs) for the estimation of the volcanic source parameters from available data (magnetic and ground deformation).Among others, the GAs are considered due to their capability of performing a much broader search over the parameter space with a greater likelihood of finding the global optimal solution even in presence of local minima in the objective function (Goldberg, 1989).This modelling inversion technique was applied on field data sets recorded during the onset of the 2002 Etna eruption.In order to provide better constraints on the intrusion mechanism, we investigate joint inversion of magnetic and deformation data during the period encompassing the opening of the eruptive fissures on the NE flank in 2002-2003.

The inversion technique
Inverse problems in volcanology involve the analysis of appropriate numerical strategies to identify the source parameters of a wide spectrum of geophysical signals recorded on active volcanic areas.As considered models are highly non-linear, the inverse numerical problem could be difficult to solve by applying local optimization techniques.In particular, nonlinear models engender in the function to be minimized inherent discontinuities and local optima, where local optimization methods are likely to get stuck.On the other hand, a trial and error procedure for solving the inversion problem may take a long time and the quality of the results might not be accurate enough.Since the parameter space is typically very large, multimodal, and poorly understood, GAs have the advantage of performing a much broader search over the space parameter using a random process with a greater likelihood of finding a near optimal solution (Holland, 1975).
Roughly speaking, the adopted GA search algorithm randomly generates a set of N initial parameters (initial population), whose values are constrained within a prefixed range.Using the chosen forward model, a population of potential solutions is computed for this set of parameters.The solutions are tested by an objective function that quantifies the degree of correlation between computed and observed data.The solutions are ranked, from the best to worst, according to the objective function and a new set of likely solutions is generated from the best solutions by the laws of evolutionary genetics.The evolutionary algorithm uses a collection of heuristic rules to modify the population of trial solutions in a way that the next generation tends to be on average better than its predecessor.Heuristic selection and replacement operators are then used to sample the best solution in the current population for reproduction at each generation (Goldberg, 1989).The population evolves from one generation to the next using probabilistic transition rules defined by genetic operators (crossover and mutation).The evolution is tuned by specifying crossover (Pc) and mutation (Pm) probabilities (Michalewicz, 1992).This procedure is repeated for a large number of generations until the best solution is reached.A simplified schematic diagram of the procedure is shown in fig. 1.
The GA search strategy is determined by two crucial issues: i) exploration of the search space to investigate new and unknown areas and ii) exploitation of the knowledge found at the previous step to improve the solution.Genetic algorithms are able to find a balance between these two contradictory requirements (exploration and exploitation of the search space) by combining elements of deterministic and stochastic search.Since the randomly generated population progressively evolves by means of a more deterministic sampling of the parameter space, the GA algorithm progressively leads the population to converge efficiently and rapidly toward the most promising solution that minimizes the objective function.

The objective function
The goal of inversion methods consists in determining the best model that produces the solution closest to the observed measurements.The model can be represented by a set of parameters vector m={m 1, ... , mk}∈M, where M is the «space parameter» and k is its dimensionality.Therefore, the best model corresponds to a parameters set, which on the basis of the forward problem formulation produces ground deformation and magnetic anomalies that well fit the observed data.For a given model m, the objective function quantifies the degree of misfit between observed data and the values computed by the forward problem.In this way, it helps to define the best parameter set m. Geodetic and magnetic data are characterized by their own physical units and number of data points.As a result, the simultaneous inversion of both data sets without equalization may yield a solution that is dominated by one set.To perform the integrated inversion of ground deformation and magnetic data, we minimize the value of chi-square that accounts for the measurements error defined by the standard deviation σ, as where M k are the measured data, Ck the calculated deformation or magnetic variations and k the number of available measurements (Bevington and Robinson, 1992).

Forward modelling
The forward models used by the GA provide a mathematical formulation that relates magnetic anomalies, displacements and stress/strain fields associated to a particular volcanic source.On the whole, we have considered spherical pressure source (Mogi model in volcanology) and rectangular sources.Mogi's model is quite appropriate for modelling the inflation/deflation of the magma reservoir that can be used to describe eruptive events associated with summit eruptions.On the contrary, rectangular sources seem more realistic for interpreting eruptive phenomena in volcanic areas such as Mt.Etna, where eruptions have their origin in dykes opening from a certain depth toward the surface.
Volcanic processes are complex geophysical systems and it is difficult to derive a forward model, unless important simplifications and approximations are taken into account (McTigue, 1987).Surface displacements in a homogeneous elastic half-space have been described by Mogi (1958) for a spherical source and by Okada (1985Okada ( , 1992) ) for a rectangular fault.Also, over recent decades, the analytical formulation of the different mechanisms, which can be the cause of volcanomagnetic signals, has advanced considerably.A series of analytical solutions have been devised and widely used in the literature for mod- / elling volcano-related variations (Adler et al., 1999).The main geophysical mechanisms leading to magnetic anomalies in active volcanic areas have been investigated: i) piezomagnetic processes related to stress-induced changes in rock magnetization, ii) electrokinetic effects arising from fluid flow within the volcano edifice in the presence of an electric double layer at the solid-liquid interface, and iii) thermal demagnetization or remagnetization phenomena.
While thermomagnetic effects are mainly concerned with temperature changes within the volcano edifice, piezomagnetism and electrofiltration process are both related to stress variations.If the volcano can be assumed to be elastic, the electrokinetic and piezomagnetic fields vary with time like the strain field and the ground deformation (Zlotnicki and Le Mouel, 1988).Hence, volcanomagnetic variations are expected to accompany crustal deformation (Johnston, 1997).Therefore, an integrated inversion of the ground deformation and magnetic data could give more insights on the state of volcano than when only one geophysical parameter is considered.
As for the piezomagnetic field, Sasai (1991) succeeded in devising an analytical expression for the Mogi model, while Utsugi et al. (2000) Fig. 2. The Okada source: a fault occurring in homogeneous and isotropic elastic half-space with a uniformly magnetized upper layer.Depth is the distance between the origin O and the upper edge of the fault; strike is the orientation of the fault with respect to the north; dip is angle of the fault plane with respect the horizontal plane; length and width are the length and width of the fault respectively; slip is the vector of the displacement.calculated the solutions for strike-slip, dip-slip, and tensile-opening of a rectangular fault with an arbitrary dip angle (fig.2).As for electrokinetic effects, the analytical form of magnetic fields by an inclined vertical source in inhomogeneous media was devised by Murakami (1989).

The 2002 Etna flank eruption
In the following we interpret and model magnetic anomalies and ground deformation data that were observed simultaneously with the onset of the 2002 Etna lateral eruption.During the night of 26-27 October 2002, two fissure systems opened on the S flank and NE Rift of the volcano, feeding explosive activity and two distinct lava flows.The mainly explosive S flank eruption lasted 3 months (26 October 2002-28 January 2003), while the vents on the NE Rift erupted only for 8 days and their activity was mainly effusive.The early stages of this lateral eruption were accompanied by local magnetic field changes and ground deformation recorded at the continuously operating stations set up on the volcanic edifice (Aloisi et al., 2003;Del Negro et al., 2004).The most significant changes to occur here are associated with the opening and propagation of the vents along the NE Rift.
The dataset used for the inversion consists of the total magnetic field measurements at DGL and PDN stations, where positive variations of about 5 and 4 nT were recorded respectively between 26-27 October (fig.3).As for ground deformation data, we inverted the static offsets recorded in the tilt components (at 6 stations) and in the GPS horizontal and vertical changes (at 8 stations) during the same period (Aloisi et al., 2003).
A critical choice is the selection of the model used to interpret the geodetic and magnetic data.The magnetic variations could have resulted from stress redistribution due to dike emplacement at a rather shallow depth, which took place in a few hours.In such a case, the piezomagnetic effect is the principal mechanism accompanying the fractures opening, since it could justify both the amplitude and the time-scale of the magnetic changes.It is worth noting that the magnetic field changes might also be expected to accompany magmatic intrusions as a result of electrokinetic effects generated by fluid flow (Fitterman, 1979).However, to explain the observed rapid and reversible changes in terms of electrokinetics would require rapid and implausibly intense fluid flow (Murakami, 1989;Adler et al., 1999).Thus, even if the electrokinetic mechanism cannot be ruled out, we favor a more straightforward explanation in terms of the piezomagnetic effect.
Starting from the assumption of a piezomagnetic fault model (fig.2), we apply Okada (1985)'s elastic dislocation theory for computing the ground displacements and Utsugi (2000)'s formulation for estimating the expected piezomagnetic field.We disregard the southern dyke since it made only a minimal contribution on the recorded pattern deformation (Aloisi et al., 2003) and it was located rather far away to affect the magnetic field at DGL and PDN stations (fig.5).Hence, we invert ground deformation and magnetic data to infer the positions and the geometrical parameters of the source leading the geophysical changes.The fault model involves ten parameters: m={δ, α, D, W, L, Us, Ut, Ud, Xc, Yc}, whose descriptions and value ranges are reported in table I.The GA was requested to determine the best solution for the fault parameters in the search space.The set of likely solutions is significantly constrained by the location of the observed surface ruptures.The values of the parameters of the magnetoelastic medium used in these calculations are shown in table II.Lame's constant was set up to λ =µ that leads for Poisson's ratio of 0.25, a reasonable approximation to the values estimated in the upper crust on Mt.Etna.
The success of the search depends on the GA control parameters, explicitly i) the size of the initial population N; ii) the crossover probability Pc; and iii) the mutation probability Pm.A large number of initial individuals ensures more possibilities to explore the overall parameter space but makes the inversion procedure more time-consuming.High values for the crossover and mutation probabilities guarantee a fine exploitation and exploration, respectively, but decrease the rate of convergence toward the solution that minimizes the objective function.Initial tests suggested using a crossover probability Pc of about 0.7, which seems to be a good compromise between exploration and exploitation as well as between efficiency and computation time.In particular, we used a population of between 50 and 150, and Pm in the range of 0.3-0.7.The best results for GA setting parameters were found to be N=150, Pc=0.7 and Pm =0.6 (fig.4) for which the objective function reaches a minimum value.The reduced χ 2 value is computed by where d = N−p is the degree of freedom, N is the number of data points, p is the number of parameters.The best solution held a reduced of 3.18.We statistically compute the error for each parameter over 50 different inversions.In all of the 50 runs of the GA, the algorithm converges to the close neighborhood of the best solution starting from different random initial populations (table I).The ground deformation and the magnetic observations are reproduced well by the obtained fault model (fig.5).The computed magnetic anomalies were about 4.7 and 4.3 nT at DGL and PDN respectively.The fault location matched well the position of the eruptive fissure observed at the surface.

Conclusions
Intensive monitoring of recent lateral eruptions on Mt.Etna has produced a large amount of ground deformation and magnetic data.Although there have been numerous attempts to independently model these data, no investigations have been carried out to jointly interpret these observations.Moreover, the consistency of the models coming from these geophysical methods has been verified only a posteriori.Such a lack of multidisciplinarity and integration of observations and theoretical studies has certainly been prejudicial to the development of a deeper understanding of the volcano dynamics.
In order to make the description of volcanic sources more robust we propose an integrated inversion from deformation and magnetic data.The inversion problem was formulated following a global optimization approach based on the use of Genetic Algorithms.The GAs are considered due to their capability of performing a much broader search over the model parameters with a greater likelihood of finding the global optimal solution even in presence of local minima in the objective function.
Attempts at modelling changes in different geophysical parameters recorded in volcanic areas often involve a great deal of effort due to the complexity of the problem.The inversion problem in geophysical modelling suffers from the ambiguity and instability of its solutions.The ambiguity arises because different combinations of source parameters may lead to similar observations.Moreover, geophysical inversion problems are notoriously unstable.If the cause of ground deformation and local magnetic changes can be ascribed to the same geophysical process, a joint inversion of different data sets would be desirable in order to yield more robust estimates of source parameters and reduce the ambiguity of solution.
We applied this procedure to model the data set recorded during the 2002-2003 Etna flank eruption.The simultaneous inversion of geodetic measurements and magnetic data provides a coherent interpretation of the observations and allows delineating with more accuracy the magma intrusion process occurred at the onset in the NE flank.Even if the fault model provides a quite reasonable representation of the dyke intrusion occurred on the northern flank, it is not able to justify all the details.The simplified assumption in the forward model of elastic half-space medium disregards effects caused by the topographic variations and the presence of a multilayered medium that could play a role in the modelling procedure.To overcome this intrinsic limitation and provide more realistic models which consider topographic effects as well as complicated distribution of medium properties, more elaborated forward models are required.In this regard, numerical solutions could provide a new insight into the expected variations that cannot be achieved through analytical models, and in particular can help in distinguishing features that can be very useful for hazard evaluation.

Fig. 3 .
Fig. 3. Plot of 10-min means of differential magnetic field at PDN and DGL stations with respect to the CSR reference station between 20 October and 03 November 2002.

Fig. 4 .
Fig. 4. Testings on GA setting parameters.The population converges toward the minimum of the objective function from its initial random distribution in the parameter space.

Fig. 5 .
Fig. 5. Schematic map of the 2002-2003 Etna eruption showing lava flows erupted from south and NE-trending fissures.Locations of GPS, tilt and magnetic stations are also indicated.Recorded and computed horizontal displacements and tilts vectors are shown as well.

Table I .
Search range for the fault model parameters.The coordinates are in UTM-WGS84 projection.Parameters for the piezomagnetic fault source modelling the deformation and magnetic data recorded during the onset of 2002 NE Etna flank eruption.The standard deviation of each parameter is computed over 50 different runs of GA.

Table II .
The magneto-elastic parameters estimated in the upper crust on Mt.Etna.