Oceanic high-frequency global seismic wave propagation 1 with realistic bathymetry 2

6 We present a new approach to simulate high-frequency seismic wave propagation in 7 and under the oceans. Based upon AxiSEM3D (Leng et al. 2019), this method sup- 8 ports a ﬂuid ocean layer, with associated water-depth phases and seaﬂoor topography 9 (bathymetry). The computational efﬁciency and ﬂexibility of this formulation means that 10 high-frequency calculations may be carried out with relatively light computational loads. 11 A validation of the ﬂuid ocean implementation is shown, as is an evaluation of the oft-used 12 ocean loading formulation, which we ﬁnd breaks down at longer periods than was previ- 13 ously believed. An initial consideration of the effects of seaﬂoor bathymetry on seismic 14 wave propagation is also given, wherein we ﬁnd that the surface waveforms are signiﬁ- 15 cantly modiﬁed in both amplitude and duration. When compared to observed data from 16 isolated island stations in the Paciﬁc, synthetics which include a global ocean and seaﬂoor 17 topography appear to more closely match the observed waveform features than synthetics 18 generated from a model with topography on the solid surface alone. We envisage that 19 such a method will be of use in understanding the new and exciting ocean-bottom and 20 ﬂoating seismometer datasets now being regularly collected. 21 N these


INTRODUCTION
given the majority of our planet's surface covered by water, and the resulting implication that an 126 arbitrary source-receiver path at teleseismic distances is more likely than not to include an oceanic 127 section. Such an approach is also justified on observational grounds -by undertaking a comparison to 128 recorded data from remote island stations (which provide valuable global coverage in sparsely sampled 129 areas and by their nature are strongly affected by the presence of an ocean), we give examples where 130 a better match to synthetics occurs when a fluid ocean is included than when it is not. 131 Such a methodology should enable higher-frequency simulations to be carried out in an oceanic 132 seismology context, with arbitrary structural complexity in the underlying solid Earth. In this paper we 133 restrict further detailed discussion to examination of bathymetric effects and water-depth phases; how-134 ever extensions of our methodology to other fluid seismology contexts (the atmosphere or localised 135 but axisymmetric oceans) are also possible. We first consider a spherical Earth model that consists of a solid domain Ω S (with density ρ and 143 elasticity tensor C) and a fluid domain Ω F (also with density ρ and bulk modulus κ). These two 144 domains can be separated by several solid-fluid interfaces, collectively denoted Σ, such as the ocean 145 floor, the core-mantle boundary and the inner core boundary. 146 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 wheren denotes the outward-pointing normal of a solid-fluid interface, f is a body force source, u is the displacement vector, and w is an arbitrary, vector-valued test function. In the fluid domain, we define a scalar potential χ such that u = ρ −1 ∇χ in Ω F , which is a descriptor for the dynamic pressure through Ω F as ∂ 2 t χ = −P . (2) Here, w is an arbitrary scalar-valued test function. The rank-two identity tensor I, which appears to be 147 somewhat redundant here, is included to enable more natural extension to an aspherical earth model 148 with undulating boundaries in eq. (7). Note that Ω S and Ω F are coupled by the two surface integrals 149 over Σ.

159
Such a pressure-free boundary condition is perfectly reflecting in an analytical sense, which is 160 justified physically on the grounds that the acoustic impedance of air is very much greater than that 161 of water so coupling from the ocean surface into the atmosphere can reasonably be neglected (though, 162 if desired, it in can be accounted for in our formulation so long as the atmosphere is included in the 163 mesh). At the seafloor boundary, the continuity of traction and normal velocity are implemented in the 164 same way as at the outer core boundary. 165 When 3D bathymetry is incorporated, eqs.
(1) and (2) remain valid, but Ω S , Ω F and Σ become  Given a reference spherical configurationΩ, and a deformed configuration Ω which is homeomorphic toΩ, the radial coordinate of a particle at r inΩ is shifted alongr by an amount τ (r) inΩ such that ξ(r) = r + τ (r)r.
Note that this transformation only shifts the boundary inward or outward alongr, which is commen-171 surate with the fact that the boundaries in question (e.g. seafloor bathymetry) have different radii at 172 different positions, but do not need to be shifted laterally.
173 ξ, where ξ :Ω → Ω (and hence also ξ :Σ → Σ), represents the undulation mapping, which is a kinematically permissible deformation with associated deformation gradient F. In spherical coordinates, where F = F(r, θ, φ), the deformation gradient can be expressed as: whilst the associated Jacobian for this transformation may be expressed as: The weak form of the undulating, three-dimensional model can thus be established in the reference configurationΩ in the solid (corresponding to eq. (1)) as whilst in the fluid domain (corresponding to eq. (2)) the formulation becomes where the equivalent material parameters are given bỹ κ(r) = κ(ξ(r))|J(τ, r)|, and the equivalent source byf andĨ andñ respectively byĨ The solutions to eqs. (6) and (7), in u and χ respectively, are related to the solutions of eqs.
This formulation is sufficient to describe all variables needed to undertake the task of finding  The choice of the highest term which must be included in the azimuthal Fourier expansion (N u ) 188 depends on the complexity of the model. In a radially symmetric (1D) model, this representation 189 becomes analytic for a second-order (quadropolar) moment tensor when N u = 2. For more complex 190 models with significant 3D structures, N u can be increased as needed to capture azimuthal structures.

195
As per eqs. (16) and (17) AxiSEM3D. An example of an explicitly meshed surface fluid layer is shown in Fig. 1. Note that 202 we use the terms 'bathymetry' and 'seafloor topography' interchangeably to describe the variation in 203 ocean depth with position.

204
AxiSEM3D supports arbitrary modifications of the seismic profile used in meshing (as generated    The source parameters used are given in Table 1. In YSpec the maximum angular degree used is

268
According to the Kristekovaćlassification, both the envelope and phase misfits are 'excelllent. Thus, 269 we conclude that the implementation of the fluid surface layer in AxiSEM3D has been verified.

270
For completeness, we also demonstrate that the effects of the oceans become negligible at long 271 seismic periods, justifying their lack of inclusion in low-frequency seismic simulations. Fig. 4 shows 272 the convergence of the fully fluid ocean synthetics to those from the PREM case without ocean (and 273 without ocean loading) at 50 seconds dominant period. Apart from the differences in the structural 274 model, the simulation parameters are otherwise identical to those given in Table 1.  Such an approach is similar in aim to that used by Zhou et al. (2016), but with two key differences.

284
Firstly, we use a global, spherical Earth with realistic velocity and density structure, rather than a 285 homogeneous half-space. Secondly, we also consider the effects of the ocean on the entire wavetrain, 286 whilst their studies are restricted to the PP phase and Rayleigh waves.   confirming that in a radially symmetric model with an isotropic ocean, the ocean's effects are con-304 fined to the source-receiver plane. Of course, the ocean depth is not a uniform 3 km across the Earth's 305 surface. For this reason, we also consider the envelope and phase misfits between the two implemen-306 tations for a variety of ocean depths (Fig. 7).   these are not expected to be significant at these seismic periods given the small size of these islands (all 358 have an above-water extent that is much less than 1 • in any direction and hence they are not properly 359 sampled anyhow our implemented bathymetric model). Source parameters are given in Table 2 and 360 receiver locations used in this study in Table 3, with seafloor-projected ray paths shown in Fig. 8. We begin by making a high-resolution comparison of the body waveforms in simulations with and 363 without bathymetry. Such synthetics are presented in Fig. 9.

364
Comparison between the top and bottom panels reinforces the finding that the effects of bathymetry 365   surface waves (leading to larger differences), and their greater sensitivity to near-surface conditions 384 mean that they dominate the difference plots.

385
Video sequences were generated using ImageJ (Abràmoff et al. 2004). We include only four perti-386 nent snapshots of the wavefield in this paper, however the full animations may be found on the 'Oxford 387 Seismology' YouTube channel: https://www.youtube.com/c/SeismologyOxford.

388
In interpreting the differences in Fig. 11, it should be noted that whilst the simulations are con-389 sistent and accurate throughout, the wavefield observed in the sector clockwise from north-west to 390 south-east is most realistic than that toward the south-west, where the scaling of the Australian conti-     Fig. 13 shows the seismograms at regular depth intervals in the water column 'above' Wake

423
Island. Of course, the seismometer on Wake Island is not actually underwater, but due to our relocation 424 of the station to the seafloor, has an ersatz water column above it.  factor. The reason that the cost decrease is not exactly the same (6% rather than 8%) is due to the fact 502 that the additional particle relabelling involves extra computation.

503
In the simulations presented here, this is clearly advantageous; but it should be noted that this ef-504 fect may be somewhat esoteric and associated with our choice of 3D models and the exact resolution 505 chosen (indeed, it occurs at 5 s but not 10 s). Hence it is not generally true that adding bathymetry to 506 simulations will decrease the computational cost. However, at high frequencies with high resolution 507 implementations of the Moho (which necessitate very small elements), the addition of a second undu-508 lating boundary may occasionally yield such cost savings, though it should be noted that the demands 509 on computer memory will be higher regardless. that phases of interest are accurately simulated whilst others are not.

544
As an example, Fig. 15 compares the generated synthetics for two different sets of simulation 545 parameters at the station II.KWAJ, with the profiles given in Table 8.3. As before, all depths are 546 relative to the ocean surface at radius 6371 km. Table 4. In each scenario, four values are specified in the AxiSEM3D N u profile. The 'Greater N u value' refers to the constant value of N u used between the surface and the depth given by 'Greater N u depth bound', whilst the 'Lesser N u value' refers to the constant value of N u used between the 'Lesser N u depth bound' and the planet's core. Linear interpolation between the Greater and Lesser N u values is used between the two boundaries. As an example, the 'Low' profile uses N u = 1000 in the top 40 km of radius, and N u = 40 below 30 km depth, and a linear scaling in between. The low case shows very little discrepancy compared to the high one, indicating that the solution is reflections from the edge of the regional-scale mesh.

553
It should be noted that the degree of discrepancy between the different traces is likely to be strongly 554 azimuth-dependant -a source-receiver path with smooth bathymetric and Moho undulation will re-555 quire comparably lower N u to achieve the same degree of convergence than one with sharp variations 556 in the interfaces. This must be accounted for when optimising the N u profile for a particular study; and 557 of course regional meshes with a reflecting boundary (as here) are not suitable for studying surface 558 wave codas due to contamination of the signal by reflections.

559
In this scenario substantial reductions of N u are not, therefore, suitable for high-resolution studies

608
With the addition of a stratified or globally varying water column structure and the use of a The development of a localised, more realistic ocean formulation is also an avenue for potential future 621 exploration. In the AxiSEM3D formulation it would be possible to create a 'localised' ocean which is 622 axisymmetric, but which nonetheless honours a realistically sloping seafloor which meets the conti-