Dispersion error analysis can help to assess the performance of finite-element discretizations of the wave equation. Although less general than the convergence estimates offered by standard finite-element error analysis, it can provide more detailed insight as well as practical guidelines in terms of the number of elements per wavelength needed for acceptable results. We present eigenvalue and eigenvector error estimates for cubic Hermite elements on an equidistant 1-D mesh and on a regular structured 2-D triangular mesh consisting of squares cut in half. The results show that in 1D, the spectrum consists of 2 modes. If these are unwrapped, the spectrum is effectively doubled. The eigenvalue or dispersion error stays below 7% across the entire spectrum. The error in the corresponding eigenvectors, however, increases rapidly once the number of elements per wavelength decreases to one. In terms of element size, the dispersion error is of order 6 and the eigenvector error of order 4. The latter is consistent with the classic finite-element error estimate. In 2D, we provide eigenvalue and eigenvector errors as a series expansion in the element size and obtain the same orders. 2-D numerical tests in the timeand frequency-domain are included.@en

Demand for hydrocarbon fuel is predicted to keep increasing in the coming decades in spite of easily accessible alternative fuels due to shifting geopolitical and economic situations. In order to find new hydrocarbon pockets, we need sharper images of earth’s subsurface. Also, the exploration of other sources of energy like geothermal will benefit from better models of the what lies underneath the surface. One way to obtain better images is to use superior numerical methods for forward modelling - Finite-element methods (FEM) are one such method, but their accuracy comes at the cost of increased compute expense. This thesis explores means to reduce this cost and adapt FEM to large-scale problems in geophysics. The Finite Difference (FD) method is the most popular numerical approximation scheme used in subsurface imaging problems. Representing the wave as the solution of individual motion and material equations is advantageous in terms of accuracy and stability and leads to the natural inclusion of density variations in the medium. This representation is referred to as the first-order formulation of the wave equation in this document. Finite Element (FE) methods are commonly derived for second-order equations because of the nature of variational formulation. Finite-element discretisations of the acoustic wave equation in the time domain often employ mass lumping to avoid the cost of inverting a large sparse mass matrix. Unfortunately, for a first-order system of equations, mass lumping destroys the superconvergence of numerical dispersion for odd degree polynomials. In chapter 3 of this thesis, we consider defect correction as a means to restore the convergence. We adapt the defect correction method to FEM by solving the consistent mass matrix with the lumped one as preconditioner. For the lowest-degree element, fourth-order accuracy in 1D can be obtained with just a single iteration of defect correction. In this chapter, we analyse the behaviour of the error in eigenvectors as a function of the normalized wavenumber in the form of leading terms in its series expansion and find that this error exceeds the dispersion error, except for the lowest degree where the eigenvector error is zero. We also present results of numerical experiments that confirm this analysis. Chapter 3 concluded that defect correction can improve the convergence property of finite-elements in the first-order system of acoustic equations in 1D; the inexpensive linear elements showed the same performance as a fourth-order scheme. However, for realistic problems we need to ensure that the same improvement holds in higher dimensions. Based on the results of the earlier chapter, we conjecture that defect correction should work for 2D problems. In the first half of chapter 4, we analyze the 2-D case. Theoretical results imply that the lowest-degree polynomial provides fourth-order accuracy with defect correction, if the grid of squares or triangles is highly regular and material properties constant. But numerical results converge more slowly than theoretical predictions. Further investigation demonstrates that this is due to the activation of error-inducing wavenumbers in the delta-source representation. In the second half of the chapter, we provide a solution to this problem in the form of a tapered-sinc source function. In chapter 5, we consider isotropic elastic wave propagation with continuous masslumped finite elements on tetrahedra with explicit time stepping. These elements require higher order polynomials in their interior to preserve accuracy after mass lumping and are recently discovered up to degree 4. Global assembly of the symmetric stiffness matrix is a natural approach but requires large memory. Local assembly on the fly, in the form of matrix-vector products per element at each time step, has a much smaller memory footprint. With dedicated expressions for local assembly, our code ran about 1.3 times faster for degree 2 and 1.9 times for degree 3 on a simple homogeneous test problem, using 24 cores. This is similar to the acoustic case. For a more realistic problem, the gain in efficiency was a factor 2.5 for degree 2 and 3 for degree 3. For the lowest degree, the linear element, the expressions for both the global and local assembly can be further simplified. In that case, global assembly is more efficient than local assembly. Among the three degrees, the element of degree 3 is the most efficient in terms of accuracy at a given cost. In chapter 6, we consider cubic Hermite elements as interpolants in place of Legendre polynomials. By nature of theirC1 continuity, they might offer a solution to the problems of ‘spurious’ wavenumbers seen in earlier chapters with conventional interpolation schemes. Results show acceptable convergence properties on homogeneous media, but the representation needs to be altered to suit discontinuities in density, which makes interesting future work.@en

Finite elements can, in some cases, outperform finite-difference methods for modelling wave propagation in complex geological models with topography. In the weak form of the finiteelement method, the delta function is a natural way to represent a point source. If, instead of the usual second-order form, the first-order form of the wave equation is considered, this is no longer true. Fourier analysis for a simple case shows that the spatial operator corresponding to the first-order form has short-wavelength null-vectors. Once excited, these modes are not seen by the spatial operator but only by the time- stepping scheme and show up as noise. A sourcewith a larger spatial extent, for instance a Gaussian or a tapered sinc, can avoid the excitation of problematic short wavelengths. A series of numerical experiments on a 2-D problem with an exact solution provides a suggestion for the best choice of parameters for these source term distributions. The tapered sinc provided the best results and the resulting accuracy can be better than that of the second-order form. The higher operation count of the former, however, does not make it more efficient in terms of accuracy for a given computational effort, at least not for the 2-D examples considered here@en

For seismic modelling, imaging and inversion, finite-difference methods are still the workhorse of the industry despite their inability to meet the increasing demand for improved accuracy in subsurface imaging. Finiteelement methods offer better accuracy but at a higher computational cost. A stress-velocity formulation with linear elements and an iterative method, defect correction, for inverting the mass matrix offers fourth-order super-convergence but is susceptible to numerical noise if waves in the wrong part of the dispersion curve are excited. We propose an improved source term that reduces that noise and investigate the accuracy of the method on structured triangular meshes as well as on unstructured rotated meshes. With an optimised source function, it is seen that the dispersive wavelengths can be avoided, giving the defect-correction approach a better performance than the mass-lumped formulation with only a marginal increase in compute effort. @en

The second-order formulation of the wave equation is often used for spectral-element discretizations. For some applications, however, a first-order formulation may be desirable. It can, in theory, provide much better accuracy in terms of numerical dispersion if the consistent mass matrix is used and the degree of the polynomial basis functions is odd. However, we find in the 1-D case that the eigenvector errors for elements of degree higher than one are larger for the first-order than for the second-order formulation. These errors measure the unwanted cross talk between the different eigenmodes. Since they are absent for the lowest degree, that linear element may perform better in the first-order formulation if the consistent mass matrix is inverted. The latter may be avoided by using one or two defect-correction iterations. Numerical experiments on triangles confirm the superior accuracy of the first-order formulation. However, with a delta-function point source, a large amount of numerical noise is generated. Although this can be avoided by a smoother source representation, its higher cost and the increased susceptibility to numerical noise make the second-order formulation more attractive. @en