FEMSA: a finite element simulation tool for quasi-static seismic deformation modeling

We set up a computational tool to numerically model static and quasi-static deformation generated by faulting sources embedded in plane or spherical domains. We use a Finite Element (FE) approach to automatically implement arbitrary faulting sources and calculate displacement and stress fields induced by slip on the fault. The package makes use of the capabilities of CalculiX, a non commercial FE software designed to solve field problems (see  for details), and is freely distributed by request. Mailing address: Dr. Manuela Volpe, Istituto Nazionale di Geofisica e Vulcanologia, Sezione Roma1, Via di Vigna Murata 605, 00143 Roma, Italy; e-mail: volpe@ingv.it


Introduction
Earthquakes are generated by abrupt motions along a fault plane.The Earth, behaving like a (nearly) elastic body, instantaneously reacts to the fault dislocation and generates elastic waves as a result of the transient stress imbalance in the rock.Since the dislocation is permanent, seismic waves decay to a residual static offset called coseismic deformation.Modeling stress and displacement fields resulting from seismic faulting is an increasingly valuable tool in earthquake seismology.The progress in geodetic techniques such as GPS and SAR has made available high precision coand post-seismic displacement measurements, which have been employed, for instance, in the inversion of the seismic source structure (Árnadóttir et al., 2005;Banerjee et al., 2005;Savage et al., 2005).
In the last two decades, many different approaches have been proposed to study the static and quasi-static deformation associated with a seismic event.On the one hand, although analytical and semi-analytical models (Okada, 1985(Okada, , 1992;;Rybicki, 1986;Roth, 1990Roth, , 1993;;Hudnut et al., 1996;Sagiya and Thatcher, 1999) are often very accurate and simple to implement, they suffer from a limited flexibility, as they give very little opportunity to account for geometrical and/or rheological complexities.Starting from the 1990's, most studies focused on deriving analytical formulations of elastic dislocations in a layered half-space or spherical domain, possibly including rheological layering (Ma andKusznir, 1992, 1994;Pollitz, 1992;Piersanti et al., 1995Piersanti et al., , 1997;;Savage, 1998;Bonafede et al., 2002;Rivalta et al., 2002).In some cases, formulae were derived for certain types of faulting or for certain types of displacements only.Some of the models are difficult to use and have only a theoretical interest.Moreover, the effects of crustal rigidity layering and the lateral variations of rheological properties are not accounted in all of the semi-analytical analyses associated with faulting, because they lead to mathematical problems which cannot be solved analytically.On the other hand, purely numerical models have been developed to overcome restrictions and improve performances with respect to analytical techniques, but in some cases they also result unable to model all the significant features or geometric complexities of interest e.g., Megna et al. (2005).Among numerical methods, however, Finite Element (FE) modeling in three dimensions allows accurate modeling of complex faulting geometry, inhomogeneous materials and realistic viscous flow, appearing an excellent tool to investigate many specific phenomena related with earthquakes.Unfortunately, at present FE analyses are mostly based on commercial or in-house codes, which in fact in-volve other practical disadvantages since they are both hardly accessible to the whole scientific community for different reasons.Commercial codes have often prohibitive prices and are not very easy to handle because of the difficulty of going into their real working, as well as for the actual lack of any user-code interaction.Moreover, most commercial codes are aimed at industrial design problems and their adaptation to geophysical problems is not straightforward (Wu, 2004).In-house codes, in several cases, turn out to be easily understandable and usable just by their developers and are often designed for specific problems, resulting hardly extendable to more general cases.Furthermore, they are not particularly robust with respect to errors or inaccuracies, as they are often based on ad The aim of the present work is to set up and test a flexible, versatile, robust and reliable numerical method to investigate crustal deformation produced by complex rupture along dislocations.The main goal is to achieve a simulation tool which is easy to implement and handle and especially based on widely used (i.e.not only in geophysics) but non commercial packages.To this purpose, we make use of a free three-dimensional FE software, CalculiX (v.1.5), distributed under the terms of the GNU General Public License (see <http://www.calculix.de>for downloading).The CalculiX package is widely used among the scientific community for a great variety of applications and provides both the mesher and the solver to set up and solve the FE model.Auxiliary codes have been developed to work as appropriate interfaces to the software.The tool has been called FEMSA, that is the acronym for «Finite Element Modeling for Seismic Applications», and it is freely distributed by request (email to volpe@ingv.it).A user manual with tutorial enclosed is also available.
The paper is organized as follows: in Section 2 we present the capabilities of CalculiX and its key features; in Section 3 we discuss our approach and give some additional theoretical and computational details concerning the FE model; in Section 4 we exhibit a series of benchmark simulations performed on an homogeneous test domain (see fig. 1a,b) in order to validate the model and check the functionality and flexibility of our approach; finally in Section 5, just to illustrate some potentialities of FEMSA, we investigate few preliminary applications on a simple rheologically heterogeneous model, representing a border zone between the oceanic and the continental crust, and on a real earthquake fault, namely the 2004 Sumatra earthquake.Concluding remarks are summarized in Section 6.

The CalculiX package
CalculiX is a package designed to solve field problems by means of the FE method.The copyright (1998) belongs to Guido Dhondt and Klaus Wittig (<http://www.calculix.de>).A useful ex-change forum to share experiences, bring forward problems and discuss improvements is available at the web site <http://groups.yahoo.com/group/calculix>.In the following, the main features of CalculiX will be briefly described.For more details the reader is referred to the specific documentation.
With CalculiX, FE models can be built, calculated and post-processed.Naming conventions and input style formats for CalculiX are based on those used by ABAQUS, a proprietary, general purpose finite element code developed and supported by Hibbitt, Karlsson & Sorensen, Inc. (HKS) (see <http://www.abaqus.com>for details).Table I summarizes some of the most commonly used element types.
The pre-and post-processor (CalculiX GraphiX: cgx) is an interactive 3D graphical interface using the OpenGL library for visualization and the GLUT library for window management and event handling.The cgx creates a structured mesh based on a description of its geometry and is able to generate and display beam, shell and brick elements in linear and quadratic form, while it can display, but not create, penta-and tetrahedral-elements.In addition, it supports common CAD formats.After the mesh is created it must be written to a file to be available for the solver.Results can be visualized by calling the cgx program again in an independent session.The post-processor functionality is mainly controlled by a pop-up menu, while the pre-processor functionality is controlled by the keyboard.
The solver (CalculiX CrunchiX: ccx) is able to perform linear and non-linear calculations on beam, shell, brick, tetrahedral and wedge elements.Each element in the structure must have a material assigned, with specific linear or nonlinear mechanical properties.Static, dynamic and thermal solutions are available.A ccx input deck basically consists of a model definition section describing the geometry and boundary conditions of the problem and one or more steps defining the forcing terms.As far as the theory behind ccx is concerned, the reader is referred to Dhondt (2004) and Cook et al. (2002).
Both programs can be used independently.Since the solver makes use of the Abaqus input format, it is possible to use other pre-processors as well.In turn, the pre-processor is able to write mesh related data for Nastran, Abaqus, Ansys, Code-Aster and for the free-cfd codes duns, ISAAC and OpenFOAM.The package can run on Unix and MS-Windows platforms.

The simulation method
Our approach basically consists in exploiting CalculiX capabilities and building interface codes to investigate postseismic deformation produced by dislocations within plane or spherical domains by means of a three dimensional FE modeling.In other words, the goal is to develop a numerical tool to automatically implement arbitrary faulting sources and perform static and quasi-static calculations to compute the horizontal and vertical displacement and stress fields induced by the fault slip.A reference analytical solution obtained with the Okada model (Okada, 1985(Okada, , 1992) can be automatically provided for comparison.
In FE modeling, faults are commonly considered like contact interfaces between deformable bodies with stick and finite frictional slip (Xing andMakinouchi, 2000, 2002;Cianetti et al., 2005), i.e. by a pair of surfaces whose nodes are connected by contact elements preventing copenetration of the deformable surfaces.During fault slip, these elements are able to account for frictional properties if a suitable friction law relating tractions to slip, slip rate or other fault state variables is given.
So far, CalculiX contact capabilities are not suitable for our aims, but such a restriction can be overcome by introducing an alternate approach based on the equivalent body force theorem.In fact, it is well known that seismic sources can be mathematically described by means of dynamically equivalent force systems replacing the actual process (Burridge and Knopoff, 1964).The displacement field produced by a fault dislocation is equivalent to the one associated with a distribution of double couples of forces placed in an unfracturated medium.Since coseismic deformations decay rapidly with distance from the source, most modeling involves vertical and horizontal ground displacements observed near the fault.Consequently, a point source approximation is not valid and a finite fault with a numerically discretized distribution of double couples is needed (Piersanti et al., 1997).This approach binds us to handle planar or quasi-planar fault surfaces, but such a limit is likely to be overcome as far as frictional contact elements will be implemented in CalculiX.
The fault plane, represented by an appropriate planar distribution of double couples simulating either a strike or a dip slip, is generated by means of a Fortran interface, which also attends to define the force field.Basically, nodes in groups of four, corresponding to force application points, are suitably picked according to the slip vector and, if needed, moved in order to obtain the correct orientation depending on the dip angle (δ).For the sake of symmetry and computational convenience, strike angles (φ) different from zero are addressed by rotating the reference frame of the displacement field produced by a fault with zero strike.As a final note, no slip direction other than pure strikeslip or pure dip-slip mechanisms is allowed.In any case an arbitrary rake angle λ ≠ 0°, 180°, 90°, 270°is easily attainable by means of a simple superposition of a pure strike and a pure dip, each having seismic moment M s=M0 ⋅ cos λ and Md=M0⋅ sin λ , respectively.
The package is built up to operate automatically.Once the fault plane with zero strike and arbitrary dip angle is generated, the same Fortran code creates an executable script to perform the following operations: i) the displacement field is analytically calculated according to the Okada model; ii) suitable boundary conditions are formulated and formalized as explained in the following; iii) the FE simulation, that is the most time consuming step, is carried out; iv) reference frame rotation is applied to the numerical solution in order to assign an arbitrary strike angle.For plane domains, the analytical solution for the system having arbitrary strike and dip is also computed for comparison.Obviously, if zero strike is expected, the procedure stops after the FE simulation.On the other hand, different strike angles can be accounted for just by point iv), performed by a suitable script generated by another Fortran interface.Such a strategy remarkably reduces the computational demand, since only one FE simulation is needed for each dip angle.This is highly desirable especially in case of inverse applica-tions.In fig.2, a block diagram of the described procedure is shown.

Benchmark simulations
In the present work, the FE test domain is represented by a 400 ×400×100 km hexahedral volume discretized into 13539 twenty-node brick elements (i.e.quadratic elements), resulting in a mesh containing 60288 nodes.Such a mesh takes about 20 min to perform a static calculation on a 64-bit Intel Xeon 3.2 GHz, 8 GB RAM.The mesh is built up to be finer within a central cubic region spanning 100×100×100 km, while outside of this region the spatial resolution decreases with increasing distance.In the central region, labeled as core region, the minimum inter-node distance is about 2.2 km.The model is shown in fig.1a,b.In the following, we will show some results and benchmarks for a homogeneous elastic medium, but the generalization to a heterogeneous system is strightforward, as described in Section 5.The Poisson coefficient is fixed at ν = 0.25 and the Lamé constants at λ = µ = =3.2⋅10 10 N/m 2 .The fault plane measures about 50×30 km and is placed within the core region, with the top of the plane located at a depth of about 10 km under the surface.The total seismic moment is set at M 0=1⋅10 19 N⋅m.Boundary conditions must be carefully chosen to avoid unwanted edge effects, the FE method being intrinsically limited to manage with finite domains.We experimented with a number of boundary conditions commonly used in the Literature.These include for example to specify zero displacement along the bottom and lateral surfaces delimiting the volume (Masterlark et al., 2001) as well as to keep nodes on these boundaries fixed in the direction perpendicular to the surface itself (Cianetti et al., 2005;Megna et al., 2005).As an example, figs.3a-e to 10a-e show the ground horizontal and vertical displacement field for vertical and tilted faulting sytems when the cited boundary conditions are imposed, compared to the Okada solutions.In the pictures, displacement directions are rendered using a regular array of vector arrows having the same length and color plots are used to exhibit and compare moduli.From the comparison between displacement directions and magnitudes it can be clearly seen that none of the imposed conditions is able to simulate an infinite domain removing any edge effect.As far as the horizontal displacements are concerned, the strongest disagreement is observed in the vector plots, where directions do not coincide near the edges, although in some cases a slight difference in the distribution of displacement magnitude can also be noticed.
The differences are even more evident for vertical displacements, both in magnitude and especially in versus.Actually, it is clear that edge ef-fects tend to vanish departing from the boundaries, but the use of not optimized boundary conditions forces us to deal with larger domains, to the detriment of the computational efficiency, especially in case of non local applications.This evidence motivated us to seek an alternative strategy as follows.Analytical models, as mentioned in the Introduction (Section 1), are able to supply the exact solution for the diplacement field generated by a faulting source within an infinite homogeneous domain, when the rhe-ological properties and the fault plane size and orientation are known.Starting from this solution, and in particular from the Okada solution, we can calculate the displacement each bottom and lateral node would have if the domain was infinite and impose these prescribed displace-  ments as inhomogeneous boundary conditions.This option is formally correct as long as rheologically homogeneous systems are treated, but it obviously represents an approximation if discontinuities and/or heterogeneities are intro-duced.In this case comparative tests should be performed in order to tune the best boundary condition set for each specific case.Simulations for several test cases were performed, taking into account several different  strike and dip slip sources and applying for all of the cases the so-called Okada boundary conditions.In the following, we will exhibit some selected results.Note that, due to our way of dealing with the strike angle, the domain ap-pears rotated in figures concerning faults with φ ≠ 0°.Of course, to obtain plots identical to the former, we would only need to start from a larger domain which can be suitably cut out after reference frame rotation.Vertical faulting represents the simplest system to study, as no mesh distortion due to node displacements is introduced and the FE simulation is carried out on the original mesh.It is worth noting that vertical dip slip has been considered for the sake of completeness, since it does not make great physical sense, being almost absent in nature.We investigated vertical strike and dip slip faults with strike angle φ = = 0°, 30°.Figures 11a-e, 12a-e, 13a-e and 14a-e show the ground horizontal and vertical displacement field produced by the strike and dip slip, respectively, compared to the Okada solutions.It can be clearly seen that the agreement is good: if some difference is observed it is just very near the source, where our mathematical  approximation, based on a finite distribution of double couples, is less accurate.Solutions on tilted faults with strike angle equal to 0°and 30°have been also calculated.Figures 15a-e and 16a-e show the ground horizontal and vertical displacement field produced by a strike slip with dip angle δ = 60°and φ = 0°, 30°, respectively.Similarly, figs.17a-e and 18a-e show the ground horizontal and vertical displacement field produced by a dip slip with dip angle δ = 20°and φ = 0°, 30°, respectively.Also in these cases the agreement between the numerical and the analytical solution is very good, indicating that numerical errors due to mesh distortion do not affect the final result.

Illustrative applications
One of the most appealing features of modeling static and quasi-static seismic deformation by means of the FE method is the possibility of dealing with rheological and/or geometrical complexities.Using FEMSA, the introduction of 3D heterogeneities within the domain basically requires the creation of a suitable geometrical model (for example to model a subducting slab in correspondence of a convergent plate boundary) and/or a differentiated definition of the mechanical properties.In the following, we present some simple examples of applications with the purpose of illustrating some FEMSA capabilities.We start considering a simple model in which only rheological heterogeneities are taken into account and perform both a static and a quasi-static analysis.
The schematic model illustrated in fig.19 represents a border zone between the oceanic and the continental crust.The usual hexahedral domain is split into three regions, i.e. 1) the oceanic crust, 2) the continental crust and 3) the mantle, having different rheological properties.Practically speaking, the mesh elements are assembled into three element sets, each having different material properties assigned.Firstly, each region is assumed to exhibit a purely elastic behaviour.In this case, the rigidities µ1, µ2 and µ3 are assigned in such a way that µ1/µ2≈ ≈0.6 and µ2/µ3≈0.8 and that their average µ= =(µ1+µ2+µ3)/3 equals the rigidity value previously used for the homogeneous system.The Poisson coefficient is fixed at ν =0.25 and the Lamé constants are set to be λ =µ, as in the homogeneous system considered before.Secondly, we introduced linear Maxwell rheology into the model, by assigning creep properties to the mantle zone, with a viscosity value η=1⋅10 19 Pa⋅s, and keeping the two crustal zones as purely elastic materials.
We considered a shallow vertical right-lateral strike slip fault within the oceanic crust.Inhomogeneous boundary conditions have been applied to calculate the coseismic deformation, while as far as the postseismic relaxation is concerned we opted to keep nodes on the boundary of the domain fixed in the direction perpendicular to the surface.The quasi-static analysis was performed for a total time period of 100 years with a time step of 20 years.In fig.20a-c the results of the two analyses are displayed and compared.Viscoelastic relaxation induces an evident amplification of both the horizontal and the vertical displacement components, the amplifying factor coming larger with increasing distance from the source.While the viscoelastic contribution shows generally the same sign as the elastic residual for the horizontal components, this is not everywhere true for the vertical displacement, which exhibits a more complex pattern of viscoelastic relaxation.All these behaviours are well known for laterally homogeneous domains and are in agreement with several previous results obtained by means of semi-analytical models (e.g., Piersanti et al., 1995;Antonioli et al., 1998;Nostro et al., 1999).
Noteworthy is that the presence of the rigidity contrast within the domain largely affects the deformation field, producing a macroscopic asymmetric pattern which appears especially pronounced in the horizontal component.As expected, the presence of the oceanic crust amplifies the displacements.The joint effects of viscoelastic relaxation and lateral heterogeneity on the horizontal components are somewhat straightforward: the asymmetry induced by lateral heterogeneity is further increased by the relaxation, which is more vigorous under the oceanic region where the crust is thinner (i.e. the viscoelastic material is nearer).Though the general asymmetric features are conserved, for the vertical compo-nents the pattern is more complex due to the trend of vertical viscoelastic relaxation that not always exhibits the same sign of the elastic offset.
Finally, we present a very simple spherical model of the elastic residual deformation associated with the giant Sumatra-Andaman earthquake of December 26, 2004.
The Sumatra earthquake resulted from complex slip on the fault along the subduction zone where the oceanic portion of the Indian Plate slides under the Eurasian Plate, by the Indonesian island of Sumatra.The direction of convergence of the subducting plate relative to the overriding plate is oriented oblique to the trench axis and the rupture occurred for 1200 km along the interplate megathrust.
Since the investigation of the deformation field produced by this event requires a very long range analysis, we built up a spherical domain consisting of a portion of spherical zone of about 1000 km thick spanning about 90×10 6 km 2 on the Earth surface (see fig. 21).A layered model having elastic rheology was introduced, based on the volume averaged mean values of the Lamé parameters according to the Preliminary Reference Earth Model (PREM) (Dziewonski and Anderson, 1981).The domain was meshed with 38348 20-nodes brick elements resulting in a mesh containing 171537 nodes.Such a mesh takes 44 min to perform a static calculation on a 64-bit Intel Xeon 3.2 GHz, 8 GB RAM.As far as the source is concerned, we Standing on these pictures, on the spatial scale involved in our analysis, the fault geometry (as far as a homogeneous moment release is considered) does not seem to be crucial in defining the displacement field generated by the earthquake.The comparison with GPS data registered from permanent stations, as shown in fig.25, is encouraging despite the roughness of this illustrative model.In particular, with respect to previous semianalytical applications (Banerjee et al., 2005;Boschi et al., 2006), we obtain a fair agreement with the GPS stations in India without needing to introduce any heterogeneity in the seismic moment release.

Conclusions
In the present work, we developed a numerical simulation tool, FEMSA (Finite Element Modeling for Seismic Applications), to automatically implement arbitrary faulting sources and perform both the analytical and numerical calculation to compute the horizontal and vertical displacement and stress fields induced by the fault slip.Our approach is based on Cal-culiX, free software designed to solve field problems by means of the FE method.The main advantages of such a method are: reliability, wide diffusion and flexibility, allowing geometrical and/or rheological heterogeneities to be included in a mechanical analysis.FEMSA is freely distributed to the community.
We carried out an optimization study on boundary conditions as well as a series of benchmark simulations on test cases, comparing our numerical results with the analytical ones.Our investigations induce to assert that the best boundary conditions are obtained when the analytically calculated displacements are prescribed on nodes located on the bottom and lateral edges of the domain.Applying this method, we obtained a very good agreement between numerical and analytical results both for strike and dip slip faults with different values of the strike and dip angles.We also verified the capability of our approach to tackle some heterogeneity within the domain, showing some simple preliminary applications.
The introduction of some real Earth complexities, like more heterogeneous domain or complex anelastic rheology, is pretty simple and can be done with relatively small effort by the user.Anyway FEMSA will be continously developed and improved.Upgraded versions of the package, accompanied by the user's manual, will be made available as new potentialities are implemented (e-mail to volpe@ingv.it).

Fig
Fig. 1a,b.a) Perspective and b) top view of the domain: the central volume is the core region.

Fig. 2 .
Fig. 2. Block diagram of the automatic procedure to numerically calculate the displacement field induced by faulting sources, as described in the text.

Fig.
Fig. 3a-e.Ground displacement field for a right-lateral strike slip with dip δ=90°and strike φ=0 °when the the bottom and lateral nodes are fixed.Horizontal displacement field (on the left): a) numerical and analytical results are displayed simultaneously as a vector plot (blue arrows: numerical solution; red arrows: analytical solution); b) numerical and c) analytical solutions are displayed as color plots.Vertical displacement field (on the right): d) numerical and e) analytical solutions are displayed as color plots.

Fig
Fig. 4a-e.Same as fig.3a-e but for a right-lateral strike slip with δ=60°and φ=0 °when the the bottom and lateral nodes are fixed.

Fig.
Fig.5a-e.Same as fig.3a-e but for a right-lateral strike slip with δ=90°and φ=0 °when the bottom and lateral nodes are fixed in the direction perpendicular to the surface.

Fig.
Fig. 6a-e.Same as fig.3a-e but for a right-lateral strike slip with δ=60°and φ =0 °when the bottom and lateral nodes are fixed in the direction perpendicular to the surface.

Fig.
Fig. 7a-e.Same as fig.3a-e but for a normal dip slip with δ=90°and φ=0 °when the the bottom and lateral nodes are fixed.

Fig.
Fig. 8a-e.Same as fig.3a-e but for a normal dip slip with δ=20°and φ=0 °when the the bottom and lateral nodes are fixed.

Fig
Fig.9a-e.Same as fig.3a-e but for a normal dip slip with δ=90°and φ=0 °when the the bottom and lateral nodes are fixed in the direction perpendicular to the surface.

Fig
Fig.10a-e.Same as fig.3a-e but for a normal dip slip with δ=20°and φ =0 °when the the bottom and lateral nodes are fixed in the direction perpendicular to the surface.

Fig. 11a -
Fig. 11a-e.Ground displacement field for a right-lateral strike slip with δ=90°and φ=0 °when the Okada displacements are applied as inhomogeneous boundary conditions.Horizontal displacement field (on the left): a) numerical and analytical results are displayed simultaneously as a vector plot (blue arrows: numerical solution; red arrows: analytical solution); b) numerical and c) analytical solutions are displayed as color plots.Vertical displacement field (on the right): d) numerical and e) analytical solutions are displayed as color plots.

Fig. 20a -
Fig. 20a-c.Ground coseismic and postseismic displacement field for a right-lateral strike slip with dip δ=90°and strike φ=0 °embedded in the heterogeneous plane domain illustrated in fig.19.Horizontal displacement field: a) coseismic (blue arrows) and postseismic after 100 years (red arrows) are displayed simultaneously as a vector plot.Vertical displacement field: b) coseismic and c) postseismic after 100 years are displayed as color plots.

Fig. 21 .
Fig. 21.Pictorial view of the spherical domain built up to investigate the 2004 Sumatra earthquake.

Fig. 22 .
Fig. 22. Ground horizontal (vector plot on the top) and vertical (color plot on the bottom) displacement field calculated for the 2004 Sumatra earthquake with a point source approximation.

Fig. 23 .
Fig. 23.Ground horizontal (vector plot on the top) and vertical (color plot on the bottom) displacement field calculated for the 2004 Sumatra earthquake with a line source (1D) approximation.

Fig. 24 .
Fig. 24.Ground horizontal (vector plot on the right) and vertical (color plot on the left) displacement field calculated for the 2004 Sumatra earthquake with a plane source (2D) approximation.

Fig. 25 .
Fig. 25.Comparison between horizontal displacements resulting from the FEM simulation and the GPS measurements recorded from permanent stations (here plotted with error ellipses matching 60% confidence).

Table I .
Brief