S eismic wave P ropagation and I maging in C omplex media: a E uropean network

1. Seismic wave Propagation and Imaging in Complex media: a European network • JOSEP DE LA PUENTE • Early Stage Researcher • Host Institution: LMU Munich • Place of Origin: Barcelona, Spain • Appointment Time: October 2004 • Project: Simulation of Earthquakes on Unstructured Meshes, Viscoelasticity, Anisotropy, ADER-DG Method. • Task Groups: TG Numerical Methods, TG Digital Library • Cooperation: Trento University

2. Project Scope: Unstructured Meshes in Seismology We are developing a method to simulate seismic wavefields on irregular meshes. ADER-DG (Arbitrary high order DERivatives-Discontinuous Galerkin) is born as a combination of the most convenient features of higher-order finite elements (FE) and finite volumes (FV). In the last year and a half it has been adapted to the problem of computational seismology. So far it is succesfully working in both serial and parallel, 2D or 3D with most important features included (free surfaces, kineatic rupture sources, liquid media, attenuation...). Still a couple of things remain to be implemented, such as dynamic rupturing, curved surface elements or anisotropy... ... so far not many works exist on the implementation of anisotropy for fully irregular meshes.

3. An overview of the ADER-DG Method • The Discontinuous Galerkin (DG) method has the following main features: • Triangular or tetrahedral cell decomposition. • Arbitrary high-order basis functions set. • Information is shared through „fluxes“. • Single-step high-order time integration. • Local dt /p-adaptivity. 1. Triangular or tetrahedral cell decomposition The space is discretized in a conforming mesh of triangles (or tetrahedra in 3D). This makes the method well suited for handling complex geometrical features (topography, material discontinuities...) while it is very fast and easy to set a mesh up.

4. 2. Arbitrary high-order basis functions set. We can use a basis of functions defined in a reference triangle to represent the variables inside each triangle of our mesh. With this basis we can exactly represent any polynomial of the same degree with the minimum possible number of coefficients. In addition the integral of all the functions of the basis in the triangle is zero, except for Φ0. 3. Information is shared through „fluxes“ EXAMPLE: The 3rd order basis (P2) The functions are not forced to be continuous between the cells. Spatial derivatives can be solved by using the „numerical flux“ concept (in essence, Gauss‘ Theorem). η ξ

5. 4. Single-step high-order integration in time (ADER). Accuracy of the scheme can be higher than in conventional methods due to a time-integration method which has the same order as the space-integration without need of storing intermediate substeps information (as in RK). 5. Allows for local timestepping and spatial adaptivity The polynomial order and the timestep for which we compute things at each cell can be set individually for optimal resolution and efficiency (local CFL). P1 P2P3 The ADER-DG method could be considered as a difficult to implement but easy to use method. Its very high accuracy has been already verified, whereas its computational performance is yet to be fully tested...

6. ... but why to be so accurate?? • LIFELINES Benchmark (2003), • Pacific Earthquake Engineering Research Center. • PARTICIPANTS: • UCBL: Dreger & Larsen (FD) • UCSB: Olsen (FD) • WCC2: Pitarka (FD) • CMUN: Bielak (FE) • MODEL: • Two layers, anelastic, point source (dislocation). • SOLUTION: • Frequency wavenumber quasi-analytical solution.

7. Anisotropy in ADER-DG Schemes A medium is considered anisotropic if its elastic properties depend on the direction on which waves propagate. Anisotropy might be due to materialheterogeneity (finely layered media, fluid-filled cracks...) or mineral properties (olivine, clay shales...) or a combination of both. Also stresses applied on the material can affect its anisotropic properties. The most usual observable of anisotropy is „shear wave splitting“, although in many circumstances it is difficult to distinguish it from purely geometrical effect (curved reflectors). Overall Earth‘s radial and lateral velocity variations are stronger than anisotropic effects, but some Earth regions (upper mantle, D‘‘ layer...) exhibit behaviours far from isotropy which can give important information on past and present physical processes going on in that areas. Isotropic: C11= C22= C33= λ+2μ C12= C13= C23= λ C44= C55= C66= μ

8. In the ADER-DG method we use the elastic equations in matrix form Where now A, B and C include the coefficients of the c elastic tensor. In our case, for computing fluxes through tetrahedral faces we will have to know the up to three different wave velocities at the direction normal to the face (and with its same orientation). Therefore we have to first rotate the material properties and later solve the eigenvalues of our matrix jacobian A, which are the wave velocities. αi are the eigenvalues (wave velocities)

9. A good way to qualitatively observe the anisotropic effects is by inputing an explosive source and looking at the wavefields it produces: The energy isosurfaces are not anymore spheres. Waves propagate clearly faster in the z-direction, and they show revolution symmetry around that same axis. One would also expect pure P wave propagation in the isotropic case, while now the wavefield is more complex even for such a simple setup. Along the z-axis appear the „cusps“ and we can observe coupled modes propagating in the XZ plane, where velocities are not anymore radial with respect to the source. The XY plane behaves in an analogous way as a purely isotropic medium. XZ plane XY plane

10. Main Results Theoretically our scheme is supposed to reach the desired convergence order. We make a periodic test with an orthorombic material aligned with the (1,1,1) direction in which 3 plane waves propagate simultaneously in the x, y and z directions. Rotated! Rusanov (apprx.) Godunov (exact)

11. Analytical solution comparison Test based on analytical solution first introduced in Carcione et al. 1992, homogeneous transversely isotropic material, source is a force acting on the symmetry axis of the material. Source: Ricker wavelet. (F0=20Hz, t0=0.1s.) Offset: 150m Mesh size: 20,480 el. Average dx = 25 m CFL = 0.5 Material used: Clay Shale (TI):

12. Outlook We have succesfully implemented anisotropy in the ADER-DG scheme and we have presented both qualitative and quantitative results of its accuracy. It remains a very high-order method. Different flux types have been explored. Immediate future: Running full-size analytical solution Coupling of anisotropy and viscoelasticity has to be tested. Realistic applications on a large scale still to be modelled. Far future: Comparison ADER-DG/SEM (CPU times, memory usage, accuracy). Not an easy topic! Coupled hexahedrical / tetrahedrical mesh ...and of course: