FEM and ANN combined approach for predicting pressure source parameters at Etna volcano

A hybrid approach for forward and inverse geophysical modeling, based on Artificial Neural Networks (ANN) and Finite Element Method (FEM), is proposed in order to properly identify the parameters of volcanic pressure sources from geophysical observations at ground surface. The neural network is trained and tested with a set of patterns obtained by the solutions of numerical models based on FEM. The geophysical changes caused by magmatic pressure sources were computed developing a 3-D FEM model with the aim to include the effects of topography and medium heterogeneities at Etna volcano. ANNs are used to interpolate the complex non linear relation between geophysical observations and source parameters both for forward and inverse modeling. The results show that the combination of neural networks and FEM is a powerful tool for a straightforward and accurate estimation of source parameters in volcanic regions.


Introduction
Inverse modeling in volcano geophysics deals with the estimate of the parameters of a source that causes significant changes in geophysical observations recorded by monitoring networks.Inverse problems are usually formulated and solved as optimization problems based on iterative procedures, minimizing an objective function that quantifies the misfit between the observed data and the estimated solutions from forward models.Analytical solutions are often used to represent the forward model because of computational convenience and fast computer implementation (e.g., Currenti et al., 2005;Nunnari et al., 2005).The intrinsic limitation of inverse methods based on analytical solutions is that irregular Correspondence to: G. Currenti (currenti@ct.ingv.it)geometries, complicated distribution of medium properties, and topography effects cannot be accounted for.Neglecting these features that strongly affect the solutions, an inaccurate estimate of source parameters is obtained (Currenti et al., 2009;Williams andWadge, 1998, 2000;Trasatti et al., 2008).To overcome this intrinsic limitation and provide more realistic models, numerical solutions based on Finite Element Method (FEM) have been investigated (Cayol and Cornet, 1998;Williams and Wadge, 2000;Currenti et al., 2007Currenti et al., , 2008Currenti et al., , 2009;;Trasatti et al., 2003;Lungarini et al., 2005).However, the use of numerical forward models in iterative methods is computationally expensive since the estimate of the objective function requires to perform a full FEM analysis at every iteration step.For this reason global optimizations based on numerical solutions are rarely investigated (Fukushima et al., 2005;Trasatti et al., 2008).As traditional optimization algorithms cannot "learn", they cannot benefit from solutions obtained previously for similar problems and each new inversion requires the minimization procedure to be re-iterated.
Recently, Artificial Neural Networks (ANNs) have been introduced to solve the inverse problem in many research applications (Haykin, 1999;Arena et al., 1998).The main advantage of inverting with ANNs consists in the availability of an approximation of the inverse model, avoiding a search for the minimum and speeding up the computation of the optimal solution that fits the observed data.ANNs have been widely used to invert geophysical models based on straightforward analytical solutions (Langer et al., 1996;Maugeri et al., 1996Maugeri et al., , 1997;;Nunnari et al., 2001) because of low computational effort.On the contrary, FEM-based numerical solutions are not often combined to ANN inversion scheme because they are time-consuming both in length of time required to design a mesh and in actual computation time.With the advent of today's powerful computer resources and the automation of mesh generation and FEM analysis, hybrid schemes based on FEM and ANN have been proposed in Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
We investigate the ability of a hybrid procedure in which ANNs are used for system identification of forward and inverse geophysical models solved by FEM.The procedure is applied to model the expected geophysical changes at Etna volcano, which constitutes a unique natural laboratory for the understanding of eruptive processes because of the high rate of eruptive events and the high quality geophysical observations gathered by monitoring networks.The 3-D numerical model allows to include not only the real topography of Mt Etna but also heterogeneous medium properties inferred from geological evidences and seismic tomography investigations (Chiarabba et al., 2000;Tibaldi and Groppelli, 2002).The ANNs are trained with numerical patterns, obtained computing geophysical changes caused by pressure sources through FEM analysis.We propose an integrated approach in which geophysical data of different nature, that can be ascribed to the same volcanic source, are jointly modeled in order to identify the source parameters with a greater degree of confidence than when only one kind of data is used (Nunnari et al., 2001).
The application of FEM and ANN based modeling is presented to solve both forward and inverse problems.Synthetic patterns are generated by FEM in order to provide a data set large enough to represent the training and testing sets of the possible models within the model space.Firstly, we estimate the volcanic source parameters performing the inversion of ground deformation, magnetic and gravity fields, which are continuously recorded by permanent monitoring stations at Mt Etna (Bonaccorso and Davis, 1999;Bonforte et al., 2008;Del Negro et al., 2004;Napoli et al., 2008;Carbone et al., 2007Carbone et al., , 2008)).We use the ANN to identify the inverse relation between the geophysical observations and the source parameters.The ANNs, once properly trained, can solve the inversion problem very fast and with an appreciable degree of accuracy.Secondly, we used the ANN to map the forward numerical model with the aim to obtain an immediate approximate solution to be used in any inversion algorithm.This approach retains the computational convenience of forward analytical solutions and is also capable of including the effects of topography and medium heterogeneity.Our results show that the association of FEM and ANN techniques could be a useful alternative for advancing the model-based assessments of geophysical observations in volcanic areas.

Forward problem: FEM solution
Because of the limited amount of volcano-related geophysical changes detected from measurements, we generate synthetic patterns to provide a large number of input/output patterns to train and test the ANN.We compute ground defor-mation, magnetic and gravity field changes at Mt Etna caused by spherical pressure sources that are quite appropriate for modeling inflation/deflation of magma reservoirs, which are likely to reside at depth greater than 2 km b.s.l.(Corsaro and Pompilio, 2004;Bonforte et al., 2008).
We develop a numerical procedure based on FEM to evaluate deformation field changes.The elastostatic problem is solved in a computational domain of a 100×100×50 km using the software PyLith (Aagaard et al., 2008).The horizontal and vertical displacements at the lateral and bottom boundaries of the domain are fixed to zero, representing the vanish displacement at the infinity.The upper boundary is stress free and represents the ground surface.These settings warrant a good accuracy of the numerical solution (Currenti et al., 2008).
Both irregular topography and medium heterogeneities are included in the FEM model to consider a realistic description of Etna volcanic edifice.The real topography of Mt Etna was generated using a digital elevation model from the 90 m Shuttle Radar Topography Mission (SRTM) data and a bathymetry model from the GEBCO database (http: //www.gebco.net/).The elastic material properties are derived from seismic tomography investigation (Patané et al., 2006).We use P-wave and S-wave seismic velocities to define the Young modulus using the following equation (Kearey and Brooks, 1991): where V p is the seismic P-wave propagation velocity and ρ is the density of the medium, which was set at an average value of 2500 kg/m 3 (Corsaro and Pompilio, 2003).The values of Poisson's ratio are obtained using the equation (Kearey and Brooks, 1991): where V s is the seismic S-wave propagation velocity.On the basis of Eqs. ( 1) and ( 2), within the computational domain the Young modulus varies from 11.5 GPa to 133 GPa, while the Poisson ratio is in the range 0.12-0.32.Pressurization of magmatic sources is also accompanied by gravity and magnetic changes, generally ascribed to the variations of density and magnetization within the magmatic source.The thermo-magnetic effect, caused by thermal demagnetization of magnetic materials, is a mechanism able to produce large magnetic changes.Other magnetic effects, such as piezomagnetism (Sasai, 1991), are generally not larger than few nanoTeslas (nT).As observed during the last Etna eruptive events (Del Negro and Currenti, 2003;Del Negro et al., 2004;Napoli et al., 2008), detectable piezomagnetic changes are induced by shallow magmatic intrusions within the volcanic edifice (Currenti et al., 2009), whereas negligible piezomagnetic variations are generated by deep pressure sources (>2 km b.s.l.).
Gravity changes due to a pressure source originate from three different contributions: (i) new mass input from remote distances into the source volume, (ii) relative volume change of the medium arising from compressibility of the surrounding rock, (iii) the displacements of density boundaries (Bonafede and Mazzanti, 1998;Okubo, 1991;Currenti et al., 2007).The last two contributions have a significant effect on the magnitude of the predicted gravity changes when the source resides at shallow depths (Bonafede and Mazzanti, 1998;Charco et al., 2006;Currenti et al., 2007).In case of deep pressure sources these contributions are below microgravity accuracy (Battaglia and Segall, 2004;Currenti et al., 2007) and are negligible compared to those produced by the input of new mass.
Under these assumptions, we focus on thermomagnetic effect and gravity changes related to mass input, which depend mainly on the magnetization and density contrasts between the source and the surrounding rocks and on the source position.Therefore, they are not affected by the presence of medium heterogeneities.Moreover, for spherical sources the topography effect is due primarily to the distance of the free surface from the magma source rather than the local shape of the free surface (Williams and Wadge, 1998).As a result, no differences are found between numerical and analytical solutions, if the topography effect is approximated by replacing in the analytical expressions the flat reference surface with the real elevation of the observation point (Currenti et al., 2007;Charco et al., 2009).At the observation point (x, y, z) the analytical solution of the thermomagnetic change T due to a source centered at (x 0 , y 0 , z 0 ) is expressed as follows (Blakely, 1995): where M m is the magnetic moment, R is the radius of the sphere, m is the magnetization, I is the magnetic inclination, and r = (x 2 is the distance between the source center and the observation point on the real topography.The gravity change g due to input of mass M g in a spherical source is expressed as follows (Blakely, 1995): where ρ is the density contrast, and G is the gravitational constant.-9 -1 Using the forward models, we generate 1050 synthetic patterns of deformation, magnetic and gravity changes for training and testing the neural network.The procedure of pattern generation is divided in several steps and is executed automatically.Firstly, the computational domain of Mt Etna is meshed into 598 948 isoparametric and arbitrarily distorted tetrahedral elements connected by 103 424 nodes.The meshed domain, generated with the LaGrit software (http://lagrit.lanl.gov),has a spatial resolution of 300 m in the summit area and around the source location (Fig. 1).Then, the parameters of the source are generated with a random distribution in the ranges reported in Table 1.The sources are distributed in a volume that contains all the pressure sources active at Mt Etna during the last decades (Bonforte et al., 2008).
Once all the source parameters are defined, the source volumes are meshed iteratively and introduced in the domain to finally obtain 1050 complete meshes, characterized by the different positions of the source.For each mesh, PyLith was automatically run to solve the elastostatic problem.Given the range of the source depth, the effect of the finite shape of the pressure source can be disregarded (Williams andWadge, 1998, 2000;Bonafede and Mazzanti, 1998;Tiampo et al., 2000;McTigue, 1987), and the radius R is kept constant to 500 m.As the numerical solutions of the deformation is proportional to the pressure change, all the computations are performed using a constant pressure of 100 MPa and the solutions are then rescaled within the range of the random value of P (Table 1).The accuracy of the numerical solution is warranted checking the convergence of the GMRES solver (Aagaard et al., 2008), stopping the iterations when the relative tolerance reaches a threshold of 10 −9 .The computing time of a FEM analysis is about 10 min.By a linear speedup on a cluster of 20 nodes the computing time for the generation of the 1050 patterns reduces from 8 days to 9 h.Finally, the solutions are interpolated at the coordinates of the continuously running GPS stations (Fig. 2).For the same source parameters, magnetic and gravity changes at the locations of the monitoring stations are computed using Eqs.( 3) and ( 4).The gravity and magnetic changes were calculated assuming values of magnetization and density contrast between 4 and 8 A/m (Currenti et al., 2009) and between 100 and 300 kg/m 3 (Corsaro and Pompilio, 2004), respectively.

Neural network model
Artificial Neural Networks (ANN) consist in a large number of simple processing elements, called neurons that are connected to each other by means of directed links, each with an associated weight (Haykin, 1999).An activation function controls the amplitude of the neuron's output, adjusting a linear combination of the weighted inputs.ANNs are typically applied for approximating a non linear input-output relationship of the form y=f (x), where x and y are the in- put and output vectors, respectively.Using the supervised learning scheme (Zurada, 1999;Haykin, 1999), the inputoutput relationship is determined by looking at the examples of many input-output pairs, which represent the knowledge of the model.In the learning process, the connection weights of the network are changed adaptively to match the actual output of the ANN with the example output.The procedure iterates until the error is small enough.
We apply a class of neural networks called Multilayer Perceptrons (MLPs) to identify both forward and inverse geophysical models.The neural network is trained using many input-output pairs generated solving the forward model.After the training process, the ANN provides a model M that approximates the geophysical model.
In the inverse problem, the model M approximates the function x=f −1 (y), in the hypothesis that it exists.The ANN inverse model M i provides an estimate of the source parameters (x) as a function of the geophysical observations (y).Once the neural network is trained, it is able to approximate the unknown function and to identify, for a set of geophysical observations, the source parameters that better reproduce these variations (Fig. 3a).
In the forward problem, the model M approximates the function f (•) providing an estimate of the geophysical observations (y) from the assigned source parameters (x).Closedform analytical solutions providing the function f (•) are available under the assumption of homogeneous half-space medium, whereas the function f (•) is unknown for complex numerical models, which consider medium heterogeneities and real topography.Therefore, the ANN is trained to identify the forward numerical model M f using input/output patterns generated from FEM simulations that provide more realistic models.The block diagram of the ANN forward model is reported in Fig. 3b. 4 Identification results

ANN based inverse model
We firstly implement the ANN inverse model M i of deformation field, and then we investigate a model M i able to invert deformation, magnetic and gravity changes together.
The MLP neural network is initially trained with a set of 800 input/output patterns, successively the inverse model is tested with a set of 250 patterns that has not been used previously for training.Several network configurations for the deformation inverse model have been tried and the better results are obtained using a three layered network with sigmoidal activation functions.The inputs are the three components of deformation field calculated at 15 stations of the continuously running GPS network and the outputs are the source parameters (position of the source and pressure change).We choose a network with a single hidden layer because we have found experimentally that more than one hidden layer leads to a considerable increase in time required for the learning process without any significant improvement in the accuracy of the solution.Ten hidden neurons are verified to be enough to obtain a compromise between performance and complexity of the network.The training of the ANN takes about 6 h.The accuracy of the inverse model is tested using the Root Mean Square Error (RMSE), whose expression is: where xi is the calculated value of the i-th parameter, x i is the corresponding true parameter, and N is the number of pattern.Since the source parameters have different scales, we also computed the mean percent error E% abs (Nunnari et al., 2001) to have a normalized measure of the misfit: where D is the range of the parameter.These performance indexes related to the inversion of deformation data are reported in Table 2.The ANN is able to approximate with good accuracy the inverse model.For the estimate of the source position the RMSE is less than 200 m and the E% abs does not reach 1%.A worse accuracy is obtained for the pressure change whose E% abs is about 3%.
To assess the accuracy of numerical ANN inverse model M i , we also train and test the network with a set of inputoutput patterns calculated analytically for a homogenous half-space model (Appendix A).Better performances are achieved for the analytical inverse model with respect to the numerical one.The E% abs for all the source parameters is lower than 0.15%.The main discrepancy between the two models is observed on the pressure change.All source parameters have quite the same value of E% abs in the analytical inversion, whereas in the numerical inverse model the error on the pressure change is higher with respect to the other parameters.This result could be ascribed to the complex distribution of elastic parameters that affects the deformation of the rock surrounding the source.Under the assumption of elastic rheology, the amplitude of the deformation field for a spheroidal pressure source is linearly related to the ratio between the pressure and the rigidity modulus ( P /µ).Therefore, it is not straightforward to distinguish whether the perturbation in ground deformation is due to changes in overpressure or in rigidity modulus.The spatial variation of elastic parameters in the 3-D numerical model makes the estimate of pressure changes more difficult.However, in the numerical model the error on the estimate of pressure change is still acceptable (Table 2).Even if the inverse model to interpolate is more complicated because of the perturbations due to medium heterogeneity and irregular topography, the numerical inversion can be achieved with good accuracy.
We train the neural network with an integrated dataset of ground deformation, magnetic and gravity fields.The neural network used for the integrated inversion is a three layered network, with 59 inputs (the three components of deformation field calculated at 15 stations, thermomagnetic www.nonlin-processes-geophys.net/17/273/2010/ Nonlin.Processes Geophys., 17, 273-282, 2010   testing set, which has not been used for the training.A good match is observed for all source parameters since the real and predicted values lie along the bisector axis with small dispersions (Fig. 4).These results demonstrate that the estimated values are well correlated with the ANN predictions.
The RMSE and E% abs for the integrated inversion (Table 2) are similar to those obtained inverting only deformation data.
The integrated approach, relying on observations of different kinds, allows to estimate the variations in different geophysical parameters, which independently give insights into the volcano activity.The large number of geodetic measurements usefully supports the inversion of magnetic and gravity observations, coming from a limited number of monitoring stations.
To demonstrate the robustness of the numerical inversion procedure, both for deformation inversion and for integrated inversion, tests were carried out adding white noise to the set of inputs.The value of the noise was calculated as follows (Nunnari et al., 2001): where Xran is a uniformly distributed random variable in the interval [-1; 1], Perc is the percentage of noise, y i is the i-th value of the considered input set (Fig. 3a).The ANN trained with noise-free patterns is tested using noise-added patterns to evaluate the robustness of the network for real observed data, which are intrinsically affected by measurement error.We evaluated the E% abs performance index for different noise levels increasing the Perc values from 1 to 10.The mean percent error E% abs is showed in Fig. 5 for numerical inversion of ground deformation, and in Fig. 6 for integrated numerical inversion, for both noise-free and noiseadded conditions.Although the E% abs performance index

ANN based forward model
The proposed ANN inverse model has the intrinsic limitation to be strictly dependent on both the configuration and the number of measurement stations.If a station does not work and no data are available for one of the input nodes, the learning process has to be performed again using a different neural network structure compliant to the new stations' configuration.To overcome this weakness, a different inversion scheme should be investigated.Instead of identifying the inverse model, the ANN can be trained to provide an approximation of the numerical solutions.Once the network is trained, the numerical solutions of the forward model are immediately available as outputs of the ANN with a computational convenience comparable to the analytical solutions (Fig. 3b).The inversion can be performed at a second time also with iterative methods since the solution of the forward model is straightforwardly computed by the ANN, which avoids the computation effort of the full FEM analysis at every iteration step.Moreover, the inverse process is more compliant and flexible to the variation of stations' configuration.
The ANN is trained to identify the forward numerical model M f only for ground deformation field, since the forward models for the gravity and thermo-magnetic fields can be approximated using closed-form analytical solutions .The ANN identifies the unknown relation between the source parameters and the deformation field in presence of medium heterogeneity and topography (Fig. 3b).The neural network is trained to interpolate the function u i = F (Xc, Yc, Zc, rx i , ry i ) at each station, where u i is the deformation at i-th station, Xc, Yc, Zc give the position of the source, rx i , ry i are the horizontal components of the relative distance between the source and the station.The neural network used for identifying the forward model is a three layered network, with 5 inputs (Xc, Yc, Zc, rx i , ry i ), and 3 outputs (u x , u y , u z ).The hidden layer is composed of 20 neurons.Sigmoidal activation functions were chosen for the first two layers, while a linear activation function was chosen for the output layer.The training of the ANN for every station requires about 10 min.
Using the input/output patterns generated with the FEM, the ANN forward model is trained for each GPS station.The model is validated computing the RMSE for the three deformation components.To provide a reference for the quality of the identification, we rescale the misfit values by the magnitude of the deformation field defining a normalized measure of the error δ m : where u m is the m-component of the computed deformation using FEM and ûm is the deformation estimated by the ANN, with m=x, y, z computed for each deformation component.
In Table 3, we show the results for four stations, two situated in the summit of volcano and the others at a lower altitude.
The results demonstrate that the neural network is able to identify with a good accuracy the forward numerical model of deformation field despite the complexity introduced by the medium heterogeneity and topography.It is worth noting that the errors at the stations located out from the summit are smaller than those obtained for the summit stations, where numerical solutions are strongly influenced by topography and heterogeneity (Fig. 2).This result highlights the significant effect that these parameters have on the solution.

Conclusions
A combined approach is proposed using neural networks and numerical models for the estimation of pressure source parameters in volcanic regions.ANN is trained to identify the unknown model using patterns obtained by numerical FEM solutions.This method permits to interpret geophysical data avoiding the intrinsic limitation of analytical solutions and providing a more realistic description of volcanic processes.
Both ANN and FEM have demonstrated to be appropriate for representing geophysical models at Mt Etna, taking the ability of ANN to learn complex non-linear behavior and that of FEM to include 3-D realistic features in models.Firstly, ANN is used in inverse scheme to identify the source parameters from geophysical observations.The results show that, notwithstanding the high nonlinearity of the considered inverse problems, it can be solved with satisfactory accuracy both for ground deformation inversion and for the integrated inversion of ground deformation, magnetic and gravity changes.Moreover, the inverse model provides an acceptable estimate of source parameters even in presence of noise, which intrinsically affects the geophysical measurements.
Secondly, neural approach is used to map the forward numerical model, reducing the high computational time usually necessary to run the FEM model in iterative inversion schemes.This approach permits to have a straightforward solution of numerical model that can substitute the analytical solution in any inversion technique, providing a more realistic description of geophysical changes in volcanic area by taking into account the real topography and medium heterogeneities.
The a priori assumption of the source shape is a limitation of our procedure.However, because of its flexibility, the procedure can be easily extended to consider more complex geometries with higher number of parameters in a way to provide a wider class of models.This would allow to make a comparison among several possible solutions and to select the model geometry that better reproduces the geophysical signals.
Results show that the neural network/finite element modeling can be efficiently used in the identification of both forward and inverse numerical models.This approach allows for an accurate estimation of source parameters in a realistic description of volcanic areas, which is of great significance for both the correct interpretation of geophysical data and for the assessment of volcanic hazard.

Appendix A
The analytical solutions of the three components of ground deformation field are computed using the simple and common Mogi source embedded in a homogeneous half-space medium (Mogi, 1958).The components of the displacements u x , u y , u z , computed at the observation point (x, y, z), with the vertical z-axis directing downwards, are the following: where λ and µ are the Lamé constants, d is the depth of the source, C = P R 3 2 , R is the radius of the sphere, P is the pressure change at the source wall, r 1 = x 2 + y 2 + (z − d) 2 , r 2 = x 2 + y 2 + (z + d) 2 .

Fig. 1 .
Fig. 1.Mesh of the computational domain.The mesh has a spatial resolution of 300 m in the summit area and around the source location and becomes coarser at greater distance.

Fig. 2 .
Fig. 2. Permanent GPS (red circles), gravity (yellow squares) and magnetic (blue triangles) stations of the monitoring networks on Mt Etna.The black rectangle corresponds to the projection on surface of the volume where the sources are located.

Fig. 4 .
Fig. 4. Predicted values of source parameters with respect to output patterns of the testing set for the integrated numerical inversion.

Fig. 5 .
Fig. 5. Mean percent error E% abs for each source parameter using noisy deformation patterns at the inputs of the ANN inverse model.

Fig. 6 .
Fig. 6.Mean percent error E% abs for each source parameter using noisy patterns at the inputs of the integrated ANN inverse model.

Table 1 .
Ranges of the random generated parameters of the source.

Table 2 .
Performance indexes RMSE and E% abs for the inversion of analytical and numerical deformation model and for the integrated numerical model.

Table 3 .
Performance indexes RMSE and normalized misfits δ m for four stations located near (A and B in Fig.2) and far (C and D in Fig.2) from the summit of volcano.