A simple diffuse interface approach on adaptive Cartesian grids for the linear elastic wave equations with complex topography

Abstract

In most classical approaches of computational geophysics for seismic wave propagation problems, complex surface topography is either accounted for by boundary-fitted unstructured meshes, or, where possible, by mapping the complex computational domain from physical space to a topologically simple domain in a reference coordinate system. However, all these conventional approaches face problems if the geometry of the problem becomes sufficiently complex. They either need a mesh generator to create unstructured boundary-fitted grids, which can become quite difficult and may require a lot of manual user interactions in order to obtain a high quality mesh, or they need the explicit computation of an appropriate mapping function from physical to reference coordinates. For sufficiently complex geometries such mappings may either not exist or their Jacobian could become close to singular. Furthermore, in both conventional approaches low quality grids will always lead to very small time steps due to the Courant-Friedrichs-Lewy (CFL) condition for explicit schemes. In this paper, we propose a completely different strategy that follows the ideas of the successful family of high resolution shock-capturing schemes, where discontinuities can actually be resolved anywhere on the grid, without having to fit them exactly. We address the problem of geometrically complex free surface boundary conditions for seismic wave propagation problems with a novel diffuse interface method (DIM) on adaptive Cartesian meshes (AMR) that consists in the introduction of a characteristic function 0≤α≤1 which identifies the location of the solid medium and the surrounding air (or vacuum) and thus implicitly defines the location of the free surface boundary. Physically, α represents the volume fraction of the solid medium present in a control volume. Our new approach completely avoids the problem of mesh generation, since all that is needed for the definition of the complex surface topography is to set a scalar color function to unity inside the regions covered by the solid and to zero outside. The governing equations are derived from ideas typically used in the mathematical description of compressible multiphase flows. An analysis of the eigenvalues of the PDE system shows that the complexity of the geometry has no influence on the admissible time step size due to the CFL condition. The model reduces to the classical linear elasticity equations inside the solid medium where the gradients of α are zero, while in the diffuse interface zone at the free surface boundary the governing PDE system becomes nonlinear. We can prove that the solution of the Riemann problem with arbitrary data and a jump in α from unity to zero yields a Godunov-state at the interface that satisfies the free-surface boundary condition exactly, i.e. the normal stress components vanish. In the general case of an interface that is not aligned with the grid and which is not infinitely thin, a systematic study on the distribution of the volume fraction function inside the interface and the sensitivity with respect to the thickness of the diffuse interface layer has been carried out. In order to reduce numerical dissipation, we use high order discontinuous Galerkin (DG) finite element schemes on adaptive AMR grids together with a second order accurate high resolution shock capturing subcell finite volume (FV) limiter in the diffuse interface region. We furthermore employ a little dissipative HLLEM Riemann solver, which is able to resolve the steady contact discontinuity associated with the volume fraction function and the spatially variable material parameters exactly. While the method is locally high order accurate in the regions without limiter, the global order of accuracy of the scheme is at most two if the limiter is activated. It is locally of order one inside the diffuse interface region, which is typical for shock-capturing schemes at shocks and contact discontinuities. We show a large set of computational results in two and three space dimensions involving complex geometries where the physical interface is not aligned with the grid or where it is is even moving. For all test cases we provide a quantitative comparison with classical approaches based on boundary-fitted unstructured meshes.