Regional wave propagation using the discontinuous Galerkin method

Introduction Conclusions References


Introduction
Understanding the influence of 3-D Earth structure on regional seismic wave propagation requires numerical algorithms that accurately simulate ground motions towards high frequencies on today's computing architectures.The schemes have to flexibly handle complicated geological formations such as the structure of subduction and rifting zones, the shape of mantle plumes, the topography of elastic discontinuities, or the heterogeneities of continental crust.In this study, we present an application of the discontinuous Galerkin (DG) method to regional scale earthquake simulations, offering considerable flexibility to accurately take into account 3-D Earth models, particularly in the meshing part of the workflow.The first regional to global scale simulations, based on spherically symmetric Earth models, were obtained by the summation of normal modes (Alterman et al., 1959).The approach was extended, to account for effects of 3-D Earth structures and localized heterogeneities using perturbation theory (e.g., Woodhouse, 1983).Besides semianalytical algorithms, numerical direct solvers evolved such as the Finite-Difference (FD) method.FD approximates solutions to the wave equation on a structured grid using low-order (e.g., Virieux, 1984), and later also high-order (e.g., Bayliss et al., 1986) numerical operators in 3-D, in order to improve the accuracy and computational efficiency.Due to restrictions in the geometrical flexibility of standard FD methods using structured grids with regular grid spacing, later also approaches using variable grids have been introduced (Pitarka, 1999;Ely et al., 2008;Kristek et al., 2010;Kozdon et al., 2012).
As an alternative, Finite-Element (FE) methods have been applied in seismology (e.g., Bao et al., 1998).The techniques make use of orthogonal basis functions to approximate the solution in finite subspaces (elements) of the modeling domain.This concept enables the use of various element types (e.g.tetrahedra, hexahedra, pyramids, prisms and combinations of them) to flexibly discretize the computational domain.
Motivated by the difficulties in mesh generation for complex geometries using structured or unstructured deformable hexahedral grids, the discontinuous Galerkin method was introduced to extend the modeling capabilities to unstructured tetrahedral grids (e.g., K äser and Dumbser, 2006;Dumbser and K äser, 2006;K äser et al., 2007;de la Puente et al., 2007;Dumbser et al., 2007, andEtienne et al., 2010).To honor complex Figures 3-D shapes, tetrahedra offer the largest flexibility, because the elements can easily be aligned to small scale structures of the Earth, and the size of the elements can locally be adapted (h-adaptivity) without overhead.Furthermore, the adaptation of the polynomial degree inside each element (p-adaptivity) enables the computational load to be focused on areas of interest.In addition to the spatial flexibility, the accumulation of errors, due to numerical dispersion, needs to be minimized.Therefore, an arbitrarily high-order derivative (ADER) scheme (Toro, 1999) was combined with the DG method for elastodynamic problems to achieve the same accuracy for the temporal as for the spatial approximation of the wavefield.The ADER-DG approach supports local time stepping schemes to balance the varying computational costs of elements with different size and polynomial order (Dumbser et al., 2007).
In this paper, we show that the ADER-DG method can serve as an alternative algorithm to established numerical methods for regional scale earthquake simulations.In particular, we investigate the applicability of the approach in terms of mesh generation and simulation accuracy.
The paper is structured as follows.In Sect. 2 we describe the proposed ADER-DG approach for the numerical simulation of wave propagation in 3-D media.In addition to the technical properties of the implementation, we explain in detail the mesh generation process for geometrically complex 3-D media.In Sect. 3 we verify the radiation of elastic waves excited by an explosive and a shear dislocation source.Synthetics of the SE and our ADER-DG approach calculated on a 1-D isotropic PREM model (Dziewonski and Anderson, 1981)

The ADER-DG approach
To give an introduction to the numerical scheme of the ADER-DG method, we will qualitatively describe the features and usability of the approach with a special focus on the model discretization.For details on the theoretical framework and the implementation of the method, the reader is referred to the work of K äser and Dumbser (2006) et al. (2007).
As every Finite-Element method, the DG method solves the weak form of a partial differential equation, here the elastic wave equation.Therefore, the wave equation is multiplied by a test function, and integrated over a finite element subspace.In the scheme the complete 3-D computational domain can be divided into a mesh of unstructured tetrahedra.The unknown solution is approximated inside each element by a linear combination of space dependent polynomial basis functions of arbitrary degree N and its purely time dependent coefficients, the so-called degrees of freedom (DOF), which advance in time.It is referred to as a Galerkin scheme, because test functions are chosen from the same basis function space.We use an orthogonal basis suggested by Dubiner (1991).
In contrast to classical FE or SE schemes, the polynomial representation of the physical solution is not forced to be continuous across element boundaries in DG.At this point, the DG approach adopts the well-known concept of fluxes to exchange information between neighboring elements, as defined in the Finite-Volume (FV) approach.The flux terms are computed via the upwinding exact Riemann solver, known as Godunov flux (Toro, 1999).Note, that differently from FV only adjacent elements communicate, and a reconstruction process is not required.Furthermore, the amount of data necessary for the flux computation reduces to a rather small number of DOF, so that in the case of a parallel run only minor MPI communication occurs, and the code is highly scalable (K äser et al., 2010) The method can be implemented efficiently by transforming each individual shaped and located tetrahedron from the physical space into a uniform reference element, since many occurring integrals can be precomputed in the reference space.The quadrature free DG discretization is advanced in time by the arbitrary high-order derivatives (ADER) approach of Toro (1999).A Cauchy-Kovalewski procedure is applied to recursively replace the temporal derivatives in the Taylor expansion by spatial derivatives.In this way, we obtain automatically the same order of accuracy for the time integration as well as for the spatial approximation, while allocating the same memory as a first-order explicit Euler time integration.The code has been proven to achieve spectral convergence (Dumbser and K äser, 2006).

Local adaptivity and load balancing
In our implementation of the ADER-DG method, unstructured tetrahedral grids can easily be refined to local structures (h-adaptivity).This increases the spatial sampling of geometrical complex structures or the wavefield itself.Depending on the problem and the refining strategy, the size of the elements can vary tremendously within the computational domain, and corollary the time increment according to the Courant-Friedrichs-Lewy (CFL) condition can change dramatically between the elements.In order to reduce the computational effort and balance the different time step lengths, the polynomial degree can be adapted in each tetrahedral element individually (padaptivity).Furthermore, computational efficiency can be enhanced using an optimal time step length for each individual element (local time stepping).This enables the use of elements with completely different size without affecting the high-order accuracy of the method (Dumbser et al., 2007).
On the other hand, asynchronous element updates due to local time stepping combined with a p-adaptive mesh can cause irregular load balancing on parallel machines.Introduction

Conclusions References
Tables Figures

Back Close
Full tasks.This aims at an equal distribution of small and large elements over all processing units to balance the work load (K äser et al., 2010).Since the zones are assigned randomly to each task, exchanging updated field values can increase MPI communication.However, compared to the numerous calculations inside each domain and the small MPI communication effort for DG, this overhead is small (Dumbser and K äser, 2006).

Mesh generation in geometrically complex 3-D media
Geoscientific models include geometrically complex 3-D structures such as rough surface topography, as well as complicated undulating interfaces at subsurface velocity contrasts, which have tremendous influence on the propagation of seismic waves.

Structured meshing
An approach of discretizing a geological subsurface model is using regular-spaced structured grids that are computationally very efficient to implement.However, this discretization strategy results in a staircase approximation of the Earth's surface or complex internal interfaces, if material properties are simply assigned to nodes or elements.Such an approximation of material discontinuities leads to first-order errors, proportional to the grid spacing or time step, that are insensitive to the approximation order of the numerical scheme (Gustafsson and Wahlund, 2004).To achieve high accuracy of the simulation, the regular grid spacing is forced to be very small, leading to an increasing computational effort especially for strong material heterogeneities (Bohlen and Saenger, 2006;Pelties et al., 2010).There are more sophisticated approaches to incorporate material interfaces using structured grids.For example, Kristek and Moczo (2003) are using values of effective grid material parameters to account for discontinuities.In contrast, Komatitsch (1997) assign the material properties directly to single integration nodes in the SE scheme, at least, if interfaces cannot be respected within Figures the computational mesh.However, the simple approach we described previously, is still widely used.
For regional wave propagation problems in 3-D, numerical algorithms evolved solving the wave equation efficiently in spherical coordinates on smaller sections of the Earth using curved hexahedral grids (Igel, 1999;Igel et al., 2002;Fichtner and Igel, 2008;Fichtner, 2009).The spherical discretization leads to singularities in the solution of the wave equation, which inhibits an application to fully global simulations.Therefore, the Cubed-Sphere approach (Ronchi, 1996) was introduced, which became very successful, in particular, in combination with the SE method (Komatitsch and Tromp, 2002).
Hereby, the projection scheme is based on mapping a cube, divided by an equidistant quadrilateral grid, to an angular-distant quadrilateral grid at the surface of a sphere.In radial direction interpolation produce hexahedra with constant size, while the lateral size decreases towards the center of the sphere due to the spherical geometry.This results in an unnecessary oversampling of the wavefield since, generally, velocity values increase with depth.To avoid this problem, the mesh is re-gridded to double the mesh spacing at different radial depths.Furthermore, grid doubling techniques have the potential to improve the efficiency as the total number of elements can be reduced.However, grid doubling is difficult to implement and can produce undesired numerical artifacts like instabilities or spurious reflections at the grid interfaces.

Unstructured meshing
Unlike structured meshes, unstructured meshes are composed by an irregular pattern of elements, which are not arranged in a logical sequence.The elements have to be defined by a list of vertex coordinates and its connectivity to the elements, which increases the algorithmic complexity of corresponding numerical approaches.Furthermore, assembling a scheme that is as accurate and efficient as a method implemented on structured grids, is challenging.But on the other hand, a strong flexibility in mesh generation of complicated 3-D structures can be gained.In general, unstructured grids for complex structures cannot be obtained automatically and a meshing tool is required.An external meshing software such as Cubit (Blacker et al., 1994) can be employed to generate unstructured hexahedral grids for example for regional seismological applications, see Casarotti et al. (2008b); Stupazzini et al. (2009); Peter et al. (2011);Cupillard et al. (2012).
The generation of high-quality unstructured hexahedral meshes, as applied in the examples given above, often requires a strong user interaction, which is time consuming and cumbersome.FE and FV methods, for example, can make use of fully unstructured grids with various element types, but are usually only in low-order formulations efficiently implementable that are very dispersive.This motivated the development of numerical methods that can handle unstructured tetrahedral meshing approaches such as the high-order accurate DG method.
Especially for arbitrary shaped 3-D models including, for example, interface and surface topography, tetrahedral grids are much more flexible to align.This tremendously reduce the meshing effort potentially at the expense of longer simulation time.

Workflow of tetrahedral mesh generation using Cubit
The general workflow to generate computational meshes can be divided into the geometry setup and the meshing process itself.For both parts, we are using the commercial software Cubit.This meshing tool contains numerical libraries to represent the underlying geometry of a computational domain using Non-uniform rational B-splines (NURBS) that can accurately describe arbitrary geometrical shapes.Furthermore, it provides a flexible tetrahedral meshing algorithm to generate high-quality meshes.
In the following subsections, the geometry and mesh generation workflow using Cubit is demonstrated for an Earth model of Europe (Sect.4), composed of the 3-D EPcrust model (Molinari and Morelli, 2011) imposed on the depth dependent ak135 velocity model (Kennett et al., 1995).Introduction

Conclusions References
Tables Figures

Back Close
Full

Geometry generation -European model
Recent studies investigated, whether internal material discontinuities of complex 3-D structural models should be geometrically respected by the mesh or not (Komatitsch et al., 2004;Lee et al., 2008;Casarotti et al., 2008a;Stupazzini et al., 2009;Cupillard et al., 2012).Not respecting material interfaces of strong contrasts requires a high resolution for ADER-DG schemes, similar as seen for structured discretizations.Consequently, the computational costs are unnecessary expensive, since the dispersion properties would allow for larger elements.In order to reduce the runtime, improve the accuracy of model approximation, and utilize the flexibility of tetrahedral meshing, we respect interfaces such as the Moho or upper-mantle discontinuities in the European model.
For 3-D models, these interfaces often are given by structured or unstructured point clouds, which have to be composed to parametric surfaces.For example, the 3-D EPcrust model is parametrized by angular-distant points given on a structured grid of (0.5 • ×0.Cubit.For large datasets these steps can be automated using the Python-interface of Cubit.
To obtain a closed volume for each layer of the model, the surfaces are connected by line segments to generate the boundary faces of the domain.Furthermore, so-called imprinting surfaces have to be defined to generate conforming meshes at each layer boundary.
As part of the workflow, the geometry generation is summarized by the first four boxes in the flow-chart Fig. 1.The frames highlighted in gray illustrate the surface reconstruction, which was performed for each model interface separately.
As already mentioned the pointsets representing subsurface structures of the Earth can be arbitrarily distributed.Therefore, it is necessary to use robust surface reconstruction algorithms.Computer-Aided-Design (CAD) tools focused on geometry generation may be an alternative to generate surfaces directly from pointsets, but still, in most cases these are commercial packages.

Tetrahedral mesh generation -European model
Mesh generation and in particular the size of elements inside different areas of the computational domain depend on three main factors: (1) the wavefield approximation, (2) the approximation of the material model, and (3) on the approximation of the geometrical model.Cubit provides geometry adaptive sizing functions to control the element size and its growth rate in a volume during meshing.This allows for a local mesh refinement or coarsening, considering a well balanced aspect ratio of the elements.In addition, quality measures can provide information on the shape of the elements.Statistics on the mesh elements is necessary, because element distortions can lead to very small time steps increasing the computational costs or even lead to instabilities interrupting the numerical simulation.Introduction

Conclusions References
Tables Figures

Back Close
Full The right part of Fig. 4 shows the result of the mesh generation for the European model.Meshing tests of different constant element sizes have shown a strong dependency of the mesh quality on the chosen edge length of the elements, even if little changes to the element size were applied.This was particularly apparent for crustal meshes, since the ratio of the element edge length and the layer thickness is small.
In case of distorted elements in the mesh, re-meshing with a different element size or manual changes of the mesh nodeset improved the mesh quality.
For a numerical simulation, boundary conditions have to be assigned to the corresponding element faces.Together with a list of the connectivity between vertices and elements, the list of boundary elements is output and converted into a mesh file used by our ADER-DG solver SeisSol.The mesh generation is summarized by the last three boxes in the flow-chart Fig. 1.
We emphasize that despite the mentioned difficulties the generation of a high-quality tetrahedral mesh took one day for the given example, which clearly indicates the flexibility of unstructured tetrahedral meshing for complex 3-D models.Furthermore, a strategy to provide a mesh generation workflow using non-commercial software, has further to be elaborated.

Benchmarks
To verify the accuracy of our ADER-DG solver SeisSol on a regional scale, synthetic seismograms are compared to solutions of the well established SE method implementation SpecFEM3D-Globe.For regional to global scale simulations SpecFEM was successfully benchmarked against quasi-analytical normal mode solutions (Komatitsch and Tromp, 2002).
In the first experiment we investigate elastic wave propagation excited by an explosive source.Afterwards, the radiation of seismic waves is studied for a shear dislocation For the simulations, in both schemes, a spatial approximation order of O = 5 is chosen.Whereas the ADER approach provides the same order of spatial and temporal accuracy, the SE method uses a time integration scheme of second order accuracy.Ground displacements calculated by SpecFEM are differentiated in time to compare them to the velocity output of SeisSol.With respect to the results of the SE method a seismogram difference can be quantified by the time-frequency misfit analysis introduced by Kristekova et al. (2006).Here, the signals are decomposed in the timefrequency plane to characterize the difference by separation of phase and envelope misfits.

Setup
In this test we analyze the numerical radiation of seismic waves excited by an explosive source using SeisSol.The point dislocation is represented by a moment tensor where the only non-zero entries are M xx = M yy = M zz = 10 18 Nm.The origin of the coordinate system is located in the center of the Earth at (x 0 , y 0 , z 0 ) = (0 km, 0 km, 0 km).The of T min = 20 s.The high spatial discretization of the mesh is also suited to model wavefields over long propagation distances.SpecFEM provides a mesh for the whole Earth adapted to the PREM model.The number of spectral elements at the Earth's surface, along one side of a Cubed-Sphere chunk (1/4 of Earths perimeter), is set to 384.Due to the spherical geometry the mesh is doubled in size three times towards the Earth's core.According to benchmark results of Komatitsch and Tromp (2002) the setup allows for the simulation of T min = 12 s as a shortest wave period.
At five stations, which are located at latitudes from 89 • N to 81 • N and longitude 0 • W, seismograms of 800 s are output.Data of SeisSol and SpecFEM are low-pass filtered at a peak frequency of 0.05 Hz using an order three Butterworth-filter.

Results
Figure 2 shows synthetics of the radial (v r ) and vertical velocity component (v v ) as a function of time.For each station the traces are normalized with respect to the v v component, which is excited most strongly by the explosive source.
The signals of body and surface waves simulated by SeisSol and SpecFEM match almost perfectly at all stations.This visual impression is confirmed by low values of phase and envelope misfits.Envelope misfits vary around 0.9 − 1.5 % while phase misfits do not exceed the level of 1 %.The lower phase misfits indicate a lower sensitivity of SeisSol to numerical dispersion effects but stronger numerical dissipation, which was already determined by K äser et al. ( 2008) who compared the algorithm to analytical solutions.Note, the synthetics of SeisSol and SpecFEM are both be prone to phase and amplitude misfits within their numerical accuracy.
In case of an explosive source excitation, transverse motion (not shown) should be recorded only at the numerical noise level, which holds for the SE method simulation.be caused by the application of unstructured grids and will be discussed more in detail in Sect.3.3.

Setup
The main purpose of this simulation is to verify accurate numerical wave radiation in SeisSol using a shear dislocation source, again applied to the 1-D isotropic PREM model.We chose a magnitude 6.0 event, which occurred on the 20th of August 2009 (06:35 GMT) located in the Norwegian sea at latitude 72.22 • N and longitude 0.84 • E at a depth of 12 km.The six independent components of the centroid moment tensor (CMT) solution are given in a spherical reference system provided by the CMT catalog of the Lamont-Doherty Earth Observatory, Columbia University, see Table 1.
For the SeisSol simulation we reduced the block model of Sec.3.1 to a spherical domain of with a depth range of 500 km.Except for the source characteristics, SpecFEM results correspond to the same setup as in the previous example.The synthetic seismograms of 350 s length are processed using the same STF and low-pass filter as specified in Sec.3.1.

Results
Figure 3 illustrates the velocity components plotted as a function of time at stations northward of the source.Here, amplitude values are normalized to the strongest transverse velocity component v t .As expected the signals of SeisSol and SpecFEM are almost identical.Except for near offset stations the envelope and phase misfits have largest values of 1%, where in general, as for the first experiment, the envelope misfit level is higher than the phase misfit.The large envelope misfit value of 2.9 % at station X73N.0, which is closest to the source, could be due to a different excitation characteristic of tetrahedral compared to hexahedral element.Introduction

Conclusions References
Tables Figures

Back Close
Full

Discussion
Accurate simulations of elastic wave propagation using the ADER-DG approach could be verified applying different source mechanisms on regional sections of the 1-D isotropic PREM model.The overall match of waveforms simulated by SeisSol compared with SpecFEM results are excellent.
However, at far offset stations a remaining, very low envelope and phase misfit of around 1 % is observed.This could be due to the difference in the approximation of material parameters within SeisSol and SpecFEM.In SeisSol a constant value for each material parameter is interpolated in one single tetrahedral element, whereas in SpecFEM material values are evaluated at each GLL point within a hexahedron leading to a higher (subcell) sampling of the PREM background model.Castro et al. (2010) demonstrated that an element-wise material sampling using constant parameters is not an intrinsic restriction of the DG approach.But, tests have shown that the spatial discretization already is determined by an accurate approximation of the wavelength due to the CFL-condition.Therefore, the sampling of the relatively low material gradient in PREM is sufficient.However, the different sampling of the material values could explain small differences in the results.
At near offset stations larger envelope misfits occurred, which are thought to be caused by the radiation of the source pulse within an asymmetric tetrahedral element.This is connected to the radiation of SH-energy obtained for an explosive source in the SeisSol simulations, which is not allowed in theory.K äser and Igel (2001) claimed that for an explosive source excited in an isotropic elastic medium the application of numerical operators on unstructured meshes can lead to S-wave artifacts at the transverse velocity component v t .Convergence tests we carried out in an isotropic elastic medium have shown a decreasing v t signal amplitude with decreasing mesh spacing, i.e. the effect is purely numerical and can be diminished by refining the mesh around the source.Introduction

Conclusions References
Tables Figures

Back Close
Full

Comparison of code performance
The performance of the numerical implementations was estimated by a theoretical calculation of de la Puente et al. (2008).The analysis is reduced to the number of matrix operations O of SeisSol and SpecFEM, which led to an estimation of O DG ≈ 2.8•O SE at a polynomial degree of N = 3.As mentioned in the thesis, one has to be aware that the performance of a software additionally depends for example on the programming, the size of the time step, number of elements, number of field variables or the cost of the time integration.Since the simulations of this study provide a different level of accuracy on different meshes and physical domains a quantitative comparison is not possible, in fair terms.
To quantitatively analyze the efficiency of numerical codes a conceptual basis has to be defined depending on an accurate definition of the problem to be solved and the criterion of comparison.For example, errors due to the order of spatial and temporal discretization in the numerical algorithm can be independently dominating depending on the modeling parameters such as the applied Earth model or the wave propagation distance.These tradeoffs must be explored separately and will be addressed in a future study to evaluate the ADER-DG and SE code efficiency.

Application of the ADER-DG method to real data: The 2009 L'Aquila earthquake
In the previous tests it could be demonstrated, that SeisSol is able to accurately simulate wave propagation on layered 1-D Earth models like PREM.In this section we qualitatively compare SeisSol synthetics with real data of the magnitude M W = 6.3 L'Aquila event (central Italy) to examine the applicability of the ADER-DG method on modeling real earthquake scenarios using 3-D Earth models.Introduction

Conclusions References
Tables Figures

Back Close
Full The earthquake source area evolved by subduction of the Adria micro-plate, by continental collision of Eurasia and Africa and by the opening of the Tyrrhenian back-arc basin.Therefore, the region of the central Apennines is tectonically very complex, and thus has one of the highest seismic hazard in Italy (Chiarabba et al., 2009).
On the 6th of April 2009 (01:32 GMT) the main shock of a seismic sequence struck the Abruzzi region due to normal faulting along a NW-SE striking fault system at a depth of about 12 km.The earthquake devastated the city of L'Aquila and surrounding villages, leading to 308 deaths, about 1500 injured and an economic loss of 15 Billion Euro.

Model setup
To consider a realistic representation of the Earth's crust we introduce lateral variations of the elastic material values using the recent 3-D model EPcrust of Molinari and Morelli (2011).This model is imposed on top of the 1-D reference Earth model ak135 of Kennett et al. (1995) to form a composite model of Europe for the ADER-DG simulations.
The EPcrust model was deduced from literature information of different global and local models using active-source seismic profiles, receiver function studies and digital maps (e.g., Bassin et al. (2000); Tesauro et al. (2008)).It is divided into 3 layers (sediments, upper and lower crust) of variable thickness, which cover the whole European Plate on a geographical grid of 0.5 • ×0.5 • .Surface topography is incorporated by the ETOPO1 model of Amante and Eakins (2008).The geologic formations are parametrized on curved layer interfaces by the isotropic P-, S-wave velocity and density.The P-wave velocity varies between 1.5 < v P < 7.2 km s −1 , while, due to scaling relations of Brocher (2005), the S-wave velocity and density range by 0.4 < v S < 4.1 km s −1 and 0.9 < ρ < 3.0 g cm −3 , respectively.
The 1-D reference model ak135 (Kennett et al., 1995) was constructed from a set of empirical traveltime curves of relocated earthquake events.It provides good Figures fits to a variety of seismic phases, but does not respect the effect of lateral heterogeneities inside the Earth.In our European model, ak135 supplies the upper-mantle (35−660 km) material values for P-wave velocity 8.0 < v P < 10.2 km s −1 , S-wave velocity 4.5 < v S < 5.6 km s −1 and density 3.3 < ρ < 4.1 g cm −3 .

Geometrical representation and mesh generation
The physical representation used for the ADER-DG simulation is spherically bounded by 20 • N−60 • N latitude and 10 • W−30 • E longitude, shown in the left sketch of Fig. 4.
In this domain, the topography of the Earth surface and the Mohorovicić discontinuity (Moho) are represented by 3-D splines.They form the top and bottom of the crustal layer with a thickness of 6.3 − 52.1 km, while interior interfaces of the upper and lower crust are not respected.The upper mantle region is composed of three layers bounded by smooth spherical interfaces according to the elastic discontinuities in ak135 at 210, 410, 660 km depth.
To generate a mesh inside the volume an appropriate mesh size has to be defined.The specification depends on geometrical restrictions of the physical domain, the prospective approximation accuracy of the material model and the spatial sampling of the wavefield.Due to K äser et al. (2008), who carried out an ADER-DG convergence study on a simple block model, we expect high simulation accuracies over regional propagation distances at a spatial sampling of 2 tetrahedral elements per smallest wavelength, if an O = 5 scheme is applied.In the mantle the mesh size is adapted according to the 1-D velocity distribution of ak135 to ensure a spatial sampling of 2 elements per shortest wavelength.For the whole mesh, this results in a total number of 3.7Melements leading to 1.164 M degrees of freedom for an O5 ADER-DG simulation.On an IBM system using 50 Intel Xeon Westmere-EX processors (10 cores, 2.4 GHz) a runtime of 24 h was necessary to output seismograms of 800 s.

Data processing
For the L'Aquila earthquake broadband data of a high signal-to-noise ratio are provided by several seismological networks such as IRIS, GEOFON or GEOSCOPE.From these networks we chose 9 stations (BFO, ESK, KIEV, KONO, MATE, PAB, PSZ, TAM, TIRR) at an epicentral distance of 330 − 2300 km azimuthal located around the source, see Fig. 6.The event-based raw data were extracted from the data centers and corrected for the instrument response.
Synthetic data is generated using SeisSol applied to the composite European model described above.For the simulation a source mechanism specified by the centroid moment tensor solution of the event is used, see Table 2.The data processing further involves the convolution by a source time function, shown in Fig. 5.To obtain the STF an inversion scheme using a Neighbourhood-algorithm is applied as described in St ähler et al. (2012).Subsequently the real and synthetic data are low-pass filtered to the peak frequency of 0.03Hz.

Results and discussion
For real and synthetic data the amplitude normalized traces of the vertical velocity component are shown in Fig. 6.Here we want to emphasize, that especially P-wave phases remarkably fit the measured phases at all observed stations.Surface wave arrivals usually show a larger misfit in phase and amplitude compared to real data.However, the overall signal duration matches well using the inverted source time function.The Introduction

Conclusions References
Tables Figures

Back Close
Full misfit between data and synthetics can be attributed mainly to the approximation of the material values inside the Earth by the applied velocity models.
Concluding we want to mention, that the effect of numerical boundary reflections is minor compared to the approximation of the material values.Exemplary this can be seen, at the near offset station MATE where boundary reflections only occur after the surface wave has passed.Here, similar phase and amplitude differences occur as for near boundary stations.

Conclusions
We presented a study to verify the application of the ADER-DG method for regional scale earthquake simulations.In particular, the flexibility of unstructured tetrahedral meshing allows the creation of high-quality computational grids to study wave propagation in 3-D media.
In two benchmark tests obtained using a 1-D isotropic PREM model, we quantitatively evaluated wave propagation over intermediate distances (< 15 • ) using an explosive and a shear dislocation source.In comparison to the Spectral-Element method, we could show that the ADER-DG method provides accurate synthetic results of low phase and envelope misfits.Due to the use of unstructured tetrahedral meshes numerical artifacts were obtained, which decrease in amplitude for successive mesh refinements.
In the second part the simulation of the L'Aquila earthquake confirms the applicability of the method using a regional scale 3-D composite model of Europe.The results expose an excellent fit of first arrival P-wave signals with real data measurements at a peak frequency of 0.03Hz.Surface waves, most sensitive to crustal structures, exhibit phase and amplitude misfits indicating the deficiencies in the representation of the material values inside the Earth using the EPcrust and ak135 model.
As described in Sec.2.2 features like p-adaptivity and local time stepping in addition to the flexible and fast unstructured tetrahedral mesh design demonstrated in this study can focus the computational effort using the ADER-DG method.This allows for Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | are compared.Finally, in Sect. 4 we show numerical results obtained on a composite model of Europe, where the 3-D crustal model EPcrust (Molinari and Morelli, 2011) is imposed on top of the 1-D model reference Earth model ak135(Kennett et al., 1995), and compare them to observations.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 5 • ) and variable depth.It is provided in a TomoJSON format, proposed by the NERIES project (Network of Research Infrastructures for European Seismology).A Matlab parser can read the file to a Matlab structure, which contains the material information and locations of the 3-D interfaces.The upper-mantle discontinuities of the ak135 model are manually projected on spherical shells using the same lateral sampling points.The geometrical representation of the European model is shown in the left part of Fig. 4. The pointset describing the topography and each individual interface in the model can be imported into Cubit from an ASCII file, including a list of vertex locations in Cartesian coordinates.In Cubit a surface reconstruction directly from these pointsets failed.Therefore, parallel spline curves along a row or column of the structured pointsets had to be created using the vertices as spatial support.From the generated lineset an accurate representation of the surfaces can be reconstructed within Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | source is located at (x s , y s , z s ) = (0 km, 0 km, 6361 km), which corresponds to a depth of 10 km beneath the North Pole.The source time function (STF) is given by a Ricker wavelet with a main period of T peak = 20 s.The physical domain of the SeisSol simulation is a cuboid of Ω = [−1000 km, 1000 km] × [−500 km, 3500 km] × [2400 km, 6400 km] cut by spherical interfaces at depths according to the elastodynamic discontinuities of PREM.In fact, we consider the complete 1-D isotropic PREM model with an uppermost fluid layer replaced by crustal material.The spherical sections Ω i are meshed separately.Here, the flexibility of the unstructured mesh is exploited to map material gradients within Ω i to corresponding changes in the element size.This keeps a constant number of n min = 3 tetrahedral elements per shortest wavelength in each subdomain to model a shortest wave period Discussion Paper | Discussion Paper | Discussion Paper | However, in the v t component of the SeisSol output clear waveform signals emerge, which gain up to 5 % of the peak amplitude at the v v component.This is expected to Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | As already mentioned, in SeisSol the material values of tetrahedral elements are averaged over all vertex values.Therefore, using an average element size of 18 km for mesh generation the interpolation of the material values smooths the underlying crustal model, which leads to a lowest surface wave velocity of 1.1 km s −1 .Assuming a seismic source signal at a peak frequency of 0.03 Hz the smallest wavelength of 36 km can be sampled correctly.The resulting mesh incorporates high-quality tetrahedral elements as illustrated in the enlarged sketches of Fig. 4Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | accurate simulations in regions such as the Alps where strong lateral heterogeneities in the material distribution as well as geometrical complexities due to surface and interface topography strongly perturbs 3-D wave propagation.In a future study we want to asses the quality of new Earth models and explore in detail the benefits of the local spatial and temporal adaptivity provided by the ADER-DG algorithm.Discussion Paper | Discussion Paper | Discussion Paper | aftershocks, Geophys.Res.Lett., 36, 1-6, doi:10.1029/2009GL039627,http://www.agu.org/pubs/crossref/2009/2009GL039627.shtml,2009.1146 Cupillard, P., Delavaud, E., Burgos, G., Festa, G., Vilotte, J.-P., Capdeville, Y., and Montagner, J.-P.: RegSEM: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale, Geophys.J. Int.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig. 1.Workflow of the mesh generation process using Cubit.

Fig. 2 .Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 2. Synthetic seismograms obtained on the 1-D isotropic PREM model by the ADER-DG (solid) and Spectral-Element method (dashed) compared for an explosive source located at the North Pole.The data are convolved with a Ricker source time function and low pass-filtered at 0.05Hz.Amplitudes of the radial v r and vertical v v velocity component are normalized with respect to v v .The fit between the different synthetic results is quantified by the phase (PM) and envelope (EM) time-frequency misfits.

Table 1 .
CMT solution of the 2009 Norwegian sea event.

Table 2 .
CMT solution of the 2009 L'Aquila event.