  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. In this paper we propose a completely different
strategy. We address the problem of geometrically complex free surface boundary
conditions with a novel diffuse interface method on adaptive Cartesian meshes
that consists in the introduction of a characteristic function $
0\leq\alpha\leq 1$ which identifies the location of the solid medium and the
surrounding air and thus implicitly defines the location of the free surface
boundary. 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. 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
$\alpha$ 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 $\alpha$ from
unity to zero yields a Godunov-state at the interface that satisfies the
free-surface boundary condition exactly. In order to reduce numerical
dissipation, we use high order DG finite element schemes on adaptive AMR grids
together with a high resolution shock capturing subcell finite volume (FV)
limiter in the diffuse interface region.
