TY - JOUR
AU - Nissen-Meyer,, Tarje
AB - SUMMARY We present a new approach to simulate high-frequency seismic wave propagation in and under the oceans. Based upon AxiSEM3D, this method supports a fluid ocean layer, with associated water-depth phases and seafloor topography (bathymetry). The computational efficiency and flexibility of this formulation means that high-frequency calculations may be carried out with relatively light computational loads. A validation of the fluid ocean implementation is shown, as is an evaluation of the oft-used ocean loading formulation, which we find breaks down at longer periods than was previously believed. An initial consideration of the effects of seafloor bathymetry on seismic wave propagation is also given, wherein we find that the surface waveforms are significantly modified in both amplitude and duration. When compared to observed data from isolated island stations in the Pacific, synthetics which include a global ocean and seafloor topography appear to more closely match the observed waveform features than synthetics generated from a model with topography on the solid surface alone. We envisage that such a method will be of use in understanding the new and exciting ocean-bottom and floating seismometer data sets now being regularly collected. Numerical modelling, Computational seismology, Guided waves, Theoretical seismology, Wave propagation 1 INTRODUCTION Oceans cover more than 70 per cent of the Earth’s surface, and have complex and nuanced effects on the propagation of seismic waves through the planet. However, with the exception of isolated stations on remote islands, global seismometer distribution is enormously skewed toward being continent-based and so biased towards the northern hemisphere. Seismology in an oceanic context has thus been somewhat neglected as compared to its land-based counterpart, and the lack of comparably high-quality global data coverage from south of the Equator can be an impediment to seismic studies. However, in recent times the deployment of ocean-bottom (OBS) and floating seismometers, together with increases in the land-based instrument density along coastlines, has begun to alleviate this problem. As a consequence there are now many large data sets which are in need of interpretation. Their complexity is such that observation alone is unlikely to prove a route to thoroughly understanding them; instead, we must make use of the synergy between observation and forward modelling. Seismologically speaking the effects of the oceans are twofold, and are particularly pronounced at high frequencies (≤10 s, where sensitivity to crustal-scale features becomes significant). First, new sources of seismic signal and noise are introduced. These include infrasounds, such as those associated with underwater volcanoes (Bohnenstiehl et al. 2013), as well as the ocean microseism, which is dynamically generated in the water column (Longuet-Higgins 1950; Ardhuin et al. 2015). Secondly, a new set of seismic phases which couple across the solid–fluid seafloor interface become supported. pWP (an upgoing P-wave reflected off the ocean surface before being retransmitted into the solid earth), and t-phases which are trapped in the SOFAR channel (Tolstoy & Ewing 1950), are examples of such. These are inherently high-frequency effects: water column reverberations like pWP have a characteristic period of a few seconds, whilst t-phases are observed at frequencies of a few Hertz and higher. In the microseism case, these are particularly useful for mapping tomographic structures (e.g. Shapiro et al. 2005; Basini et al. 2013), tracking storms (Davy et al. 2014) and making inference about climatic trends (e.g. Grob et al. 2011; Stutzmann et al. 2009). Use of the pWP phase can substantially improve source localization in Wadati-Benioff zones (Robert Engdah et al. 1998), whilst t-phases may be detected in situ in the water column, far from land and thus are useful for tomographic inversion in otherwise sparsely sampled areas (Blackman et al. 1995). Additionally, as t-phases experience a reduced degree of geometric spreading (cylindrical, |$\sim \frac{1}{r}$|⁠, rather than spherical, |$\sim \frac{1}{r^2}$|⁠), they may offer greater sensitivity to more distant or smaller seismic events than would be obtained from the corresponding P or S waves at the same epicentral distance; furthermore, the lower and better-constrained sound speed in water can yield tighter constraints on source location than are possible from signals in the solid Earth alone (Dziak et al. 2004). Other potential uses of hydroacoustic phases include early-warning systems for tsunamis (Sasorova et al. 2005; Lay et al. 2019), acoustic thermometry for remote monitoring of long-term changes in ocean temperature (Dushaw et al. 1999), detection of glacial and ice-calving events (Chapp et al. 2005), as well as monitoring of cetacean populations (e.g. McDonald et al. 1995; Dréo et al. 2019) and for illicit nuclear tests in violation of the Comprehensive Nuclear Test Ban Treaty (Mitchell 2002; Okal 2008). However, there are numerous areas of basic physics which remain to be explored within the context of ocean seismology. These include how P waves which are refracted to near-normal angles at the seafloor are converted to horizontally travelling t waves in the SOFAR channel (de Groot-Hedlin & Orcutt 2001), how primary microseismic noise with a transverse component is generated (Nishida & Takagi 2016), and how best to isolate or account for water column reverberations, especially over areas rough bathymetry (Blackman et al. 1995). The latter of these can significantly complicate understanding of earthquake dynamics (Yue et al.2017; Qian et al. 2019), especially of the largest earthquakes which occur at subduction zones with rough bathymetry (Lay & Rhode 2019; Wu et al. 2020). Such puzzles are unlikely to be resolved by observation alone, given that comprehensive seismic data from the oceans’ depths are still sparse as compared to those from on land. Making use of the interplay between observation and modelling is thus one route to exploring such questions. Multiple tools for exploring the modelling side of these problems exist, including SPECFEM3D Cartesian (Peter et al. 2011), the hybrid Direct Solution-Spectral Element Method of Wu et al. (2018) and Salvus (Afanasiev et al. 2019). However, modelling in an oceanic context is not without its own challenges, even at low frequencies. Forward modelling in seismology can be done through finding numerical solutions to the weak form of the equations of motion, which in solid media implicitly include the Neumann-type traction-free surface boundary condition (Igel 2016). In the case of a fluid surface layer, a Dirichlet boundary condition is needed and in our method must be explicitly solved. Furthermore, the low P-wave speeds in the oceans (≈1450 ms−1) yield acoustic waves with shorter wavelengths than in the underlying crust, thus necessitating the use of small elements. These in turn require shorter time steps to ensure sufficient temporal sampling in explicitly time-stepped methods, increasing the computational expense. The implementation of realistic seafloor bathymetry is also especially difficult when combined with the need to ensure that a hexahedral mesh remains conformal, with elements which are not overly deformed or skewed. These challenges have led to a number of simplifications being made in the current norm of seismic modelling, such that the simulations become computationally tractable. These have involved either approximating the ocean as a mass loading in a global formulation (where we consider ‘global’ to mean that the effects of the Earth’s sphericity are important), or introducing a realistic ocean layer but restricting the simulation to a local geometry or to fully axisymmetric domains with flat seafloors. Komatitsch & Tromp (2002) adopt the first approach, wherein the bulk modulus at the free surface is modified to account for the weight of the water column without explicitly including the ocean in the mesh. This ‘ocean loading’ formulation (also known as the ‘water column approximation’) is valid only at long periods (where the wavelength of the seismic waves is significantly larger than the ocean depth), and does not reproduce hydroacoustic phases. It does however demonstrate that the ocean delays the arrival of the Rayleigh wave train and considerably changes its dispersion characteristic, whilst leaving the Love waves unchanged. Zhou et al. (2016) present an initial evaluation of this approximation for two specific phases, as discussed further in Section 4. Conversely, Cristini & Komatitsch (2012), Bottero et al. (2016) and Mazoyer et al. (2013) use local-scale simulations in SPECFEM3D, which reproduce hydroacoustic phases but are restricted to axisymmetric ‘2.5-D’ formulations. Such simulations reproduce in-plane scattering induced by local-scale features like seamounts, but are not suitable for use at planetary scales or for simulating processes where out-of-plane scattering is thought to be important, such as Love wave generation at the seafloor. The need for global seismic modelling with realistic ocean layers and which supports their associated hydroacoustic phases at the high frequencies used in modern seismology is, therefore, clear. In this paper, we present such an implementation. Our work is based on the spectral-element methodology AxiSEM (Nissen-Meyer et al. 2014), which was expanded to include full 3-D structures (Leng et al. 2016) with undulating discontinuities and ellipticity (Leng et al. 2019). A realistic ocean layer in AxiSEM3D is now introduced, with a pressure-free surface boundary condition. A benchmark against the code YSpec (Al-Attar & Woodhouse 2008) is shown for the flat seafloor case, which introduces new water-depth seismic phases, as well as convergence to the ocean loading formulation at long periods. We show that in a global context, substantial differences between the ocean loading and realistic ocean formulations become apparent below dominant periods of ∼20 s. We also evaluate the effects of bathymetry on seismic waveforms, with a particular focus on Rayleigh waves. A clear limitation of this implementation is that it requires a consistent surface boundary condition along any line of azimuth; that is, an ocean must be either cover the planet entirely (the ‘global’ ocean) or lie in a ring along the surface (what we term a ‘local’ ocean, which is effectively doughnut-shaped). This paper deals with only the first of these cases, and though with careful choice of simulation geometry realistic simulations may still be performed, an arbitrary ‘patched’ ocean, with some solid surface areas and other fluid ones cannot be supported at present on a global scale in AxiSEM3D. Thus, this method is in no way presented as a complete replacement for synthetics generated in a solid-surface implementation, but rather as a complimentary approach which is widely applicable given the majority of our planet’s surface covered by water, and the resulting implication that an arbitrary source–receiver path at teleseismic distances is more likely than not to include an oceanic section. Such an approach is also justified on observational grounds—by undertaking a comparison to recorded data from remote island stations (which provide valuable global coverage in sparsely sampled areas and by their nature are strongly affected by the presence of an ocean), we give examples where a better match to synthetics occurs when a fluid ocean is included than when it is not. Such a methodology should enable higher-frequency simulations to be carried out in an oceanic seismology context, with arbitrary structural complexity in the underlying solid Earth. In this paper, we restrict further detailed discussion to examination of bathymetric effects and water-depth phases; however extensions of our methodology to other fluid seismology contexts (the atmosphere or localized but axisymmetric oceans) are also possible. 2 METHODOLOGY This section describes the elastodynamic theory of wave propagation in an earth model which has both a fluid surface ocean and an undulating seafloor, as well as the implementation of such a setup in AxiSEM3D. As compared to Leng et al. (2016, 2019), the main difference is in the stress-free (and hence pressure-free) boundary condition on the fluid ocean surface. 2.1 Theory We first consider a spherical earth model that consists of a solid domain ΩS (with density ρ and elasticity tensor |$\mathbf {C}$|⁠) and a fluid domain ΩF (also with density ρ and bulk modulus κ). These two domains can be separated by several solid–fluid interfaces, collectively denoted Σ, such as the ocean floor, the core–mantle boundary and the inner core boundary. In the solid domain, the weak formulation of the equations of motion (ignoring attenuation and long-period effects such as gravitation and rotation) may be written as $$\begin{eqnarray*} &&\int _{\Omega _{S}} \rho \mathbf {w} \cdot \, \partial _t^2\mathbf {u} \, d^3\mathbf {r} + \int _{\Omega _{S}} \nabla \mathbf {w}:\mathbf {C}:\nabla \mathbf {u} \, d^3\mathbf {r} \nonumber \\ &&= \int _{\Omega _{S}} \mathbf {w \cdot f} \, d^3\mathbf {r} - \int _{\Sigma } \partial _t^2 \chi \, \mathbf {\hat{n}\cdot w}d^2\mathbf {r}, \end{eqnarray*}$$(1) where |$\mathbf {\hat{n}}$| denotes the outward-pointing normal of a solid–fluid interface, |$\mathbf {f}$| is a body force source, |$\mathbf {u}$| is the displacement vector and |$\mathbf {w}$| is an arbitrary, vector-valued test function. In the fluid domain, we define a scalar potential χ such that |$\mathbf {u} = \rho ^{-1} \nabla \chi$| in ΩF, which is a descriptor for the dynamic pressure through ΩF as |$\partial _t^2\chi = -P$|⁠. $$\begin{eqnarray*} &&\int _{\Omega _{F}} \kappa ^{-1} w \, \partial _t^2\chi \, d^3\mathbf {r} + \int _{\Omega _{F}} \rho ^{-1} \nabla w \cdot \mathbf {I} \cdot \nabla \chi \, d^3\mathbf {r} \nonumber \\ &&= \int _{\Sigma } \mathbf {u}\cdot w \mathbf {\hat{n}} \ \, d^2\mathbf {r}. \end{eqnarray*}$$(2) Here, w is an arbitrary scalar-valued test function. The rank-two identity tensor |$\mathbf {I}$|⁠, which appears to be somewhat redundant here, is included to enable more natural extension to an aspherical earth model with undulating boundaries in eq. (7). Note that ΩS and ΩF are coupled by the two surface integrals over Σ. Next, we consider the stress-free boundary condition on the Earth’s surface. Without a fluid ocean, this takes the form of a Neumann boundary condition where |$\mathbf {C}:\nabla \mathbf {u}\cdot \mathbf {\hat{n}} =\mathbf {0}$|⁠, which is automatically satisfied by eq. (1). However, in the fluid ocean (where the stress-free boundary condition becomes a pressure-free one as the off-diagonal elements of the stress tensor vanish), it becomes a Dirichlet boundary condition. In this case, |$P(\mathbf {x},t) = - \partial _t^2 \chi (\mathbf {x},t) =\, 0$|⁠, which is not automatically satisfied by eq. (2) and must be explicitly prescribed. From the requirement for a steady-state solution it follows that ∂tχ and χ must also be identically zero at all times on the surface; these conditions are thus imposed at the boundary at each time step. Such boundary conditions have been used previously in a seismological context, for example by Peter et al. (2011) and Bottero et al. (2016). Such a pressure-free boundary condition is perfectly reflecting in an analytical sense, which is justified physically on the grounds that the acoustic impedance of air is very much greater than that of water so coupling from the ocean surface into the atmosphere can reasonably be neglected (though, if desired, it in can be accounted for in our formulation so long as the atmosphere is included in the mesh). At the seafloor boundary, the continuity of traction and normal velocity are implemented in the same way as at the outer core boundary. When 3-D bathymetry is incorporated, eqs (1) and (2) remain valid, but ΩS, ΩF and Σ become aspherical. As AxiSEM3D requires an axisymmetric computational domain which is spherical at a global scale, these equations cannot then be solved directly. Leng et al. (2019) describes the implementation of the particle relabelling transformation of Al-Attar & Crawford (2016) to handle such interface undulations. Given a reference spherical configuration |$\tilde{\Omega }$|⁠, and a deformed configuration Ω which is homeomorphic to |$\tilde{\Omega }$|⁠, the radial coordinate of a particle at |$\mathbf {r}$| in |$\tilde{\Omega }$| is shifted along |$\hat{\mathbf {r}}$| by an amount |$\tau (\mathbf {r})$| in |$\tilde{\Omega }$| such that $$\begin{eqnarray*} \xi (\mathbf {r}) = \mathbf {r} + \tau (\mathbf {r})\hat{\mathbf {r}}. \end{eqnarray*}$$(3) Note that this transformation only shifts the boundary inward or outward along |$\hat{\mathbf {r}}$|⁠, which is commensurate with the fact that the boundaries in question (e.g. seafloor bathymetry) have different radii at different positions, but do not need to be shifted laterally. ξ, where |$\xi :\tilde{\Omega }\rightarrow \Omega$| (and hence also |$\xi :\tilde{\Sigma }\rightarrow \Sigma$|⁠), represents the undulation mapping, which is a kinematically permissible deformation with associated deformation gradient |$\mathbf {F}$|⁠. In spherical coordinates, where |$\mathbf {F} = \mathbf {F}(r,\theta ,\phi )$|⁠, the deformation gradient can be expressed as: $$\begin{eqnarray*} \mathbf {F}\left(\tau ,\mathbf {r}\right)=\left[\nabla \left(\tau \hat{\mathbf {r}} \right)\right]^\text{T}= {\begin{pmatrix}\hat{\mathbf {r}} && \hat{\mathbf {\theta }} && \hat{\mathbf {\phi }} \end{pmatrix}} {\begin{pmatrix}\partial _r\tau && \dfrac{\partial _\theta \tau }{r} && \dfrac{\partial _\phi \tau }{r \sin \theta } \\ 0 && \dfrac{\tau }{r} && 0 \\ 0 && 0 && \dfrac{\tau }{r} \end{pmatrix}} {\begin{pmatrix}\hat{\mathbf {r}} \\ \hat{\mathbf {\theta }} \\ \hat{\mathbf {\phi }} \end{pmatrix}}, \end{eqnarray*}$$(4) whilst the associated Jacobian for this transformation may be expressed as: $$\begin{eqnarray*} \mathbf {J}(\tau ,\mathbf {r}) = \mathbf {I} + \mathbf {F}(\tau ,\mathbf {r}). \end{eqnarray*}$$(5) The weak form of the undulating, 3-D model can thus be established in the reference configuration |$\tilde{\Omega }$| in the solid [corresponding to eq. (1)] as $$\begin{eqnarray*} &&\int _{\tilde{\Omega }_{S}} \tilde{\rho }\mathbf {\tilde{w}} \cdot \, \partial _t^2\mathbf {\tilde{u}} \, d^3\mathbf {r} + \int _{\tilde{\Omega }_{S}} \nabla \mathbf {\tilde{w}}:\mathbf {\tilde{C}}:\nabla \mathbf {\tilde{u}} \, d^3\mathbf {r} \nonumber \\ &&= \int _{\tilde{\Omega }_{S}} \mathbf {\tilde{w} \cdot \tilde{f}} \, d^3\mathbf {r} - \int _{\tilde{\Sigma }} \partial _t^2\tilde{\chi }\mathbf {\tilde{n} \cdot \tilde{w}}d^2\mathbf {r}, \end{eqnarray*}$$(6) whilst in the fluid domain [corresponding to eq. (2)] the formulation becomes $$\begin{eqnarray*} &&\int _{\tilde{\Omega }_{F}} \tilde{\kappa }^{-1} \tilde{w} \, \partial _t^2\tilde{\chi }\, d^3\mathbf {r} + \int _{\tilde{\Omega }_{F}} \tilde{\rho }^{-1} \nabla \tilde{w} \cdot \mathbf {\tilde{I}} \cdot \nabla \tilde{\chi }\, d^3\mathbf {r} \nonumber \\ &&= \int _{\tilde{\Sigma }} \mathbf {\tilde{u}} \cdot \tilde{w} \mathbf {\tilde{n}} d^2 \mathbf {r}, \end{eqnarray*}$$(7) where the equivalent material parameters are given by $$\begin{eqnarray*} \tilde{\rho }(\mathbf {r}) = \rho (\xi (\mathbf {r})) |\mathbf {J}(\tau ,\mathbf {r})|, \end{eqnarray*}$$(8) $$\begin{eqnarray*} \tilde{\kappa }(\mathbf {r}) = \kappa (\xi (\mathbf {r})) |\mathbf {J}(\tau ,\mathbf {r})|, \end{eqnarray*}$$(9) $$\begin{eqnarray*} \mathbf {\tilde{C}} (\mathbf {r}) = \mathbf {J}^{-1}(\tau ,\mathbf {r}) \cdot \mathbf {C}(\xi (\mathbf {r})) \cdot \mathbf {J}^T(\tau ,\mathbf {r})|\mathbf {J}(\tau ,\mathbf {r})|, \end{eqnarray*}$$(10) and the equivalent source by $$\begin{eqnarray*} \tilde{f}(\mathbf {r}) = f(\xi (\mathbf {r})) |\mathbf {J}(\tau ,\mathbf {r})|, \end{eqnarray*}$$(11) and |$\mathbf {\tilde{I}}$| and |$\mathbf {\tilde{n}}$|⁠, respectively, by $$\begin{eqnarray*} \mathbf {\tilde{I}} = \mathbf {J}^{-1}(\tau ,\mathbf {r}) \, \cdot \, \mathbf {J}^{-T}(\tau ,\mathbf {r}), \end{eqnarray*}$$(12) $$\begin{eqnarray*} \mathbf {\tilde{n} (r)} = \mathbf {J}^{-T}(\tau ,\mathbf {r}) \cdot \mathbf {\hat{n}}(\xi (\mathbf {r})) |\mathbf {J}(\tau ,\mathbf {r})|. \end{eqnarray*}$$(13) The solutions to eqs (6) and (7), in |$\mathbf {u}$| and χ, respectively, are related to the solutions of eqs (1) and (2) by $$\begin{eqnarray*} \mathbf {\tilde{u}} (\mathbf {r},t) = \mathbf {u}(\xi (\mathbf {r},t)), \end{eqnarray*}$$(14) and $$\begin{eqnarray*} \mathbf {\tilde{\chi }} (\mathbf {r},t) = \mathbf {\chi }(\xi (\mathbf {r},t)). \end{eqnarray*}$$(15) This formulation is sufficient to describe all variables needed to undertake the task of finding solutions to the elastodynamic equations of motion in an aspherical Earth which possesses undulating boundaries. 2.2 Implementation in AxiSEM3D Eqs (1) and (2) may be solved in a fully 3-D mesh, as is the case in SPECFEM (Komatitsch & Tromp 2002). AxiSEM3D, however, makes use of the axisymmetric formulation of Nissen-Meyer et al. (2007). The smoothness of most global tomographic models (which have significantly smaller gradients in seismic parameters in the two lateral directions than along their radius) means that the global seismic wavefield is, in general, also significantly smoother in the azimuthal direction than it is in a meridional plane. As such the use of a pseudospectral parametrization in the azimuthal direction which exploits this wavefield smoothness and accounts for the periodicity of the solution over the interval (0, 2π) can offer a substantial computational speed-up as compared to methods like SPECFEM3D Globe which rely on a full ‘cubed-sphere’ mesh. This speedup scales with increasing frequency—at 10 s, AxiSEM3D is approximately two orders of magnitude faster (Leng et al. 2019). The choice of the highest term which must be included in the azimuthal Fourier expansion (Nu) depends on the complexity of the model. In a radially symmetric (1-D) model, this representation becomes analytic for a second-order (quadropolar) moment tensor when Nu = 2. For more complex models with significant 3-D structures, Nu can be increased as needed to capture azimuthal structures. Working in a cylindrical basis (⁠|$\mathbf {\hat{s}},\mathbf {\hat{\phi }},\mathbf {\hat{z}} = \mathbf {\hat{g}_{1}},\mathbf {\hat{g}_{2}},\mathbf {\hat{g}_{3}}$|⁠), the velocity |$\mathbf {u}$| and potential χ may be expressed as: $$\begin{eqnarray*} \mathbf {u} = \mathbf {\hat{g}_i}(\phi )\, u_i(s,\phi ,z;t) = \sum _{\mid {\alpha }\mid \le N_{u}} u_i^\alpha (s,z;t)\Psi ^\alpha (\phi )\mathbf {\hat{g}_i}(\phi ) \end{eqnarray*}$$(16) and $$\begin{eqnarray*} \chi = \chi _i(s,\phi ,z;t) = \sum _{\mid {\beta }\mid \le N_{u}} \chi _i^\beta (s,z;t)\Psi ^\beta (\phi ), \end{eqnarray*}$$(17) respectively, where |$\Psi ^\alpha (\phi ) = e^{\sqrt{-1}\alpha \phi }$| are the terms with order α ≤ |Nu|, and similarly for Ψβ(ϕ). As per eqs (16) and (17), the decomposition of the solution into Fourier modes means that a 2-D mesh may be used (see Nissen-Meyer et al. 2007; Leng et al. 2016). 3-D structural complexity may be added as required, subject to an increase in the Fourier order Nu, whilst topography may be implemented on structural discontinuities (e.g. the seafloor and Moho) through the particle relabelling technique described in Section 2.1. Thus, full 3-D structures may be reliably accounted for in simulations involving AxiSEM3D. An example of an explicitly meshed surface fluid layer is shown in Fig. 1. Note that we use the terms ‘bathymetry’ and ‘seafloor topography’ interchangeably to describe the variation in ocean depth with position. Figure 1. Open in new tabDownload slide A 2-D mesh as used in AxiSEM3D, with the colour scale corresponding to the value of Vs. Panel (a) shows the entire mesh in the ocean case, from surface to core, whilst panel (b) shows a detail of the surface layers in the no ocean case and panel (c) is the equivalent detail for the ocean containing mesh. The surface is leftmost in all cases. The differences between panels (b) and (c) are slight at this scale: the ocean mesh has a thin blue layer at its outer edge (where Vs = 0) and an extra mesh refinement layer midway through the displayed section to account for the velocity contrast at the seafloor, whilst the mesh without an ocean has neither. Note that in all regions of the mesh, the smallest non-zero velocity (i.e. Vs in the solid and Vp in the fluid) is the dominant determinant of the element size. Figure 1. Open in new tabDownload slide A 2-D mesh as used in AxiSEM3D, with the colour scale corresponding to the value of Vs. Panel (a) shows the entire mesh in the ocean case, from surface to core, whilst panel (b) shows a detail of the surface layers in the no ocean case and panel (c) is the equivalent detail for the ocean containing mesh. The surface is leftmost in all cases. The differences between panels (b) and (c) are slight at this scale: the ocean mesh has a thin blue layer at its outer edge (where Vs = 0) and an extra mesh refinement layer midway through the displayed section to account for the velocity contrast at the seafloor, whilst the mesh without an ocean has neither. Note that in all regions of the mesh, the smallest non-zero velocity (i.e. Vs in the solid and Vp in the fluid) is the dominant determinant of the element size. AxiSEM3D supports arbitrary modifications of the seismic profile used in meshing (as generated using the salvus_mesher_lite package, detailed in Afanasiev et al. (2019)). This enables inclusion of heterogeneities, anisotropy (van Driel & Nissen-Meyer 2014a) and attenuation (van Driel & Nissen-Meyer 2014b). In this paper, all simulations were performed without ellipticity, rotation or gravitation as we use only intermediate periods (up to 100 s), and with attenuation enabled in the solid earth only due to the comparatively high Q value of acoustic wave propagation in water. Simulations are performed on an anisotropic PREM mesh (6368 km radius without ocean, and 6371 km with ocean). The ocean is modelled as a 3-km-thick homogeneous layer with vp = 1450 ms−1 and ρ = 1040 kgm−3, though this may be arbitrarily modified if required. Two elements per wavelength are used in all AxiSEM3D simulations. 3 VALIDATION 3.1 Nomenclature For clarity, here we will briefly summarize the terms used to describe the different ocean configurations used in this paper. In Sections 3 and 4, all models are spherically symmetric (i.e. use 1-D seismic profiles only). ‘No Ocean’ refers to a PREM model of radius 6368 km with a solid surface, ‘Ocean Load’ is identical but accounts for the weight of a 3-km-thick ocean layer without explicitly meshing it, whilst ‘Fluid Ocean’ indicates a PREM model with a global, explicitly meshed 3-km-thick homogeneous ocean layer and outer radius 6371 km. From Section 5 onward, non-spherically symmetric (fully 3-D) models are used. The exact details are specified in Section 5, but in short, ‘No Ocean’ indicates a 6368 km radius model with undulation on the solid surface (surface topography), ‘Fluid Ocean without Bathymetry’ indicates a model with global ocean, radius 6371 km and no undulation on the seafloor (i.e. no bathymetry), and ‘Fluid Ocean with Bathymetry’ is otherwise identical except for the inclusion of seafloor undulation. As it is not the main topic of this paper, the ocean loading approximation is not considered beyond the end of Section 4. 3.2 Benchmark parameters In order to verify the reliability of the AxiSEM3D simulations, we perform a benchmark against YSpec (Al-Attar & Woodhouse 2008). YSpec is a semi-analytical method which uses direct radial integration to compute full waveform synthetics, chosen specifically because it is not a spectral element method, and hence provides a robust benchmark. Unlike AxiSEM3D, YSpec does support calculations with full self-gravity, but these are not used here. For reproducibility, it should be noted that we edited the YSpec source code to remove an automatic ellipticity correction. A radially symmetric model with an ocean layer is chosen as this is the most complex configuration for which a reference solution with a global ocean is available. It should be noted that we are not therefore undertaking a full, global benchmark with 3-D structures [which in the fluid ocean case is not possible with open-source software at high frequency, and in the no ocean case was done by Leng et al. (2016)]. Rather, we seek to confirm that the addition of a fluid ocean layer does not detract from the reliability of the synthetics. The source parameters used are given in Table 1. In YSpec the maximum angular degree used is 22 000 and the maximum frequency is 600 mHz, whilst in AxiSEM3D the dominant simulation period is 2 s. Note we define the ‘dominant period’ as the minimum globally resolved period in the mesh, whilst Komatitsch & Tromp (2002) define it as the corner frequency above which no energy is observed in the simulation. The stations in the benchmark are located at the surface in the no ocean case, and on the seafloor in the ocean case (i.e. at radius of 6368 km in both). In AxiSEM3D, the order of the Gauss–Lobato–Legendre quadrature is 4. Table 1. Source parameters. To enable exact reproducibility of our benchmark, we give the full parameters of the source used here. Latitude 90° Longitude 0 ° Depth 12 km (no ocean), 15 km (ocean) Receiver position Surface (no ocean), Seafloor (ocean) Source-time function Heaviside (YSpec), Error (AxiSEM3D) Moment tensor Mrr = 1 x 1026 dyne-cm Latitude 90° Longitude 0 ° Depth 12 km (no ocean), 15 km (ocean) Receiver position Surface (no ocean), Seafloor (ocean) Source-time function Heaviside (YSpec), Error (AxiSEM3D) Moment tensor Mrr = 1 x 1026 dyne-cm Open in new tab Table 1. Source parameters. To enable exact reproducibility of our benchmark, we give the full parameters of the source used here. Latitude 90° Longitude 0 ° Depth 12 km (no ocean), 15 km (ocean) Receiver position Surface (no ocean), Seafloor (ocean) Source-time function Heaviside (YSpec), Error (AxiSEM3D) Moment tensor Mrr = 1 x 1026 dyne-cm Latitude 90° Longitude 0 ° Depth 12 km (no ocean), 15 km (ocean) Receiver position Surface (no ocean), Seafloor (ocean) Source-time function Heaviside (YSpec), Error (AxiSEM3D) Moment tensor Mrr = 1 x 1026 dyne-cm Open in new tab The source is represented through a near-instantaneous impulse, with a source time function which is a Heaviside in YSpec and a narrow error function (half-width 0.5 s) in AxiSEM3D. The result is that a Green’s function is extracted in both cases, with the output Butterworth-Bandpass filtered (filter order 4) between 2 and 100 s. The remaining differences between these Heaviside and error source-time functions are accounted for by a small (∼1 s) temporal shift and convolution of both traces with a 2 s half-width Gaussian function. 3.3 Benchmark results The modifications induced by the addition of the ocean layer are apparent in Fig. 2, whilst Fig. 3 shows the time-distance record section of the observed waveforms in the fluid ocean case. The Rayleigh wave train is substantially extended, lasting many times its original length, whilst the peak surface wave amplitude is decreased and the dispersion characteristic is modified. The changes to the body waveforms are more subtle, but so-called ‘dog-leg multiples’ (reverberations in the ocean column) become visible, and may be identified by their characteristic period of ∼5 s which corresponds to the vertical path length within the water column. Identical behaviour is seen in the radial component of the seismograms, whilst the transverse component remains unchanged in the presence of a homogenous ocean with uniform seafloor. Figure 2. Open in new tabDownload slide The modification to the waveforms from the addition of the ocean at 2 s resolution. The top row shows the seismograms in a PREM model with no ocean, whilst the bottom shows seismograms with a 3-km-deep ocean layer in an otherwise identical setting. The vertical component is shown in all cases. The whole seismogram (including dispersive surface waves) is shown in the left-hand column at 10° from the source; whilst the right column shows the modification to the P-wave train at 40° associated with the trapping of water depth phases in the ocean. For reference, the TauP arrival times of a 3 km s–1 surface wave and P phase are shown, respectively in the left- and right-hand panels. The lowermost row shows the residuals (uaxisem − uyspec), and it should be noted that the scales in the bottom two panels are different to those above. Figure 2. Open in new tabDownload slide The modification to the waveforms from the addition of the ocean at 2 s resolution. The top row shows the seismograms in a PREM model with no ocean, whilst the bottom shows seismograms with a 3-km-deep ocean layer in an otherwise identical setting. The vertical component is shown in all cases. The whole seismogram (including dispersive surface waves) is shown in the left-hand column at 10° from the source; whilst the right column shows the modification to the P-wave train at 40° associated with the trapping of water depth phases in the ocean. For reference, the TauP arrival times of a 3 km s–1 surface wave and P phase are shown, respectively in the left- and right-hand panels. The lowermost row shows the residuals (uaxisem − uyspec), and it should be noted that the scales in the bottom two panels are different to those above. Figure 3. Open in new tabDownload slide Benchmark results for the fluid ocean case between AxiSEM3D and YSpec, bandpass filtered between 2 and 100 s in the vertical component. The red trace shows the reference solution (YSpec), whilst the black trace shows the AxiSEM3D synthetics. Figure 3. Open in new tabDownload slide Benchmark results for the fluid ocean case between AxiSEM3D and YSpec, bandpass filtered between 2 and 100 s in the vertical component. The red trace shows the reference solution (YSpec), whilst the black trace shows the AxiSEM3D synthetics. Strong agreement between AxiSEM3D and YSpec is demonstrated in both cases. To quantify this, we consider the time-frequency misfits between the two seismograms (after Kristeková et al. (2009)) using the ObsPy package (Beyreuther et al. 2010). The phase misfit remains increases slightly, from 0.02 to 0.03, when the ocean layer is added, whilst the envelope misfit remains unchanged at 0.01. According to the Kristekova ćlassification, both the envelope and phase misfits are ‘excellent’. Thus, we conclude that the implementation of the fluid surface layer in AxiSEM3D has been verified. For completeness, we also demonstrate that the effects of the oceans become negligible at long seismic periods, justifying their lack of inclusion in low-frequency seismic simulations. Fig. 4 shows the convergence of the fully fluid ocean synthetics to those from the PREM case without ocean (and without ocean loading) at 50 s dominant period. Apart from the differences in the structural model, the simulation parameters are otherwise identical to those given in Table 1. Figure 4. Open in new tabDownload slide Comparison between the fluid ocean formulation and the PREM case without ocean, both simulated in AxiSEM3D, in the vertical component at 20° from the source. Synthetics are bandpass filtered between 50 and 100 s. No additional convolution with a Gaussian function is needed in this case as both source-time functions are identical. Figure 4. Open in new tabDownload slide Comparison between the fluid ocean formulation and the PREM case without ocean, both simulated in AxiSEM3D, in the vertical component at 20° from the source. Synthetics are bandpass filtered between 50 and 100 s. No additional convolution with a Gaussian function is needed in this case as both source-time functions are identical. 4 EVALUATION OF OCEAN LOADING Having established the accuracy of AxiSEM3D simulations with a fluid ocean layer through means of a benchmark, we now consider the reliability of the ocean-loading formulation (Komatitsch & Tromp 2002). This is done by comparing AxiSEM3D seismograms in a PREM model overlain by a realistic ocean to those in PREM with an ocean load at the surface. The ocean loading formulation itself was validated in AxiSEM3D against SPECFEM by Leng et al. (2016). This is done for a variety of ocean depths and seismic periods at 20° from the source, which is distant enough to capture the modified dispersion but close enough that the surface-wave train is well-recorded in the simulation interval. Such an approach is similar in aim to that used by Zhou et al. (2016), but with two key differences. First, we use a global, spherical Earth with realistic velocity and density structure, rather than a homogeneous half-space. Secondly, we also consider the effects of the ocean on the entire wave train, whilst their studies are restricted to the PP phase and Rayleigh waves. Fig. 5 shows the evolution of the waveforms with increasing distance from the source for an ocean depth of 3 km in the vertical component, whilst Fig. 6 shows the reducing discrepancy between two formulations at increasing seismic periods. At greater than 30 s dominant period (lower right), no difference is resolvable visually. At less than 20 s, the extension to the Rayleigh wave coda is clear. By 10 s, the peak amplitudes are substantially overpredicted in the ocean loading formulation, whilst the coda length is much curtailed. At less than 7.5 s minimum period, differences in the body wave arrivals become apparent, whilst even the gross structures in the surface waveforms are missed, with only arrival times being predicted correctly. Animated versions of these comparisons, showing the transition in waveforms with decreasing period, are presented in the Supporting Information. Figure 5. Open in new tabDownload slide Comparison between the ocean loading and fluid ocean formulations for a 3 km deep ocean layer in the vertical component, both in AxiSEM3D. Seismograms are bandpass filtered between 2 and 100 s. Figure 5. Open in new tabDownload slide Comparison between the ocean loading and fluid ocean formulations for a 3 km deep ocean layer in the vertical component, both in AxiSEM3D. Seismograms are bandpass filtered between 2 and 100 s. Figure 6. Open in new tabDownload slide Comparison between the fluid ocean and ocean loading formulations for a 3 km depth ocean, in the vertical component, at a variety of different minimum periods 20° from the source. Note that the vertical scale is different for each row of figures. Figure 6. Open in new tabDownload slide Comparison between the fluid ocean and ocean loading formulations for a 3 km depth ocean, in the vertical component, at a variety of different minimum periods 20° from the source. Note that the vertical scale is different for each row of figures. Thus, we conclude that the ocean loading formulation works as expected above 20 s period for the entire wave train, but becomes unreliable at higher temporal resolutions. This result is consistent with the findings of Zhou et al. (2016), who identify modification to the Rayleigh wave dispersion characteristic and changes in the amplitude and arrival times of the body waves as being key features not predicted by the loading formulation. For reference, neither formulation has any effect on the transverse component of the waveforms, confirming that in a radially symmetric model with an isotropic ocean, the ocean’s effects are confined to the source–receiver plane. Of course, the ocean depth is not a uniform 3 km across the Earth’s surface. For this reason, we also consider the envelope and phase misfits between the two implementations for a variety of ocean depths (Fig. 7). Figure 7. Open in new tabDownload slide Phase misfits (top panel) and envelope misfits (bottom panel) between the ocean loading and realistic ocean implementations, with the vertical component sampled at 20° epicentral distance. Figure 7. Open in new tabDownload slide Phase misfits (top panel) and envelope misfits (bottom panel) between the ocean loading and realistic ocean implementations, with the vertical component sampled at 20° epicentral distance. A greater misfit is observed for the deeper oceans, as is expected where the ocean depth becomes more comparable to the seismic wavelength. A steep ‘shoulder’ is also clear in the intermediate period range, wherein the phase and envelope misfits both decrease rapidly with increasing dominant period, suggesting strong convergence to the ocean loading formulation at longer periods. In general, it appears that the degree of misfit is controlled by the ratio of minimum resolved seismic wavelength to ocean depth, though as the misfit calculation is non-linear and unreliable for large differences, no simple qualitative relationship for the transition is apparent. The degree of misfit also increases with increasing epicentral distance from the source, which we attribute to the increased length of the surface wave train at greater distances and the larger number of propagated wavelengths for in-ocean phases like pWP. 5 BATHYMETRY IN A GLOBAL OCEAN CONTEXT 5.1 Bathymetric implementation Having considered the effects of the addition of a constant-thickness fluid ocean layer on synthetic seismograms, we now consider the effects of undulating the seafloor solid–fluid boundary. In our current formulation, the ocean must cover the entire planet to accommodate homeomorphic boundary undulation, and ocean depths <500 m (continents and their shelves) are linearly scaled to lie between 1500 and 500 m below sea level. This ensures that the time step does not become unfeasibly small in thin ocean layers overlying land areas. Naturally, this restricts the domain of applicability to cover ocean basins and purely oceanic ray paths. The Earth’s current continental configuration is such that many oft-used paths exist, for example at regional-to-teleseismic distances across the Pacific or Atlantic Oceans. The ocean density and sound speed are constant throughout this model, but can be arbitrarily made to vary by latitude, longitude, and depth to reproduce geographical variations in seismic parameters caused by variations in water temperature, pressure and salinity. Further work in this area will allow for global-scale meshes which support, for example, the propagation of high-frequency t-phases through the SOFAR channel. To illustrate the effects of bathymetry we consider a fluid ocean with and without bathymetry. Simulations are conducted at 5 s dominant period, apart from the animations which are at 10 s due to the enormous memory needed to produce them. All simulations use Crust 1.0 Moho (Laske et al. 2013), ETOPO1 bathymetry [sampled at 1° resolution, Amante & Eakins (2009)], and the SEMUCB-WM1 volumetric tomographic model of French & Romanowicz (2014) (sampled at 0.5° intervals and with δvp = 0.5δvs). The source used is the Mw 6.6 earthquake which occurred in New Britain, Papua New Guinea, on 2015 April 30. Arrival times for TauP (Crotwell et al. 1999) in PREM, that is without the fluid ocean, are also shown, with a small shift to account for the implementation of the source–time function in AxiSEM3D. To capture the complexity of the 3-D models used, Nu is set to 1500 in the uppermost 100 km of the crust and mantle, and 200 below 400 km depth (depths are relative to 6371 km in all cases to ensure realism and consistency between the ocean and no ocean models). Between these depths linear interpolation of Nu is used. It should be noted that if smoother structural models are used, lower Nu may be used, decreasing computational cost significantly. In undertaking particle relabelling on multiple boundaries (e.g. both the Moho and the seafloor), care must be taken to ensure that elements do not become too skewed (overly large bathymetric gradients being accommodated by too few elements), or too small (the Moho and the seafloor moving towards each other such that the element shrinks and the time step becomes too small). With the coarse bathymetric sampling used here this did not pose an issue, but in simulations using stronger bathymetry this may be an issue. At higher frequencies, with more elements across which to spread a given gradient, this may be less of a concern. Green’s Functions are generated for island stations in the Pacific Ocean, and convolved with the earthquake average source–time function from the SCARDEC database (Vallée & Douet 2016). Note that whilst the stations are located on land, we relocate them to the seafloor in our global ocean formulation. Thus, any receiver-side conversions on the islands’ submarine slopes are neglected; however these are not expected to be significant at these seismic periods given the small size of these islands (all have an above-water extent that is much less than 1° in any direction and hence they are not properly sampled anyhow our implemented bathymetric model). Source parameters are given in Table 2 and receiver locations used in this study in Table 3, with seafloor-projected ray paths shown in Fig. 8. Figure 8. Open in new tabDownload slide Source–receiver paths from the 2015 April 30 Papua New Guinea earthquake, plot generated in Google Earth Pro using Landsat/Copernicus/IBCAO images, SIO/NOAA/US Navy/NGA/GEBCO data). Figure 8. Open in new tabDownload slide Source–receiver paths from the 2015 April 30 Papua New Guinea earthquake, plot generated in Google Earth Pro using Landsat/Copernicus/IBCAO images, SIO/NOAA/US Navy/NGA/GEBCO data). Table 2. Source parameters. Date . 2015 April 30 . Time 10:45 UTC Latitude 5.6°S Longitude 151.7 °W Depth 38.3 km (below mean sea level) Mw 6.6 Half-duration 11.3 s Date . 2015 April 30 . Time 10:45 UTC Latitude 5.6°S Longitude 151.7 °W Depth 38.3 km (below mean sea level) Mw 6.6 Half-duration 11.3 s Open in new tab Table 2. Source parameters. Date . 2015 April 30 . Time 10:45 UTC Latitude 5.6°S Longitude 151.7 °W Depth 38.3 km (below mean sea level) Mw 6.6 Half-duration 11.3 s Date . 2015 April 30 . Time 10:45 UTC Latitude 5.6°S Longitude 151.7 °W Depth 38.3 km (below mean sea level) Mw 6.6 Half-duration 11.3 s Open in new tab Table 3. Receiver locations. Name . Location . Latitude . Longitude . Distance . II.KWAJ Kwajalein Atoll, Marshall Islands 8.8°N 167.6°E 21.6° IU.WAKE Wake Island, US Minor Outlying Islands 19.3°N 166.7°E 28.9° IU.RAO Kermadec Islands, New Zealand 29.4°S 177.9°W 37.3° IU.JOHN Johnston Atoll, US Minor Outlying Islands 16.7°N 169.5°W 44.6° Name . Location . Latitude . Longitude . Distance . II.KWAJ Kwajalein Atoll, Marshall Islands 8.8°N 167.6°E 21.6° IU.WAKE Wake Island, US Minor Outlying Islands 19.3°N 166.7°E 28.9° IU.RAO Kermadec Islands, New Zealand 29.4°S 177.9°W 37.3° IU.JOHN Johnston Atoll, US Minor Outlying Islands 16.7°N 169.5°W 44.6° Open in new tab Table 3. Receiver locations. Name . Location . Latitude . Longitude . Distance . II.KWAJ Kwajalein Atoll, Marshall Islands 8.8°N 167.6°E 21.6° IU.WAKE Wake Island, US Minor Outlying Islands 19.3°N 166.7°E 28.9° IU.RAO Kermadec Islands, New Zealand 29.4°S 177.9°W 37.3° IU.JOHN Johnston Atoll, US Minor Outlying Islands 16.7°N 169.5°W 44.6° Name . Location . Latitude . Longitude . Distance . II.KWAJ Kwajalein Atoll, Marshall Islands 8.8°N 167.6°E 21.6° IU.WAKE Wake Island, US Minor Outlying Islands 19.3°N 166.7°E 28.9° IU.RAO Kermadec Islands, New Zealand 29.4°S 177.9°W 37.3° IU.JOHN Johnston Atoll, US Minor Outlying Islands 16.7°N 169.5°W 44.6° Open in new tab 5.2 The effects of bathymetry on waveforms We begin by making a high-resolution comparison of the body waveforms in simulations with and without bathymetry. Such synthetics are presented in Fig. 9. Figure 9. Open in new tabDownload slide Synthetic seismograms for station IU.JOHN on Johnston Atoll, Hawaii. The upper panel compares the full, bathymetric ocean with the no ocean (PREM and Crust 1.0) case, whilst the lower panel shows the comparison between the bathymetry and no bathymetry cases. Green’s functions are convolved with the SCARDEC source–time function and bandpass filtered at 5 s. The vertical component is shown and TauP times for the PREM case are overlaid. Figure 9. Open in new tabDownload slide Synthetic seismograms for station IU.JOHN on Johnston Atoll, Hawaii. The upper panel compares the full, bathymetric ocean with the no ocean (PREM and Crust 1.0) case, whilst the lower panel shows the comparison between the bathymetry and no bathymetry cases. Green’s functions are convolved with the SCARDEC source–time function and bandpass filtered at 5 s. The vertical component is shown and TauP times for the PREM case are overlaid. Comparison between the top and bottom panels reinforces the finding that the effects of bathymetry are less appreciable than those which come about from the addition of a fluid ocean layer at a resolution of 5 s, but nonetheless are notable. Examination of the bottom panel reveals that the coda associated with the ringing of the body waves in the ocean column is modified by the roughening of the seafloor, and the modifications to PP (which bounces off the seafloor) are more substantial than those to P, which does not. The most substantial changes are in the surface waves, where the ratio of peak-to-mean amplitude is reduced, and a much more substantial coda present. These effects are likely more pronounced at higher frequencies, as the shorter-wavelength surface waves are more sensitive to the small-scale bathymetric structures present along the seafloor. Of course, such effects are strongly path-dependent. Fig. 10 illustrates this by considering the effects of the ocean and bathymetry on body waves at four separate stations. As expected, the most significant differences are observed in at stations where the source–receiver path has strong bathymetric gradients (e.g. IU.RAO), whilst stations along reasonably flat source–receiver paths (e.g. IU.JOHN) do not. Figure 10. Open in new tabDownload slide Detail of body wave arrivals in the vertical component, bandpass filtered between 5 and 100 s. TauP traveltimes for the phases P and PP in the PREM case are overlaid. Traces are normalized to the peak P-wave amplitude. Figure 10. Open in new tabDownload slide Detail of body wave arrivals in the vertical component, bandpass filtered between 5 and 100 s. TauP traveltimes for the phases P and PP in the PREM case are overlaid. Traces are normalized to the peak P-wave amplitude. We now present an animated comparison of the differences between the wavefields recorded at the seafloor in the case of the global fluid ocean with bathymetry and without bathymetry; with a focus on the surface waves. These visualizations are presented in Fig. 11. Note that in these simulations, conducted at 10 s dominant period, an 11 s half-duration Gaussian source–time function is used to ensure smooth interpolation and visualization of the wavefield. The larger amplitude of the surface waves (leading to larger differences), and their greater sensitivity to near-surface conditions mean that they dominate the difference plots. Figure 11. Open in new tabDownload slide A visualization of the simulated wavefield from the 2015 April 30 Papua New Guinea earthquake, as observed on the seafloor. The upper row of figures shows synthetics with a fluid ocean but without bathymetry, whilst the middle row is otherwise identical apart from the inclusion of bathymetry. The bathymetric colourbar applies to all three rows, and for ease of reference displays the physical elevation in ETOPO1 rather than the scaled elevation which we implement (for a detailed comparison to the scaling which we apply, see the Supporting Information). The L2 norm of the two displacement vectors (||uwb|| and ||unb|| , respectively) is shown in the top two rows, where red represents the largest displacements and white the smallest. The data range over which the colour is scaled in both is [0, 2×10−5 m] to ensure that the plots are not dominated by the surface wave amplitudes. For clarity, this means that displacements in the range [2×10−5 m,∼4×10−5 m] all appear as the the deepest shade of red, where the latter value is the maximum amplitude observed in the simulations. The lowermost row is the L2 norm of the difference in displacements between the bathymetric and non-bathymetric simulations, that is ||uwb − unb||. In this row, the colour scale saturates at twice the maximum amplitude (4×10−5 m). An additional linear ‘suppression’ filter is applied to the wavefield opacities in the range [0, ∼1×10−5 m]. This has the effect of hiding small differences where ||uwb − unb|| ∼ 0 (i.e. making them ‘transparent’) and ensures that the non-causal differences in numerical noise between the two simulations do not obscure the entire globe [will be a full landscape page]. Figure 11. Open in new tabDownload slide A visualization of the simulated wavefield from the 2015 April 30 Papua New Guinea earthquake, as observed on the seafloor. The upper row of figures shows synthetics with a fluid ocean but without bathymetry, whilst the middle row is otherwise identical apart from the inclusion of bathymetry. The bathymetric colourbar applies to all three rows, and for ease of reference displays the physical elevation in ETOPO1 rather than the scaled elevation which we implement (for a detailed comparison to the scaling which we apply, see the Supporting Information). The L2 norm of the two displacement vectors (||uwb|| and ||unb|| , respectively) is shown in the top two rows, where red represents the largest displacements and white the smallest. The data range over which the colour is scaled in both is [0, 2×10−5 m] to ensure that the plots are not dominated by the surface wave amplitudes. For clarity, this means that displacements in the range [2×10−5 m,∼4×10−5 m] all appear as the the deepest shade of red, where the latter value is the maximum amplitude observed in the simulations. The lowermost row is the L2 norm of the difference in displacements between the bathymetric and non-bathymetric simulations, that is ||uwb − unb||. In this row, the colour scale saturates at twice the maximum amplitude (4×10−5 m). An additional linear ‘suppression’ filter is applied to the wavefield opacities in the range [0, ∼1×10−5 m]. This has the effect of hiding small differences where ||uwb − unb|| ∼ 0 (i.e. making them ‘transparent’) and ensures that the non-causal differences in numerical noise between the two simulations do not obscure the entire globe [will be a full landscape page]. Video sequences were generated using ImageJ (Abràmoff et al. 2004). We include only four pertinent snapshots of the wavefield in this paper, however the full animations may be found on the ‘Oxford Seismology’ YouTube channel: https://www.youtube.com/c/SeismologyOxford. In interpreting the differences in Fig. 11, it should be noted that whilst the simulations are consistent and accurate throughout, the wavefield observed in the sector clockwise from northwest to southeast is most realistic than that toward the southwest, where the scaling of the Australian continent to lie underwater occurs. The difference wavefield is not physically meaningful, but by tracking its ‘propagation’ the regions of bathymetry which significantly influence seismic wave propagation may be examined. One hundred seconds after the event, two regions of significant difference are apparent. The northern portion corresponds to the passage of the surface waves over the island of New Britain and into the shallow Bismarck Sea, whist the southern portion emerges as the wave fronts pass over the steep topography of the New Britain Trench and into the deeper Bismarck Sea. By 300 s, these differences have become more pronounced as elastically scattered remnants of the surface wave train appear to remain trapped in the bathymetric simulation for a substantial duration. The absence of these ‘fingers’ of remnant wavefield in the no bathymetry plots indicates that they are associated with the seafloor topography in these regions, rather than the Moho or crustal structure. At 500 s these fingers remain pronounced, with the strongest appearing around the regions of rough bathymetry surrounding the Caroline Seamounts in the Federated States of Micronesia. Other substantial regions of difference appear in the Coral Sea off the eastern coast of Australia, which may be due to wavefield trapping in the ‘valley’ between the Australian coast and the shallower seas of the Lord Howe Rise, in the Celebes Sea, and to a lesser degree in the Mariana Trench region of the Philippine Sea. By 1100 s, the largest differences are in the surface waves over the Australian continent, which are not realistically modelled and hence we do not consider their physical origin. Nonetheless, other substantial differences can be seen in regions with significant seafloor topography, including the northern part of the Tonga-Kermadec Ridge and the eastern Melanisian Basin. Hence, it is clear from these animations that bathymetry has a significant effect upon the surface wave train, with particularly strong differences appearing in synthetics near regions of strong variation in seafloor topography, as would be expected. Whilst the effects on the body waves again appear to be relatively small in comparison, this suggests that analysis of observed surface waves should involve consideration of the roughness of the seafloor over the source–receiver propagation path. 6 PROPAGATION THROUGH THE WATER COLUMN In AxiSEM3D, it is also possible to add receivers into the water column, which may be of interest in understanding the data recorded by MERMAIDS or at stations of the International Monitoring System for nuclear non-proliferation. To illustrate how waveforms within the water column vary as a function of depth, Fig. 12 shows the seismograms at regular depth intervals in the water column ‘above’ Wake Island. Of course, the seismometer on Wake Island is not actually underwater, but due to our relocation of the station to the seafloor, has an ersatz water column above it. Figure 12. Open in new tabDownload slide Seismograms in the water column above station IU.WAKE (Wake Island) at regularly spaced intervals. The vertical component, which is continuous across the seafloor interface, is shown. Figure 12. Open in new tabDownload slide Seismograms in the water column above station IU.WAKE (Wake Island) at regularly spaced intervals. The vertical component, which is continuous across the seafloor interface, is shown. The largest signals by amplitude in the vertical component occur at the surface, with notable changes in waveform shaped descending through the water column. A detailed interpretation of these waveforms is beyond the scope of this paper, though the most significant body wave codas appear higher up in the water column. The stations in the middling depths appear to have the smoothest waveforms, which may be due to complex reflections and scattering being more significant at the seafloor and surface. 7 THE EFFECTS OF THE OCEAN ON OBSERVED DATA We now investigate whether the modifications to the waveforms induced by the addition of a fluid ocean and bathymetry are actually notable in observations. The comparison is made to both our new fluid ocean implementation with bathymetry, and also to the ‘current standard’, a model with a solid surface with topography upon it. Given the results presented earlier regarding the inapplicability of the ocean loading formulation at intermediate-to-high frequencies, we do not consider application of an ocean load to avoid conflating physical effects (the effects of the ocean’s additional mass) with numerical errors (arising from the load’s inaccuracy at the resolutions in question). As detailed previously, Moho undulation and tomographic models are also included. As the 5 s dominant period used in the generation of synthetic waveforms falls in the middle of the ocean microseism band (Longuet-Higgins 1950), we filter both data and synthetics using a log-Gabor filter with a central frequency of 25 s and focus on the surface waves. Although this removes some of the higher frequency content of both synthetics and data, the prevalence of ocean noise at these island stations at shorter periods means that meaningful comparison would be extremely challenging. The choice of a log-Gabor filter over a bandpass filter is made as its response is more suited to removing the numerical artefacts at above the maximum mesh-resolved frequency whilst also suppressing the intermediate-frequency ocean microseismic noise. It should be emphasized that the 25 s central frequency of the filter’s instantaneous impulse response in this section is not directly comparable to the 25 s panel that in Fig. 6—there it is the high-corner frequency of the bandpass filter above which little-to-no energy is transmitted, whilst here it is the frequency of maximal filter response for a relatively wide passband filter. Fig. 13 shows the data overlaid with synthetics in the no ocean case on the left, and in the fluid ocean with bathymetry case on the right. There is a notable change in the qualitative fit between data and synthetics. Given that at 25 s resolution we expect the most significant changes to appear in the surface waves, we focus on these in this section. Figure 13. Open in new tabDownload slide Comparisons between observed data (dashed lines) and synthetics (solid lines) for four island stations, with both data and synthetics log-Gabor filtered at 25 s and the synthetics convolved with the derived source-time function. Each trace is normalized to its peak P-wave amplitude, with the vertical component shown. Figure 13. Open in new tabDownload slide Comparisons between observed data (dashed lines) and synthetics (solid lines) for four island stations, with both data and synthetics log-Gabor filtered at 25 s and the synthetics convolved with the derived source-time function. Each trace is normalized to its peak P-wave amplitude, with the vertical component shown. It can be seen that the longer duration of the surface wave train, together with what appears to be a scattering-induced coda, are notable in the fluid ocean case. The length of the initial Rayleigh wave train is also better matched. The most substantial discrepancy between data and synthetics is at station IU.RAO in New Zealand’s Kermadec Islands. This is likely attributable to the station’s location on the edge of the Tonga Ridge, a region of significant tectonic and bathymetric complexity which is not accurately or reliably reproduced at the coarse sampling level of the existing crustal models used here. Thus, from a qualitative inspection, we conclude that the contributions of a fluid ocean with realistic bathymetry account for more of the observed waveform features visible in data than are predicted from synthetics which include only topography on the solid surface without a fluid ocean, at least for this particular set of source–receiver pairs. In other words we suggest that this implementation is a more reliable representation of the actual wave propagation physics. 8 DISCUSSION 8.1 Computational cost The efficiency of AxiSEM3D’s method is such that high-frequency simulations as described in this paper may be carried out using an intermediate-size cluster. The 30-min duration simulations at 10 s (as shown in the animations), with both crustal and tomographic models, require approximately 6000 core-hours on a Cray XC30 architecture. The addition of a 3 km fluid ocean with a flat seafloor increases the cost by approximately 90 per cent, in part to the modest increase in the number of elements needed (from ∼124 300 to ∼142 500) to accommodate the ocean’s volume, but mostly due to the 30 per cent decrease in the minimum time step to accommodate its lower sound speed. The latter effect is more significant because AxiSEM3D uses the same time step across the mesh, and if the limiting global time step is set in the ocean, it acts as a constraint on the time step everywhere within the simulation volume. At 5 s resolution (as in the seismogram plots), the corresponding cost for a simulation without an ocean layer is approximately 46 000 core-hours, with a relatively smaller increase (65 per cent) when a fluid layer is added. This illustrates the addition of the ocean layer is less of an additional computational burden at higher frequencies. The increase in cost associated with the addition of bathymetry is strongly model-dependent, and evaluating it is somewhat more involved. At 10 s resolution in the models discussed above, replacing the flat seafloor with an undulating one increases the cost by approximately 10 per cent. In the 5 s case there is actually a cost decrease of approximately 6 per cent associated with the addition of bathymetry. This may be somewhat counterintuitive, given that the particle relabelling transformation is computationally intensive and in-and-of itself always increase the computational cost, but can be understood by considering the deformation of elements at the location where the limiting time step is set. In the 5 s case without bathymetry (but with undulation along the crust–mantle boundary), the minimum time step is set at the Moho at a position approximately 111° epicentral distance from the source. Because of the axisymmetric way in which we formulate the problem, a more exact geographic location cannot be identified. The addition of bathymetry deforms the elements between the Moho and seafloor again, in order to ensure that both interfaces are properly honoured. In the examples considered here, the element in which the limiting time step was originally being set is stretched to a degree that it is no longer has the smallest time step globally. Instead, the new limiting time step is set at a location 106° from the source, though still at the Moho. This new element has a time step of 0.0236 s, ∼8 per cent higher than the original value of 0.0218 s. Just as the decrease in minimum time step (of approximately 30 per cent) with the addition of the ocean increased the overall simulation cost by a comparable value, this increase in dt decreases the simulation cost by a corresponding factor. The reason that the cost decrease is not exactly the same (6 per cent rather than 8 per cent) is due to the fact that the additional particle relabelling involves extra computation. In the simulations presented here, this is clearly advantageous; but it should be noted that this effect may be somewhat esoteric and associated with our choice of 3-D models and the exact resolution chosen (indeed, it occurs at 5 s but not 10 s). Hence it is not generally true that adding bathymetry to simulations will decrease the computational cost. However, at high frequencies with high resolution implementations of the Moho (which necessitate very small elements), the addition of a second undulating boundary may occasionally yield such cost savings, though it should be noted that the demands on computer memory will be higher regardless. 8.2 Cost scaling For a fixed structural model in AxiSEM3D it is possible to derive approximate cost scalings with frequency, which can be compared to other tools such as SPECFEM. In the scenarios considered in this paper, the seismic wavelength is substantially smaller than the characteristic structural scale associated with the most densely sampled feature (the bathymetry). As such, an approximate scaling relationship with frequency ω and azimuthal order Nu[derived by Leng et al. (2016, 2019)] is given by O(ω3 · Nu). This linear behaviour in Nu holds regardless of whether the increase is global (across the entirety of the planet’s radius) or confined to only the upper crust and mantle, though a localized increase in Nu is of course cheaper than a global one. Thus, AxiSEM3D offers significant flexibility and can accommodate suitable 3-D heterogeneities of arbitrary complexity with a favourable cost scaling, so long as the seismic wavelength is shorter than the structural scales. If the model becomes extremely complicated such that the seismic wavelength is comparable to the characteristic scale of 3-D features, the scaling tends toward O(ω4 · Nu), more akin to that of full 3-D methods such as SPECFEM. For reference, Leng et al. (2016) conduct a full cost comparison against SPECFEM3D-GLOBE, and at 10 s with a full 3-D tomographic model (similar to the scenario here but without the fluid ocean and at lower resolution), the speed-up is ∼100. The advantage increases with frequency but decreases with model complexity. We also provide an estimate of the computational resources required for a global t-wave simulation, in reference to the scenarios discussed in Section 1 above—though it should be noted again that this estimate is strongly model dependent. Assuming that simulations were carried out at around 1 Hz (likely the lowest frequency suitable for t-phase modelling), the flexibility of AxiSEM3D is such that the Nu profile could be substantially reduced below the Moho, a regional scale mesh used, and mantle tomographic models removed. Depending on the scales of propagation required, we estimate that this would have a computational cost of 105−106 processor-hours. Increasing the frequency to 10 Hz would increase this by a factor of approximately ∼100. 8.3 Potential cost savings All simulations conducted thus far in this paper have made use of a ‘full’ global mesh and an azimuthal expansion of Nu which is more than high enough to capture both body and surface wave interactions with 3-D structures. Using the flexibility of AxiSEM3D, we now consider the effect of the cost-saving measures discussed in Section 8.2. Such modifications may include using a regional mesh (one which does not solve the wave equation on a complete sphere using a ‘D’ shaped mesh, but rather uses a wedge-shaped ‘chunk’ to simulate a smaller area), or localizing the regions of high Nu such that phases of interest are accurately simulated whilst others are not. As an example, Fig. 14 compares the generated synthetics for two different sets of simulation parameters at the station II.KWAJ, with the profiles given in Table 4. As before, all depths are relative to the ocean surface at radius 6371 km. Figure 14. Open in new tabDownload slide Synthetics for the station II.KWAJ in the Marshall Islands, showing the sensitivity to a reduction in the value of Nu. Note the largest discrepancies in the surface-multiple body waves, and in the surface wave coda. Figure 14. Open in new tabDownload slide Synthetics for the station II.KWAJ in the Marshall Islands, showing the sensitivity to a reduction in the value of Nu. Note the largest discrepancies in the surface-multiple body waves, and in the surface wave coda. Table 4. In each scenario, four values are specified in the AxiSEM3D Nu profile. The ‘Greater Nu value’ refers to the constant value of Nu used between the surface and the depth given by ‘Greater Nu depth bound’, whilst the ‘Lesser Nu value’ refers to the constant value of Nu used between the ‘Lesser Nu depth bound’ and the planet’s core. Linear interpolation between the Greater and Lesser Nu values is used between the two boundaries. As an example, the ‘Low’ profile uses Nu = 1000 in the top 40 km of radius, and Nu = 40 below 30 km depth, and a linear scaling in between. Regime . ‘High cost’ . ‘Low cost’ . Mesh colatitude 180° 90° Greater Nu value 1500 1000 Lesser Nu value 200 40 Greater Nu depth bound (km) 100 30 Lesser Nu depth bound (km) 400 400 Regime . ‘High cost’ . ‘Low cost’ . Mesh colatitude 180° 90° Greater Nu value 1500 1000 Lesser Nu value 200 40 Greater Nu depth bound (km) 100 30 Lesser Nu depth bound (km) 400 400 Open in new tab Table 4. In each scenario, four values are specified in the AxiSEM3D Nu profile. The ‘Greater Nu value’ refers to the constant value of Nu used between the surface and the depth given by ‘Greater Nu depth bound’, whilst the ‘Lesser Nu value’ refers to the constant value of Nu used between the ‘Lesser Nu depth bound’ and the planet’s core. Linear interpolation between the Greater and Lesser Nu values is used between the two boundaries. As an example, the ‘Low’ profile uses Nu = 1000 in the top 40 km of radius, and Nu = 40 below 30 km depth, and a linear scaling in between. Regime . ‘High cost’ . ‘Low cost’ . Mesh colatitude 180° 90° Greater Nu value 1500 1000 Lesser Nu value 200 40 Greater Nu depth bound (km) 100 30 Lesser Nu depth bound (km) 400 400 Regime . ‘High cost’ . ‘Low cost’ . Mesh colatitude 180° 90° Greater Nu value 1500 1000 Lesser Nu value 200 40 Greater Nu depth bound (km) 100 30 Lesser Nu depth bound (km) 400 400 Open in new tab The low case shows very little discrepancy compared to the high one, indicating that the solution is well converged and the 3-D structures present along this particular source–receiver path are adequately sampled. The only apparent differences are in the deeply penetrating body waves, as may be expected from their non-negligible sensitivity at depths where the sampling of Nu has been significantly reduced; and in the surface wave coda. The differences in the coda are likely caused by the interacting reflections from the edge of the regional-scale mesh. It should be noted that the degree of discrepancy between the different traces is likely to be strongly azimuth-dependant—a source–receiver path with smooth bathymetric and Moho undulation will require comparably lower Nu to achieve the same degree of convergence than one with sharp variations in the interfaces. This must be accounted for when optimizing the Nu profile for a particular study; and of course regional meshes with a reflecting boundary (as here) are not suitable for studying surface wave codas due to contamination of the signal by reflections. In this scenario substantial reductions of Nu are not, therefore, suitable for high-resolution studies of body wave surface multiples (such as PP or PPP), or generation of long-duration synthetics (where spurious reflections from the edge of the truncated mesh become problematic), but with careful optimization of the Nu profile significant computational savings can be made if only specific phases are being investigated. For example, the ‘Low’ profile reduces the computational cost by 97 per cent at 5 s. This enormous speedup is associated both with the reduction in Nu and the decrease in the size of the computational domain, the latter of which also increases the time step as the previous time step limiting element in the mesh (at 106° distance) is not included. 8.4 Limitations As mentioned previously, the use of a global time step means that the addition of an ocean significantly increases the computational cost of simulations; the positive caveat being that for a given model the degree of temporal oversampling decreases with increasing frequency. As a result, higher-frequency simulations are more expensive than lower-frequency ones, and simulations with an ocean layer are more expensive than ones without; but the efficiency savings from moving to higher frequencies in AxiSEM3D are greater in the case with an ocean than the case without. The addition of localized time stepping (e.g. Rietmann et al. (2015)) is one route through which the computational cost can be further decreased without sacrificing simulation accuracy. The fact that the ocean must cover the entire planet is a clear current limitation of this method. The comparison to observations in Section 7 indicates that this is still an advance upon current numerical modelling methods, which can be justified physically on the grounds that the increase in the complexity of wave propagation physics captured with a fluid ocean, can, with careful use, outweigh the reductions in reliability associated with the relocation of stations on isolated islands to the seafloor. A more realistic representation of ocean–continent boundaries would be far more flexible in its application, and hence remains a long-term development goal. 8.5 Extensions in a global ocean model In simulations constrained to use a global ocean model, the simplest increase in the realism (and hence of reliability) of the synthetics would be a reduction in the minimum ocean depth (here 500 m) over the continents. This would mean that at a given seismic period, the sensitivity to the ersatz fluid layer would be reduced. Quantifying the exact degree to which this would change the synthetics at a given location is difficult, given that the effects of bathymetry are significant and highly spatially inhomogeneous, and that the addition of the fluid layer affects each component and seismic phase differently. Nonetheless, some insight may be gleaned from Fig. 12. A rough estimate based on the convergence of the fluid surface layer case to the no-fluid surface layer case with increasing period (see Supporting Information) suggests that where the ocean depth is at least an order of magnitude less than the seismic wavelength, the effects of the ersatz fluid layer can be completely neglected. At 5 s resolution, this would equate to the ocean being no deeper than ∼200 m over continents, which would decrease the time step by factor ∼2.5 and hence increase the computational costs reported in Section 8.1 by roughly the same factor. This option may be especially attractive at higher frequencies, where as discussed in Section 8.4 the inclusion of the ocean layer is less of an ‘inefficiency’. Nonetheless, if we restrict ourselves to purely oceanic source–receiver paths, various interesting phenomena not considered here can still be simulated on global scales. These could include scenarios either without or with a source in the fluid ocean. In the first case, an exploration of the effects of ocean bathymetry on studies of source properties using intermediate frequency (10–20 s) body waves would pertinent, given that such effects are of particular concern when attempting to understand the dynamics of the largest earthquakes, which occur at subduction zones with rough bathymetry. Alternatively, a source in the fluid ocean could be used to study the generation of the primary ocean microseism by considering the propagation of acoustic waves from a near-surface pressure source down towards an area of rough bathymetry. With the addition of a stratified or globally varying water column structure and the use of a regional-scale mesh with an optimized Nu profile, simulation of t-phase propagation in the SOFAR channel is in theory achievable, if the computational resources needed to reach the necessary high frequencies (∼5 Hz) are available. Sources here could either be in the water column, for example representing illicit nuclear tests; or in the solid earth to study conversion of earthquake-generated P waves. Care would need to be taken here to ensure that bathymetry along the propagation paths here is realistic, as t-phases are more sensitive to small-scale bathymetric variations than body waves, and hence our relocation of small islands to be underwater may misrepresent the degree of ‘propagation blockage’ which such waves experience. AxiSEM3D is ideally suited to such studies, as the value of Nu in the crust can be adjusted to ensure sufficient sampling of the bathymetry without adding unnecessary computational cost in the deeper mantle. 8.6 Development of a localized ocean formulation The development of a localized, more realistic ocean formulation is also an avenue for potential future exploration. In the AxiSEM3D formulation it would be possible to create a ‘localized’ ocean which is axisymmetric, but which nonetheless honours a realistically sloping seafloor which meets the continental surface at its edges. This can be envisaged as a doughnut-shaped ocean ‘trench’ circling the entire planet about its equator, for example. This would require manual editing of the mesh, but would enable more realistic examination of the effects of the bathymetric profiles of continental shelves, if appropriate source–receiver paths are chosen. Such a setup might, for example, be of interest in understanding the generation of the transverse component of ocean microseismic noise. Going one step further, development of a completely realistic ‘patched’ ocean would be extremely challenging due to the use of global basis functions in the AxiSEM3D Fourier expansion. These basis functions impose the constraint that the same physical principles, such as the boundary conditions, must be the same on any line of azimuth. Hence, the arbitrary switching between fluid and solid areas which would be required for a patched ocean is not at present achievable. It ought be possible to replace the global pseudospectral representation with a localized one, for example one in which the basis functions are discontinuous wavelets or Slepian functions (e.g. Simons 2010). Such a basis would not use globally constrained interpolating functions, and hence could model different boundary conditions at different azimuths. 9 CONCLUSION We have presented a novel, open-source method for simulating global seismic wavefields in a model with a fluid ocean layer and 3-D bathymetry. This method is embedded in the spectral-element method AxiSEM3D, introduced by Leng et al. (2016), which uses a pseudospectral representation in the azimuthal direction to significantly reduce computational cost. Our implementation is benchmarked for a radially symmetric model against the code YSpec to ensure its validity. An evaluation of the ocean-loading formulation previously used in global seismology is conducted, and we find that the approximation is not valid for the entire wave train unless the period reaches longer than 20 s for a 3 km ocean depth. Such findings concur with previous evaluations of this formulation which used a simplified testing setup (Zhou et al. 2016). In a non-radially symmetric (i.e. fully 3-D) model, we are able to investigate the effects of bathymetry through use of a global ocean formulation, and show that seafloor topography introduces substantial modifications to the surface waves. The modification to the body waves, whilst appreciable at 5 s dominant period, is much more slight. A suggestive comparison to data was also made, wherein we find that more of the gross features of the observed surface wave train are reproduced when a fluid, bathymetric ocean is added to the structural model. Whilst this study considers only specific oceanic propagation paths, the prevalence of oceanic source–receiver paths in modern seismic analysis and the importance of interpreting data from remote stations in otherwise poorly sampled areas suggests this implementation can be a useful advance in computational seismic modelling. SUPPORTING INFORMATION depths.png,not_overplot.avi,overplot.avi Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. supplement ACKNOWLEDGEMENTS We are grateful to David Al-Attar, Martin van Driel, Maria Tsekhmistrenko, Kasra Hosseini, Claudia Haindl and Thomas Garth for their technical help on this topic. Calculations were performed on the UK National Supercomputer (ARCHER), under TNM’s NERC/EPSRC grant. KL received support from NERC grant NE/R012199/1. BF is funded by the NERC Doctoral Training Partnership in Environmental Research at the University of Oxford. Additional support from the STFC-AURORA grant ST/S001379/1 is also gratefully received. We also thank Dr Wenbo Wu, Dr Christian Boehm and Dr Sidao Ni for their helpful comments in reviewing this manuscript. REFERENCES Abràmoff M.D. , Magalhães P.J., Ram S.J., 2004 . Image processing with imageJ , Biophoton. Int. , 11 ( 7 ), 36 – 42 . OpenURL Placeholder Text WorldCat Afanasiev M. , Boehm C., Van Driel M., Krischer L., Rietmann M., May D.A., Knepley M.G., Fichtner A., 2019 . Modular and flexible spectral-element waveform modelling in two and three dimensions , J. geophys. Int. , 216 ( 3 ), 1675 – 1692 . Google Scholar Crossref Search ADS WorldCat Al-Attar D. , Crawford O., 2016 . Particle relabelling transformations in elastodynamics , J. geophys. Int. , 205 ( 1 ), 575 – 593 . Google Scholar Crossref Search ADS WorldCat Al-Attar D. , Woodhouse J.H., 2008 . Calculation of seismic displacement fields in self-gravitating earth models-applications of minors vectors and symplectic structure , J. geophys. Int. , 175 ( 3 ), 1176 – 1208 . Google Scholar Crossref Search ADS WorldCat Amante C. , Eakins B., 2009 . ETOPO1 1 arc-minute global relief model: procedures, data sources and analysis , NOAA Technical Memorandum, NESDIS NGDC-24 . OpenURL Placeholder Text WorldCat Ardhuin F. , Gualtieri L., Stutzmann E., 2015 . How ocean waves rock the Earth: two mechanisms explain microseisms with periods 3 to 300s , Geophys. Res. Lett. , 42 ( 3 ), 765 – 772 . Google Scholar Crossref Search ADS WorldCat Basini P. , Nissen-Meyer T., Boschi L., Casarotti E., Verbeke J., Schenk O., Giardini D., 2013 . The influence of nonuniform ambient noise on crustal tomography in Europe , Geochem. Geophys. Geosyst. , 14 ( 5 ), 1471 – 1492 . Google Scholar Crossref Search ADS WorldCat Beyreuther M. , Barsch R., Krischer L., Megies T., Behr Y., Wassermann J., 2010 . ObsPy: a python toolbox for seismology , Seismol. Res. Lett. , 81 ( 3 ), 530 – 533 . Google Scholar Crossref Search ADS WorldCat Blackman D.K. , Orcutt J.A., Forsyth D.W., 1995 . Recording teleseismic earthquakes using ocean-bottom seismographs at mid-ocean ridges , Bull. seism. Soc. Am. , 85 ( 6 ), 1648 – 1664 . OpenURL Placeholder Text WorldCat Bohnenstiehl D.W.R. , Dziak R.P., Matsumoto H., Lau T.K., 2013 . Underwater acoustic records from the March 2009 eruption of Hunga Ha’apai-Hunga Tonga volcano in the Kingdom of Tonga , J. Volc. Geotherm. Res. , 216 ( 3 ), 12 – 24 . Google Scholar Crossref Search ADS WorldCat Bottero A. , Cristini P., Komatitsch D., Asch M., 2016 . An axisymmetric time-domain spectral-element method for full-wave simulations: application to ocean acoustics , J. acoust. Soc. Am. , 140 ( 5 ), 3520 – 3530 . Google Scholar Crossref Search ADS PubMed WorldCat Chapp E. , Bohnenstiehl D.R., Tolstoy M., 2005 . Sound-channel observations of ice-generated tremor in the Indian Ocean , Geochem. Geophys. Geosyst. , 6 ( 6 ), doi:10.1029/2004GC000889. OpenURL Placeholder Text WorldCat Cristini P. , Komatitsch D., 2012 . Some illustrative examples of the use of a spectral-element method in ocean acoustics , J. acoust. Soc. Am. , 131 ( 3 ), 229 – 235 . Google Scholar Crossref Search ADS WorldCat Crotwell H.P. , Owens T.J., Ritsema J., 1999 . The TauP Toolkit: flexible seismic travel-time and ray-path utilities , Seismol. Res. Lett. , 70 ( 2 ), 154 – 160 . Google Scholar Crossref Search ADS WorldCat Davy C. , Barruol G., Fontaine F.R., Sigloch K., Stutzmann E., 2014 . Tracking major storms from microseismic and hydroacoustic observations on the seafloor , Geophys. Res. Lett. , 41 ( 24 ), 8825 – 8831 . Google Scholar Crossref Search ADS WorldCat de Groot-Hedlin C.D. , Orcutt J.A., 2001 . Excitation of T -phases by seafloor scattering , J. acoust. Soc. Am. , 109 ( 5 ), 1944 – 1954 . Google Scholar Crossref Search ADS PubMed WorldCat Dréo R. , Bouffaut L., Leroy E., Barruol G., Samaran F., 2019 . Baleen whale distribution and seasonal occurrence revealed by an ocean bottom seismometer network in the Western Indian Ocean , Deep-Sea Res. Part II: Top. Stud. Oceanogr. , 161 , 132 – 144 . Google Scholar Crossref Search ADS WorldCat Dushaw B.D. et al. ., 1999 . Multimegameter-range acoustic data obtained by bottom-mounted hydrophone arrays for measurement of ocean temperature , IEEE J. Ocean. Eng. , 24 ( 2 ), 202 – 214 . Google Scholar Crossref Search ADS WorldCat Dziak R.P. et al. ., 2004 . P- and T-wave detection thresholds, Pn velocity estimate, and detection of lower mantle and core P-waves on ocean sound-channel hydrophones at the Mid-Atlantic Ridge , Bull. seism. Soc. Am. , 94 ( 2 ), 665 – 677 . Google Scholar Crossref Search ADS WorldCat French S.W. , Romanowicz B.A., 2014 . Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography , J. geophys. Int. , 199 ( 3 ), 1303 – 1327 . Google Scholar Crossref Search ADS WorldCat Grob M. , Maggi A., Stutzmann E., 2011 . Observations of the seasonality of the Antarctic microseismic signal, and its association to sea ice variability , Geophys. Res. Lett. , 38 ( 11 ), doi:10.1029/2011GL047525. OpenURL Placeholder Text WorldCat Igel H. , 2016 . Computational Seismology , Oxford Univ. Press . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Komatitsch D. , Tromp J., 2002 . Spectral-element simulations of global seismic wave propagation – II. Three-dimensional models, oceans, rotation and self-gravitation , J. geophys. Int. , 150 ( 1 ), 303 – 318 . Google Scholar Crossref Search ADS WorldCat Kristeková M. , Kristek J., Moczo P., 2009 . Time-frequency misfit and goodness-of-fit criteria for quantitative comparison of time signals , J. geophys. Int. , 178 ( 2 ), 813 – 825 . Google Scholar Crossref Search ADS WorldCat Laske G. , Masters G., Ma Z., Pasyanos M., 2013 . Update on CRUST1.0—a 1-degree global model of Earth’s crust , in Geophys. Res. Abstracts , EGU . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lay T. , Rhode A., 2019 . Evaluating the updip extent of large megathrust ruptures using pcoda levels , Geophys. Res. Lett. , 46 ( 10 ), 5198 – 5206 . Google Scholar Crossref Search ADS WorldCat Lay T. , Liu C., Kanamori H., 2019 . Enhancing tsunami warning using P wave coda , J. geophys. Res. , 124 ( 10 ), 10 583 – 10 609 . Google Scholar Crossref Search ADS WorldCat Leng K. , Nissen-Meyer T., van Driel M., 2016 . Efficient global wave propagation adapted to 3-D structural complexity: a pseudospectral/spectral-element approach , J. geophys. Int. , 207 ( 3 ), 1700 – 1721 . Google Scholar Crossref Search ADS WorldCat Leng K. , Nissen-Meyer T., van Driel M., Hosseini K., Al-Attar D., 2019 . AxiSEM3D: broadband seismic wavefields in 3-D global Earth models with undulating discontinuities , J. geophys. Int. , 217 ( 3 ), 2125 – 2146 . Google Scholar Crossref Search ADS WorldCat Longuet-Higgins M.S. , 1950 . A theory of the origin of microseisms , Phil. Trans. R. Soc., A , 243 ( 857 ), 1 – 35 . OpenURL Placeholder Text WorldCat Mazoyer C. , Guillon L., Guennou C., Jamet G., Royer J.-Y., 2013 . T-wave generation and propagation: a comparison between data and spectral element modeling , J. acoust. Soc. Am., 134 ( 4 ), 3376 – 3385 . OpenURL Placeholder Text WorldCat McDonald M.A. , Hildebrand J.A., Webb S.C., 1995 . Blue and fin whales observed on a seafloor array , J. acoust. Soc. Am. , 98 ( 2 ), 712 – 721 . Google Scholar Crossref Search ADS PubMed WorldCat Mitchell B.J. , 2002 . Monitoring the comprehensive nuclear-test-ban treaty , Pure appl. Geophys. , 159 ( 4 ), 619 – 620 . Google Scholar Crossref Search ADS WorldCat Nishida K. , Takagi R., 2016 . Teleseismic S wave microseisms , Science , 353 ( 6302 ), 919 – 921 . Google Scholar Crossref Search ADS PubMed WorldCat Nissen-Meyer T. , Fournier A., Dahlen F.A., 2007 . A two-dimensional spectral-element method for computing spherical-earth seismograms – I. Moment-tensor source , J. geophys. Int. , 168 ( 3 ), 1067 – 1092 . Google Scholar Crossref Search ADS WorldCat Nissen-Meyer T. , Van Driel M., Stähler S.C., Hosseini K., Hempel S., Auer L., Colombi A., Fournier A., 2014 . AxiSEM: Broadband 3-D seismic wavefields in axisymmetric media , Solid Earth , 5 , 425 – 445 . Google Scholar Crossref Search ADS WorldCat Okal E.A. , 2008 . The generation of T waves by earthquakes , Adv. Geophys. , 49 ( 07 ), 1 – 65 . OpenURL Placeholder Text WorldCat Peter D. et al. ., 2011 . Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes , J. geophys. Int. , 186 ( 2 ), 721 – 739 . Google Scholar Crossref Search ADS WorldCat Qian Y. , Wei S., Wu W., Zeng H., Coudurier-Curveur A., Ni S., 2019 . Teleseismic waveform complexities caused by near trench structures and their impacts on earthquake source study: application to the 2015 illapel aftershocks (Central Chile) , J. geophys. Res. , 124 ( 1 ), 870 – 889 . Google Scholar Crossref Search ADS WorldCat Rietmann M. , Peter D., Schenk O., Ucar B., Grote M., 2015 . Load-balanced local time stepping for large-scale wave propagation , in Proceedings - 2015 IEEE 29th International Parallel and Distributed Processing Symposium, IPDPS 2015 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Robert Engdah E. , Van Hilst R.D., Buland R., 1998 . Global teleseismic earthquake relocation with improved travel times and procedures for depth determination , Bull. seism. Soc. Am. , 88 ( 3 ), 722 – 743 . OpenURL Placeholder Text WorldCat Sasorova E.V. , Levin B.W., Morozov V.E., Didenkulov I.N., 2005 . Hydro-acoustic monitoring on the Kamchatka Shelf: a possibility of early location of oceanic earthquake and local tsunami warning , in Tsunamis. Advances in Natural and Technological Hazards Research , Vol. 23 , pp. 305 – 317 ., ed. Satake K., Springer . Google Scholar Crossref Search ADS WorldCat Shapiro N.M. , Campillo M., Stehly L., Ritzwoller M.H., 2005 . High-resolution surface-wave tomography from ambient seismic noise , Science , 307 ( 5715 ), 1615 – 1618 . Google Scholar Crossref Search ADS PubMed WorldCat Simons F.J. , 2010 . Slepian functions and their use in signal estimation and spectral analysis , in Handbook of Geomathematics , Springer . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Stutzmann E. , Schimmel M., Patau G., Maggi A., 2009 . Global climate imprint on seismic noise , Geochem. Geophys. Geosyst. , 10 ( 11 ), doi:10.1029/2009GC002619. OpenURL Placeholder Text WorldCat Tolstoy I. , Ewing M., 1950 . The T phase of shallow-focus earthquakes , Bull. seism. Soc. Am. , 40 ( 1 ), 25 – 51 . OpenURL Placeholder Text WorldCat Vallée M. , Douet V., 2016 . A new database of source time functions (STFs) extracted from the SCARDEC method , Phys. Earth Planet. Inter. , 257 , 149 – 147 . Google Scholar Crossref Search ADS WorldCat van Driel M. , Nissen-Meyer T., 2014a . Seismic wave propagation in fully anisotropic axisymmetric media , J. geophys. Int. , 199 ( 2 ), 880 – 893 . Google Scholar Crossref Search ADS WorldCat van Driel M. , Nissen-Meyer T., 2014b . Optimized viscoelastic wave propagation for weakly dissipative media , J. geophys. Int. , 1999 ( 2 ), 1078 – 1093 . Google Scholar Crossref Search ADS WorldCat Wu W. , Ni S., Zhan Z., Wei S., 2018 . An SEM-DSM three-dimensional hybrid method for modelling teleseismic waves with complicated source-side structures , J. geophys. Int. , 215 , 133 – 154 . Google Scholar Crossref Search ADS WorldCat Wu Z. , Lay T., Ye L., 2020 . Shallow megathrust slip during large earthquakes that have high P coda levels , J. geophys. Res. , 125 ( 1 ), doi:10.1029/2019JB018709. OpenURL Placeholder Text WorldCat Yue H. , Castellanos J.C., Yu C., Meng L., Zhan Z., 2017 . Localized water reverberation phases and its impact on backprojection images , Geophys. Res. Lett. , 44 ( 19 ), 9573 – 9580 . Google Scholar Crossref Search ADS WorldCat Zhou Y. , Ni S., Chu R., Yao H., 2016 . Accuracy of the water column approximation in numerically simulating propagation of teleseismic PP waves and Rayleigh waves , J. geophys. Int. , 206 ( 2 ), 1315 – 1326 . Google Scholar Crossref Search ADS WorldCat © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Oceanic high-frequency global seismic wave propagation with realistic bathymetry
JF - Geophysical Journal International
DO - 10.1093/gji/ggaa248
DA - 2020-08-01
UR - https://www.deepdyve.com/lp/oxford-university-press/oceanic-high-frequency-global-seismic-wave-propagation-with-realistic-IPioCMMYPR
SP - 1178
VL - 222
IS - 2
DP - DeepDyve
ER -