1) Istituto Nazionale di Geofisica e Vulcanologia – Sezione Catania 2) GNS Science – New Zealand

1. Automated procedure for InSAR data inversion using Finite Element Method Currenti G.1, Del Negro C. 1, Scandura D. 1, Williams C.A. 2 1) Istituto Nazionale di Geofisica e Vulcanologia – Sezione Catania 2) GNS Science – New Zealand USEReST, Naples, 11-13 November 2008

2. Deformation model: Fault slip distribution Dislocation source parameters determined from multiple types of static deformation data such as GPS displacements, InSAR imagery, and surface offset measurements suggest that slip along a fault is not uniform and is best described as a distribution of dislocation sources. displacement slip Green’s function The procedure requires: 1) subdividing the faults in a finite number of patches; 2) computing the Green's Function for all the patches and the measurement points; 3) solving a linear inversion problem to determine the slip distribution.

3. Young modulus 11.5-133 GPa Poisson ratio 0.12-0.32 FEM generated Green’s function At each node of the mesh different values of Young modulus and Poisson ratio are assigned on the basis of seismic tomography investigations (Patanè et al., 2006). Medium Heterogeneity Real Topography PyLith PyLith is a finite element parallel code for the solution of dynamic and quasi-static deformation problems. www.geodynamcis.org LaGrit Mesh refinement and smoothing are available to provide a mesh with more resolution in areas of interest. www.meshing.lanl.gov

4. FEM Inversion System of linear equations n x m patches 3 x n x m unknown variables! (dip-slip, strike-slip, tensile) j ≥ 3 x n x m Non-uniquiness of the solution can arise from limited data and poor knowledge of the internal structure. This issue can be faced using: (i) as much data as possible and (ii) numerical technique to reduce the algebraic ambiguity.

5. Geodetic Data Integration EDM, GPS, leveling data assure high measurements accuracy but the coverage area is usually limited. InSAR data, despite the lower accuracy, can provide a better overview of the deformation pattern thanks to the wide coverage. A Wd covariance matrix is introduced to weight the data depending on measurements uncertainties. Linear least-squares methods for this problem require to incorporate regularization techniques in order to stabilize the problem and to reduce the set of likely solution.

6. The minimization of the quadratic functional φsubjected to bound constraints can be solved by using a Quadratic Programming algorithm based on an active set strategy (Gill et al., 1991): where: FEM Inversion & Regularization The inverse problem can be re-formulated as an optimization problem aimed at finding the unknown slip values s that minimize a data misfit and a smoothing functional: As smoothing functional, the Laplacian Operator was used to avoid large variations between neighboring dislocations.

7. φd φw L-curve criterion The L-curves of the data misfit versus the model norm as a function of the regularization parameter. The best value of regularization parameter lies on the corner of the L-curve. Input Fault Slip Distribution

8. Ground deformation data (SAR, GPS, EDM) 1 – Domain and Fault discretization: LaGrit s Procedure Scheme 2 – FEM generated Green’s function: PyLith 3 – Solving a linear inversion problem: QP algorithm

9. Analytical vs Numerical Solution We considered a 3D FEM model reproducing a rectangular dislocation source in a homogeneous and isotropic half-space and compare it with the Okada’s solution. A misfit index was used to quantify the discrepancy between the numerical and the analytical solution: Currenti et al. PEPI 2008

10. Tests on the Accuracy We compare the ground deformation achieved as the sum of the displacements generated by each patch and those obtained by the overall fault. In these computations a uniform slip distribution on the patch is assumed. The error has a RMSE (root mean square error) of the order of 10-7 m, which is well below the measurement error. The accuracy of the solution is dependent on the number of the observation points. A uniform distribution of observation points within a regular grid was assumed, centered on the ground projection of the sources. The RMSE between the assumed and the inverted slip is computed as the number of measurement points increases.

11. A case study: Northeastern flank movement at Etna volcano in 2002 On 22 September 2002, 1 month before the beginning of the flank eruption on the NE Rift, an M-3.7 earthquake struck the northeastern part of Mt. Etna, on the westernmost part of the Pernicana fault. The GPS surveys carried out in September and July 2002 shows a ground deformation pattern that affects the northeastern flank clearly shaped by the Pernicana fault. Differential interferogram for ascending scene pair 31 July 2002 to 09 October 2002. The scale indicates the phase variation along the LOS (negative values correspond to the approaching of the surface to the sensor) Bonforte et al., Bull Volc. 2007

12. Source from GPS data Inversion A numerical model was set up using the dislocation sources from GPS data inversion. Topography and medium heterogeneity are also taken into account. Bonforte et al., Bull Volc. 2007

13. 3D FEM of Mt. Etna: slip on Northeastern Flank Mesh Domain:100 x 100 x 50 km Nodes: 129253 Tetrahedra Elements: 744886 Computing Time: 15 minutes

14. A deformation pattern is clearly seen around the Pernicana volcano edifice structure. The dislocation sources were subdivided in patches to obtain the slip distribution from the inversion procedure. Sar Data and Source Discretization Line of sight change map calculated by the unwrapped interferograms, processed using ROI_PAC software, developed by Jet Propulsion Laboratory & Caltech; Courtesy of F. Guglielmino ed A. Spata.

15. FEM simulations 87 1 2 744 Green Functions Computations 248 fault patches Dip-slip, strike-slip and tensile dislocation were assigned to the nodes lying along the fault surface. For each patches, in which the fault is discretized, Green’s functions were computed using PyLith. For each each simulations, the accuracy was warranted checking the convergency of the FEM solution. The iteration of GMRES solver is stopped when a threshold of 10-9 is reached or the number of iterations is higher than 200. 248 x 3 = 744 FEM simulations By a linear speedup on a cluster of 20 nodes the computing time reduces from 10 days to 10 hours.

16. Inversion Results SAR data The kinematic of the faults was constrained to those derived from GPS data inversion. The solution in correspondence of the corner of the L-curve provides a deformation pattern which resembles the SAR data and also matches the GPS obtained from the analytical solution. Inverted Solution

17. Slip Distribution Strike-slip, dip-slip and tensile dislocation distribution.

18. Model Comparison The GPS data inversion provides a deformation pattern which misses the clear anomaly around the Pernicana fault system. A heterogeneous distribution of the slip along the structure is able to better justify the SAR data. The FEM model based on the numerical inversion provides a more complex pattern which is likely to be expected. SAR data Numerical Model Analytical Model

19. Conclusions An automated procedure for geodetic data inversion is proposed to estimate slip distribution along fault interfaces. 3D finite element models (FEMs) are implemented to compute synthetic Green’s functions for static displacement. FEM-generated Green’s functions computed using PyLith, a source-free parallel finite element code, are combined with a Quadratic Programming algorithm to invert ground deformation data and obtain an estimate of slip distributions along seismogenic structures. The procedure was applied to study the ground deformation preceding the 2002-03 Etna eruption. SAR images showed a significant deformation pattern of the north-eastern flank of the volcano involving the main local volcanic edifice features on the northeastern flank. The numerical model highlighted a heterogeneity slip distribution along the Pernicana fault with a predominant strike-slip mechanism associated with a dip-slip movement in the western part. The time consuming computation of the procedure is related to the generation of the Green’s functions. For the main volcano edifice features and the seismogenic structures they has to be computed only once and stored in a database.

20. Numerical Solution Several computations were performed to better understand the effects induced by the topography and the medium heterogeneity. The topography significantly alters the general pattern of the ground deformation especially near the volcano summit. On the contrary, the medium heterogeneity does not strongly affect the expected deformation. Currenti et al. PEPI 2008

21. Smoothing Functional The Laplacian Operator is introduced to avoid large variations between neighboring dislocations The matrix WL is constructed so that the n-th row in WL contains coefficients in the equation above for columns corresponding to the neighboring source.