TY - JOUR
AU - Yablokov,, Aleksander
AB - SUMMARY A promising way to perform seismological studies in the Arctic region is deploying seismic stations on ice floes. The pioneering works by the Alfred Wegener Institute Bremerhaven have demonstrated the efficiency of such floating networks to explore local and regional seismicity and to build 3-D seismic models. However, problems remain, related to the identification of different types of seismic waves, particularly S waves. Here, we perform 2-D and 3-D numerical simulations of seismic waves emitted by an earthquake to explore the possibility of recording different phases on the sea surface. We use different types of simple shear source models, namely strike-slip, vertical displacement and normal faults. In the calculated wave field, we obtain three major types of seismic waves recorded on the sea surface—Pw, Sw and SPw (w denotes an acoustic wave in the water layer)—and numerous multiple waves. The clarity of the recorded phases strongly depends on the type of wave, source mechanism, epicentral distance, thickness of the water layer and depth of the source. For example, the Pw phase is clearest for the strike-slip mechanism, less clear for the normal fault and almost invisible for the vertical displacement. The Sw phase is observable in all of these cases; however, it can be confused with the SPw phase that arrives earlier. In addition, at some distances, the Sw wave interferes with the multiple Pw2 wave and therefore is hardly detectable. In summary, the numerical simulations in a model with a water layer have demonstrated several non-obvious features of wave propagation that should be taken into account when analysing experimental data recorded on ice floes. Computational seismology, Earthquake dynamics, Seismic tomography, Theoretical seismology, Wave propagation 1 INTRODUCTION During the last decade, due to various political and economic reasons, the Arctic region has become a priority target for exploration by several countries that perform multidisciplinary scientific studies. However, the harsh natural conditions in the Arctic make any experiments logistically very difficult and expensive. For example, to study the underground structure beneath the Arctic Ocean, offshore seismic profiling faces many difficulties related to ice drift and severe weather conditions. This situation requires using large icebreakers and other resource-consuming facilities and strongly limits the times and locations for executing experiments. Therefore, a large part of geophysical investigations in the Arctic region is based on remote data, such as satellite measurements of gravity and magnetic fields (e.g. Gaina et al.2014) or data from global seismological networks. For example, there are several regional models of the deep mantle structures beneath the Arctic region, which were constructed based on global databases using the body wave (Jakovlev et al.2012) and surface wave tomography methods (Schaefer & Lebedev 2013). Such approaches are efficient to obtain a general view of the geodynamics of the Arctic as a whole (Koulakov et al.2013). However, they give almost no information about crustal structures that are crucial for understanding local- and regional-scale tectonic processes. Obtaining detailed information about the crust requires installing dense temporary networks of seismic stations, some of which should be placed offshore, which is especially problematic in Arctic conditions. Performing traditional surveys using active and passive seismic sources is exceptionally expensive; therefore, alternative methods for geophysical observations specially adapted for Arctic conditions have been actively developed in recent decades. One such technology is deploying seismic networks on floating ice, which has the potential to considerably improve seismological investigations in Arctic areas (Laderach & Schlindwein 2011). According to worldwide seismological catalogues, most seismicity in the Arctic region is confined along a narrow zone of the Gakkel Ridge, a slow-spreading centre in the Arctic Ocean (Schlindwein et al.2007, 2015). Most other areas are presumed to be aseismic. At the same time, the seismicity in the circum-Arctic area is mostly recorded by seismic stations sparsely distributed in the surrounding lands. Because the distances between these stations may reach hundreds and even thousands of kilometres, they can record only events with magnitudes larger than 3.5–4. Such networks are not able to detect the background seismicity of weaker magnitudes, which is the major indicator for delineating active tectonic boundaries. A lack of information about weak seismicity may be a reason for incorrect estimates of seismic hazards in the Arctic area. In turn, this misestimation may cause some accidents during industrial development and natural resource exploration at high latitudes. A detailed investigation of seismicity is also important for many fundamental scientific problems. For example, recent experiments aimed at recording the background seismicity in the area of the Gakkel Ridge found evidence for ongoing magmatic activity (Muller & Jokat 2000, Schlindwein et al.2005, Schlindwein & Riedel 2010). This discovery provided much information about the occurrence of slow-spreading processes that are accompanied by much more magmatic activity than previously expected. As most of the surface of the Arctic Ocean is covered by ice, deployment of ocean bottom seismometers is very difficult. In practice, there is a risk of losing the instruments and data; therefore, such deployment is rarely performed in Arctic areas. In this sense, the technologies of seismic investigations based on the deployment of seismic networks on ice floes are very promising. In the case of ice drift, such stations permanently change their locations, which allows larger areas to be covered by the network. Pioneering works on deploying seismic stations on ice floes were performed by the team led by Vera Schlindwein from the AWI Bremerhaven (Schlindwein et al.2007; Laderach & Schlindwein 2011). They successfully conducted several field campaigns in which they recorded a large number of microseismic events along the Gakkel Ridge and even found some events related to volcanic activity (Schlindwein & Riedel 2010). Furthermore, the local earthquake data collected by their floating networks were used to build seismic tomography models of the upper crust in the area of volcanic activity on the Gakkel Ridge (Korger & Schlindwein 2013). This relatively inexpensive technology appears to be very efficient for performing seismological studies. However, deploying stations on ice located above a thick water layer makes recording the seismic shear waves problematic. The S waves do not propagate through liquid water; however, they can be revealed by the recording of Sw phases corresponding to S waves propagating in the solid crust, which are converted into acoustic waves when passing through the sea bottom (Fig. 1). The availability of such waves is very important for improving the accuracy of local earthquake locations and performing high-quality tomography. Furthermore, the existence of both P and S phases for tomography makes it possible to build Vp, Vs and especially Vp/Vs models that provide much more information about petrophysical parameters of the crust than a single P-wave velocity model. Figure 1. View largeDownload slide Conceptual scheme for recording seismic waves using seismometers on floating ice. Triangles depict seismic stations, and the star represents an earthquake. Ray paths are depicted by coloured solid or dotted lines (for the P and S waves, respectively). The records of the corresponding seismic phases are schematically shown above. Figure 1. View largeDownload slide Conceptual scheme for recording seismic waves using seismometers on floating ice. Triangles depict seismic stations, and the star represents an earthquake. Ray paths are depicted by coloured solid or dotted lines (for the P and S waves, respectively). The records of the corresponding seismic phases are schematically shown above. Schlindwein et al. (2007) reported Sw phases from both regional and local events recorded by the network on the ice floe. However, these phases are detected in only a limited number of seismograms. This limitation explains why most of the source locations and the seismic tomography model by Korger & Schlindwein (2013) used only the P phases. Clearly, the existence of detectable converted Sw waves depends on several factors, such as epicentral distance, focal depth and focal mechanism. In addition, for some relationships of epicentral distance versus thickness of the water layer, the converted Sw wave might be hidden by multiple waves. In this study, we perform numerical modelling of wave fields to investigate the propagation of seismic waves from local earthquakes in the model corresponding to the realistic observation system on the ice floe (Fig. 1), as described by Schlindwein et al. (2007). The main purpose of this work is to prove that the S-wave phases can be robustly retrieved by such networks and can be used for different seismological studies. The numerical simulation of the seismic field generated by an earthquake is a classic mechanical problem, which is described in dozens of monographs and articles (e.g. Graves 1996; Pitarka et al.1998; Olsen et al.2006; Ben-Menahem & Singh 2012; Golubev et al.2013; Hanyga 2016 and many others). In exploration geophysics, there are many examples of numerical simulation of wave fields generated by airguns on the sea surface that emit acoustic waves in the water layer and seismic waves in the solid Earth (e.g. Van Vossen et al.2002). In particular, Favorskaya et al. (2016) and Favorskaya & Petrov (2016) modelled seismic wave propagation in a realistic set-up corresponding to an oil exploration survey in the Arctic using acoustic sources on the sea surface. The problem of generating acoustic waves in the water layer after strong earthquakes has previously been considered in several studies in the context of tsunami generation (e.g. Nakamura et al. 2011, 2012; Maeda & Furumura 2013). However, we are not aware of any study that considered the details of wave fields recorded on the sea surface generated by earthquakes with different focal mechanisms. Such modelling seems to us crucial for prompt interpretation of seismograms recorded by seismic networks deployed on ice floes in the Arctic. 2 MATHEMATICAL DESCRIPTION OF THE MODEL 2.1 Source description In this study, we describe an earthquake using the standard seismic source theory, which is described in Aki & Richards (1980). Consider the equation of motion: \begin{eqnarray*} \rho {(\bf x)}\frac{{{\partial ^2}{u_i}}}{{\partial {t^2}}} - \frac{{\partial {\sigma _{ik}}}}{{\partial {x_k}}} = - {M_{ij}}(t)\frac{{\partial \delta ({\bf x} - {{\bf x}_s})}}{{\partial {x_j}}}, \end{eqnarray*} (1) where |${\bf{u}}({\bf{x}})$| is the displacement vector in the space point |${\bf{x}}$|⁠, |$\rho (\bf x)$| is the density and |${\sigma _{ij}}$| are the components of the Cauchy stress tensor. The right-hand side of eq. (1) represents a point source with the seismic moment |${M_{ij}}(t)$|⁠. The strain tensor is defined as follows: \begin{eqnarray*} {\varepsilon _{ij}} = \frac{1}{2}\left( {\frac{{\partial {u_j}}}{{\partial {x_i}}} + \frac{{\partial {u_i}}}{{\partial {x_j}}}} \right). \end{eqnarray*} (2) The equation of motion (eq. 1) and eq. (2), with Hooke's law \begin{eqnarray*} {\sigma _{ij}} = {c_{ijkl}}{\varepsilon _{kl}}, \end{eqnarray*} (3) compose the system of elasticity equations. In the source point (see Fig. 2a), consider the displacement discontinuity |$[{\bf{u}}]$| over a small (compared with the excited wavelength) fault surface with the normal vector |$\bf{ v}$| and area |$d\Sigma $|⁠. The displacement discontinuity is averaged over the fault surface. The corresponding seismic moment is given by the following expression: \begin{eqnarray*} {M_{pq}} = [{u_i}]{v_j}{c_{ijpq}}d\sum \end{eqnarray*} (4) Figure 2. View largeDownload slide Model indications used in the description of the algorithm. (a) Model of a simple shear source. (b) 2-D standard staggered grid structure. (c) Representative elastic volume and power vectors corresponding to stress components. Figure 2. View largeDownload slide Model indications used in the description of the algorithm. (a) Model of a simple shear source. (b) 2-D standard staggered grid structure. (c) Representative elastic volume and power vectors corresponding to stress components. Such double-couple seismic sources are conventionally used for earthquake modelling (Aki & Richards 1980). The dip, strike and slip angles |$\varphi $|⁠, |$\delta $| and |$\lambda $|⁠, respectively, as shown in Fig. 2(a), are used to define the fault plane and displacement orientations. The displacement vector and the normal to the fault plane are defined by the following expressions: \begin{eqnarray*} [{u_1}] &=& {\bf{u}} (\cos \lambda \cos \phi + \cos \delta \sin \lambda \sin \phi ), \nonumber \\ {[{u_2}]} &=& {\bf{u}} (\cos \lambda \sin \phi - \cos \delta \sin \lambda \cos \phi ), \nonumber \\ {[{u_3}]} &=& - {\bf{u}} \sin \delta \sin \lambda , \nonumber \\ {v_1} &=& - \sin \delta \sin \phi , \nonumber \\ {v_2} &=& \sin \delta \cos \phi , \nonumber \\ {v_3} &=& - \cos \delta, \end{eqnarray*} (5) where |${\bf{u}}$| is the displacement vector length; the indices 1, 2 and 3 correspond to the x, y and z axes, respectively. 2.2 Numerical solution of the elastic equations The algorithm for numerical simulation of seismic waves in this study used the general principle of the standard staggered grid (SSTG) scheme initially developed for 2-D isotropic elastic equations (Virieux 1986) and then expanded to the 3-D case (Bohlen 2002). The following velocity–stress first-order representation of the elastic equations is used (considering the y = 0 plane): \begin{eqnarray*} \frac{{\partial {v_x}}}{{\partial t}} &=& \frac{1}{\rho }\left( {\frac{{\partial {\sigma _{xx}}}}{{\partial x}} + \frac{{\partial {\sigma _{xz}}}}{{\partial z}}} \right),\nonumber \\ \frac{{\partial {v_z}}}{{\partial t}} &=& \frac{1}{\rho }\left( {\frac{{\partial {\sigma _{zx}}}}{{\partial x}} + \frac{{\partial {\sigma _{zz}}}}{{\partial z}}} \right),\nonumber \\ \frac{{\partial {\sigma _{xx}}}}{{\partial t}} &=& (\lambda + 2\mu )\frac{{\partial {v_x}}}{{\partial x}} + \lambda \frac{{\partial {v_z}}}{{\partial z}},\nonumber \\ \frac{{\partial {\sigma _{zz}}}}{{\partial t}} &=& (\lambda + 2\mu )\frac{{\partial {v_z}}}{{\partial z}} + \lambda \frac{{\partial {v_x}}}{{\partial x}},\nonumber \\ \frac{{\partial {\sigma _{xz}}}}{{\partial t}} &=& \mu (\frac{{\partial {v_x}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial x}}), \end{eqnarray*} (6) where |${\bf{v}} = ({v_x},{v_z})$| is the velocity of the particle displacement. Regular spatial and time discretizations are used: |${x_i} = {x_0} + i\Delta x,{\rm{ }}{{\rm{z}}_i} = {z_0} + j\Delta z$|⁠, |${t_k} = k\Delta t$|⁠. In addition to the nodes with integer indices |$(i,j,k)$|⁠, the half-integer index nodes |$(i \pm \frac{1}{2},j \pm \frac{1}{2},k \pm \frac{1}{2})$| are used. The nodes of two staggered grids, which are shown in Fig. 2(b), correspond to different components of the velocity vector and stress tensor. The second-order central differences are used for spatial and time derivative approximations. The resulting finite-difference equations are as follows: \begin{eqnarray*} U_{i,j}^{k + \frac{1}{2}} &=& U_{i,j}^{k - \frac{1}{2}} + {B_{i,j}}\frac{{\Delta t}}{{\Delta x}}\left( {\Sigma _{i + \frac{1}{2},j}^k - \Sigma _{i - \frac{1}{2},j}^k} \right) + {B_{i,j}}\frac{{\Delta t}}{{\Delta z}}\left( {\Xi _{i + \frac{1}{2},j}^k - \Xi _{i - \frac{1}{2},j}^k} \right), \nonumber \\ V_{i + \frac{1}{2},j + \frac{1}{2}}^{k + \frac{1}{2}} &=& V_{i + \frac{1}{2},j + \frac{1}{2}}^{k - \frac{1}{2}} + {B_{i + \frac{1}{2},j + \frac{1}{2}}}\frac{{\Delta t}}{{\Delta x}}\left( {\Xi _{i + 1,j + \frac{1}{2}}^k - \Xi _{i,j + \frac{1}{2}}^k} \right) + {B_{i + \frac{1}{2},j + \frac{1}{2}}}\frac{{\Delta t}}{{\Delta z}}\left( {T_{i + \frac{1}{2},j + 1}^k - T_{i + \frac{1}{2},j}^k} \right),\nonumber \\ \Sigma _{i + \frac{1}{2},j}^{k + 1} &=& \Sigma _{i + \frac{1}{2},j}^k + {(L + 2M)_{i + \frac{1}{2},j}}\frac{{\Delta t}}{{\Delta x}}\left( {U_{i + 1,j}^{k + \frac{1}{2}} - U_{i,j}^{k + \frac{1}{2}}} \right) + {L_{i + \frac{1}{2},j}}\frac{{\Delta t}}{{\Delta z}}\left( {V_{i + \frac{1}{2},j + 1}^{k + \frac{1}{2}} - V_{i,j}^{k + \frac{1}{2}}} \right),\nonumber \\ T_{i + \frac{1}{2},j}^{k + 1} &=& T_{i + \frac{1}{2},j}^k + {(L + 2M)_{i + \frac{1}{2},j}}\frac{{\Delta t}}{{\Delta z}}\left( {V_{i,j + 1}^{k + \frac{1}{2}} - V_{i,j}^{k + \frac{1}{2}}} \right) + {L_{i + \frac{1}{2},j}}\frac{{\Delta t}}{{\Delta x}}\left( {U_{i + 1,j}^{k + \frac{1}{2}} - U_{i,j}^{k + \frac{1}{2}}} \right),\nonumber \\ \Xi _{i,j + \frac{1}{2}}^{k + 1} &=& \Xi _{i,j + \frac{1}{2}}^k + {M_{i,j + \frac{1}{2}}}\frac{{\Delta t}}{{\Delta z}}\left( {U_{i,j + 1}^{k + \frac{1}{2}} - U_{i,j}^{k + \frac{1}{2}}} \right) + {M_{i,j + \frac{1}{2}}}\frac{{\Delta t}}{{\Delta x}}\left( {V_{i + 1,j}^{k + \frac{1}{2}} - V_{i,j}^{k + \frac{1}{2}}} \right), \end{eqnarray*} (7) where |$(U,V) = ( {{u_x},{u_z}} )$| are the values of the components of the particle velocity vector and |$( {\Sigma ,\Xi ,T} ) = ( {{\sigma _{xx}},{\sigma _{zz}},{\sigma _{xz}}} )$| are the values of the components of the stress tensor in the corresponding spatial-time grid nodes, B is the inverse of the density and |$L,M$| correspond to the Lame parameters |$\lambda ,\mu $|⁠. As usual, the discretization steps along spatial variables are the same: |$\Delta x = \Delta z = h$|⁠. For the considered second-order scheme, a good choice is 12 points per minimum wavelength: \begin{eqnarray*} h \le \frac{{V_S^{\min }}}{{12 \cdot {\omega _{\max }}}}. \end{eqnarray*} (8) The time discretization step |$\Delta t$| is determined based on the Courant–Friedrichs–Lewy condition. For the SSTG scheme, |$\Delta t$| can be defined as \begin{eqnarray*} \Delta t \le \frac{{0,6 \cdot h}}{{V_P^{\max }}}. \end{eqnarray*} (9) The standard technique to suppress artificial boundary reflections is perfectly matched layers (Kristek et al.2009). One can also introduce the free surface, namely the sea surface in the considered case (sea ice not being taken into account). The SSTG scheme is suitable for numerical seismic wave propagation in a fluid–solid configuration (Van Vossen et al.2002). A seafloor is placed in a layer of nodes with the half-integer indices |$j + 1/2$|⁠, which correspond to the horizontal component of particle velocity |${u_z}$| and the |${\sigma _{xz}}$| stress component (see Fig. 2b). The seismic moment tensor is symmetric. A couple of forces balances the symmetrical couple of forces with the non-zero moment (⁠|${M_{ij}} = {M_{ji}}$|⁠). Consider the elementary volume, which is in a state of equilibrium. The resulting force and momentum induced by the stress components |${\sigma _{ij}} = {\sigma _{ji}}$|(shown in Fig. 2c) are zero. The inclusion of a point source, which is defined by the seismic moment tensor, inside the elementary volume is equivalent to corresponding additives to the components of the stress tensor (see Fig. 2c). This process is how the dip–strike–slip point source can be implemented in the SSTG scheme (Coutant et al.1995). 3 NUMERICAL EXPERIMENTS 3.1 Seismic waves in 2-D Prior to simulating the 3-D wave fields, we perform 2-D modelling of wave propagation in the model shown in Fig. 3(a). We consider a model of solid rock half-space overlain by a 3-km-thick water layer. In water, we define a constant P-wave velocity (Vp) of 1.5 km s−1 and a density of 1000 kg m−3. Shear waves cannot propagate within this layer. In the lower solid half-space, we define a constant gradient for the P-wave velocity that linearly increases from 4 km s−1 at a depth of 3 km below sea level (bsl; sea bottom) to 6 km s−1 at a depth of 15 km. In the solid half-space, the S-wave velocity (Vs) is derived from Vp based on a constant Vp/Vs ratio equal to 1.7. The density values in the solid rock layer are determined as half of the value of the P-wave velocity. To avoid reflections from the bottom and side boundaries, we define absorbing layers with extremely high attenuation at the boundaries: z = 23 km, x = 0 km and x = 60 km. The source point is located at z = 8 km (or 5 km below the sea bottom) and x = 30 km. We present simulation results for a source corresponding to normal faulting with dip angle of the slip plane of 45 degrees. To model the source pulse, we use the Ricker wavelet with 4.5 Hz dominant frequency, which corresponds to a typical frequency of seismic signal from a local earthquake located at some tens of kilometres of epicentral distance. This frequency is higher than that used in the main 3-D model because the main purpose of this 2-D modelling is to identify the major seismic phases that would be less clear in the case of a lower frequency source. Figure 3. View largeDownload slide Two models with different 1-D P-wave velocity distributions used for simulations of seismic fields. The white dots depict the source points. The receivers are located on the sea surface. Figure 3. View largeDownload slide Two models with different 1-D P-wave velocity distributions used for simulations of seismic fields. The white dots depict the source points. The receivers are located on the sea surface. The snapshots of the resulting 2-D wavefield with indications of the main phases are shown in Fig. 4. The animation of wave propagation for this model is presented in the Supporting Information. We can observe several waves that correspond to different combinations of the P and S waves in the solid medium, acoustic waves in the water layer (denoted as w) and multiples. The schematic ray geometries for the major observed phases are summarized in Fig. 5. In previous studies, the water segment of converted waves was indicated with ‘P’ (e.g. ‘SP’ in Schlindein et al.2007). However, in this study, we denote the acoustic wave in the water layer with the symbol ‘w’ (e.g. ‘Pw’ or ‘Sw’) to distinguish it from cases of multiple reflections of the P wave from the sea bottom boundary in the solid Earth and to avoid other confusions. In Fig. 4, we observe clear Pw and Sw phases corresponding to the S and P waves in the solid half-space, respectively, converted to acoustic waves in the water layer. In addition to these waves, we observe a clear phase indicated as SPw that propagates faster than the Sw phase. Looking at the wave field in 2-D, we propose that the SPw phase might correspond to the S- to P-wave conversion after reflection in the solid half-space, refraction and conversion to the acoustic wave in the water layer. The same phases are observed in the 3-D cases described below. Figure 4. View largeDownload slide Snapshot of the 2-D modelling results in the model shown in Fig. 3(a) at 4.8 and 7.5 s after the source initiation. The source is a normal fault with a dip-slip angle of 45 degrees. Several phases of the major waves are indicated. The seafloor is shown by a black line. Figure 4. View largeDownload slide Snapshot of the 2-D modelling results in the model shown in Fig. 3(a) at 4.8 and 7.5 s after the source initiation. The source is a normal fault with a dip-slip angle of 45 degrees. Several phases of the major waves are indicated. The seafloor is shown by a black line. Figure 5. View largeDownload slide Schematic representation of the ray geometry of the main phases identified in modelling results. The solid lines depict the P waves, and the dotted lines show the S waves. The grey layer depicts the sea water. The star is the source, and triangles are the receivers. Figure 5. View largeDownload slide Schematic representation of the ray geometry of the main phases identified in modelling results. The solid lines depict the P waves, and the dotted lines show the S waves. The grey layer depicts the sea water. The star is the source, and triangles are the receivers. 3.2 Seismic wave simulations in the 3-D model We simulate the 3-D seismic waves generated by earthquakes as represented by a simple double-couple model (Aki & Richards 1980). We considered a variety of different source mechanisms, different cases of source locations and different velocity distributions. In addition, we made numerical simulations for a model with varied bathymetry. First of all, we perform numerical modelling for three types of sources with different focal mechanisms: strike-slip, vertical displacement and normal faulting for the same velocity model as considered in the case of the 2-D modelling shown in Fig. 3(a). The source pulse is described by the Ricker wavelet with a dominant frequency of 1 Hz. Here, we present the vertical displacement recorded at the sea surface along one circular line located 20 km from the epicentre and four radial profiles, as shown in Fig. 6. The results of simulations for three types of source mechanisms are shown in Figs 7–9. In the synthetic seismograms, we highlight the arrivals of major recorded waves. Figure 6. View largeDownload slide Source position (red dot) and receiver profiles (with the direction of receiver numbering). Figure 6. View largeDownload slide Source position (red dot) and receiver profiles (with the direction of receiver numbering). Figure 7. View largeDownload slide Results of the 3-D seismic wave simulations for the strike-slip source mechanism for the model set-up shown in Fig. 3(a). Vertical displacements at the sea surface are presented along one circular line and four profiles with locations shown in Fig. 6. The major primary phases of Pw, Sw and SPw waves are highlighted with red lines. Figure 7. View largeDownload slide Results of the 3-D seismic wave simulations for the strike-slip source mechanism for the model set-up shown in Fig. 3(a). Vertical displacements at the sea surface are presented along one circular line and four profiles with locations shown in Fig. 6. The major primary phases of Pw, Sw and SPw waves are highlighted with red lines. Figure 8. View largeDownload slide Same as Fig. 7 but for the source mechanism of vertical displacement. Figure 8. View largeDownload slide Same as Fig. 7 but for the source mechanism of vertical displacement. Figure 9. View largeDownload slide Same as Fig. 7 but for the source mechanism of normal faulting. Figure 9. View largeDownload slide Same as Fig. 7 but for the source mechanism of normal faulting. The results of the wavefield simulation for the case of the strike-slip event for the velocity model shown in Fig. 3(a) are presented in Fig. 7. It can be seen that horizontal particle motions in the source lead to the origin of both Pw and Sw waves of similar amplitudes. The azimuthal distribution of polarity of both types of waves is determined by the compressional and extensional quadrants of the source mechanism. In the directions of 0, 90, 180 and 270 degrees, the amplitudes of all recorded Pw and Sw waves are zero. At distances larger than 20 km, we can observe that the Sw wave interferes with the multiple Pw2 wave, which may make identification of the correct arrival time of the Sw wave at a single station difficult. In addition to the Pw and Sw phases, we observe another phase, which is identified as SPw and has the same apparent velocity as the Pw wave. We propose that this phase is due to reflection of the S wave from the sea bottom and conversion to the P wave (Fig. 5g). This phase becomes visible from an epicentral distance of 15 km and then arrives at earlier times than the Sw wave. It was unexpected to observe that the amplitude of the SPw phase at a distance of 20 km appeared to be stronger than that of the Sw phase. In this case, explaining why we identify a very clear SPw phase but do not observe the PPw phase with similar ray geometry and composed of only P waves is not obvious. In the later recording times after arrivals of the Pw, Sw and SPw phases, we clearly observe a series of multiple waves related to these three types of primary waves. In a circular line at 20 km from the epicentre, the Pwn and Swn multiples remain visible throughout the recording time interval with very little decay, whereas the SPwn phase attenuates after the first reflection and is almost invisible at later times. Note that every multiple wave has the opposite polarity of particle motions (compression or extension) with respect to the previous wave. The case of a source with vertical displacement (Fig. 8) in the velocity model shown in Fig. 3(a) results in considerably different wavefields than in the previous case. The Pw wave is robustly observed near the epicentre but strongly decays at distances of more than 10 km. At the circular line with a 20 km radius, the Pw wave is much weaker than the Sw wave. Correspondingly, Pwn multiples can hardly be identified, whereas the Swn multiples are very clear and are the only phases observed at later times. The Sw-wave polarity depends on the azimuth with respect to the displacement in the source. The negative phase is observed from 0 to 180 degrees, where the ground moves downwards, and the positive from 180–360 degrees of azimuth on the side of upward displacement. The SPw phase appears to be much weaker than in the case of the strike-slip source mechanism. This difference can be explained by the dominant vertical particle displacement, which contributes very little to the reflected Pwave. Fig. 9 shows the wavefield emitted by a source with the normal faulting mechanism (slip at a dip angle of 45 degrees) in the velocity model shown in Fig. 3(a). This source generates a Pw wave with only positive polarity (compression), which is strongest at azimuths of 90 and 270 degrees and has zero amplitudes at 0 and 180 degrees. In the circular line, the Pwn multiples have alternating polarities and gradually attenuate with increasing n. The SPw phase, which arrives next after the Pw phase, has the maximum amplitude. The first multiple SPw2 is much weaker, and the other multiples of this wave are almost not recognizable. The Sw phase is initially weaker than the SPw phase, but the multiples of Swn almost do not decay, and they become the strongest waves at later times. Note that both Swn and SPwn have significant amplitudes of anomalies at the azimuths of 0 and 180 degrees, whereas the primary Pw wave has zero amplitudes at these azimuths. To assess the effects of the source depth and velocity distributions, we have created another model, shown in Fig. 3b. In this case, the source was located at 11 km depth bsl or 8 km below sea bottom. The P-wave velocity in this model had the constant gradient and varied from 4.5 km s−1 at 3 km bsl and 8 km s−1 at 15 km bsl. The shear velocity and density were computed using the same relationships as in the first model. Here, we used the same frequency of 1 Hz as in the previous cases. In Fig. 10, we present the result of the wavefield modelling based on the normal faulting mechanism. Here, we show the wavefield in one straight line x2 oriented S–N and one circular line at the radius of 20 km of epicentral distance. Comparing with the corresponding sections in Fig. 9 reveals some important differences. In Fig. 10, we observe only Pw and Sw waves and their multiples. No SPw wave is observed here. Likely, this wave would appear at larger distances, which are outside our computing area. Another difference of this case is that the multiple waves at 20 km of epicentral distance appear to be much weaker compared to the case shown in Fig. 9. In the linear profile, we see that in the case with higher velocity and a deeper source, the clearest multiple waves are detected at 15 km of epicentral distance. Figure 10. View largeDownload slide Results of the 3-D seismic wave simulations for the source mechanism of normal faulting for a higher velocity gradient model, shown in Fig. 3b. Vertical displacements at the sea surface are presented along the x2 profile and the circular line (see Fig. 6). The major primary phases of Pw and Sw waves are highlighted with red lines. Figure 10. View largeDownload slide Results of the 3-D seismic wave simulations for the source mechanism of normal faulting for a higher velocity gradient model, shown in Fig. 3b. Vertical displacements at the sea surface are presented along the x2 profile and the circular line (see Fig. 6). The major primary phases of Pw and Sw waves are highlighted with red lines. To make the simulations more realistic, we have created a model with varied sea bottom topography (Fig. 11). The location of the source and velocity gradient was the same as in the case of the model shown in Fig. 2(a) for which most of the results were computed. The topography was created artificially and had the deviations from 2.7 km to 3.3 km depth bsl. In this case, we used a higher frequency of 2 Hz for the source generating the seismic signal. The source mechanism in this case was normal faulting, the same as in the case presented in Fig. 9. The results of the wavefield simulation are shown in Fig. 12 in one S–N oriented profile and in a circular line at the epicentral distance of 20 km. It can be seen that the presence of a higher frequency source and varied bathymetry does change considerably the shapes of the main waves (Pw, SPw and Sw). In the circular line, we observe some time delays in the azimuthal direction caused by the bathymetry variations. We can also observe in Fig. 12 that the multiple waves are much less clear. After two or three reflections they become hardly visible, whereas in Fig. 9, we observe the multiples in the whole range of travel times. Less clear multiple waves in the case of the model with varied bathymetry might be caused by the scattering effect of the topography on the converted and reflected waves in the water layer. Figure 11. View largeDownload slide Model set-up of the case with varied bathymetry. (a) Shaded relief of the seafloor relief, source position (red dot) and receiver profiles; (b) and (c) 2-D velocity distributions in two vertical sections located along the x1 and x2 profiles; (d) and (e) exaggerated bathymetry along the x1 and x2 profiles. Figure 11. View largeDownload slide Model set-up of the case with varied bathymetry. (a) Shaded relief of the seafloor relief, source position (red dot) and receiver profiles; (b) and (c) 2-D velocity distributions in two vertical sections located along the x1 and x2 profiles; (d) and (e) exaggerated bathymetry along the x1 and x2 profiles. Figure 12. View largeDownload slide Results of the 3-D seismic wave simulations for the source mechanism of normal faulting for the model shown in Fig. 11. Vertical displacements at the sea surface are presented along the x2 profile and circular line (see Fig. 11a). The major primary phases of Pw, SPw and Sw waves are highlighted with red lines. Figure 12. View largeDownload slide Results of the 3-D seismic wave simulations for the source mechanism of normal faulting for the model shown in Fig. 11. Vertical displacements at the sea surface are presented along the x2 profile and circular line (see Fig. 11a). The major primary phases of Pw, SPw and Sw waves are highlighted with red lines. 4 DISCUSSION AND CONCLUSIONS The results of numerical modelling of the wavefields generated by different types of seismic sources produce shear waves that can be recorded at the sea surface as the Sw phase. These data are very important for obtaining information about the S waves in the solid rock below the sea layer that can be used, for example, to perform seismic tomography inversion for the S-wave velocity distributions. However, there is a risk of confusion with the SPw phase, which may arrive earlier and appear to be stronger than the Sw phase. Misidentification of this phase may cause regular bias in the S-wave travel times and may lead to unrealistically high values of the S-wave velocities in the crust. Another unexpected issue is that the Pw waves are much more sensitive to the type of source and the azimuthal coverage than the Sw waves. For example, in the case of vertical displacement (Fig. 8), the first arrival of the P wave is almost invisible. In practice, when analysing such seismograms, one may pick the Sw phase as the first arrival and associate it with the Pw phase. The same confusion may occur for the normal slip event (Fig. 9) when the station is close to the azimuths of 0 and 180 degrees. The only robust identification of the Pw phase as the first arrival can be performed for the strike-slip source mechanism (Fig. 7), for which this phase has the maximum amplitudes at all azimuths. A large problem of time picking might be related to the interference of several wave fronts, such as those of the Sw and multiple Pw2. Our results have shown that in our model, these waves interfere at distances larger than 20 km. Obviously, this distance depends largely on the velocity model, source depth and thickness of the water layer, as was demonstrated in tests shown in Figs 10 and 12. When processing experimental data, it is necessary to take this effect into account to avoid misidentification of wave phases. To perform approximate estimates for the interference zone, simple kinematic modelling would be sufficient. We confess that the frequencies of 1 and 2 Hz (which were limited by our computing facilities) are too low to adequately simulate the signal of the local earthquake with the expected major frequency range between 3 and 6 Hz. At the same time, the comparison of the 3-D results for 1 and 2 Hz sources with the 2-D simulations at 4.5 Hz frequency shows that the arrivals of the major waves in these cases are similar and provide identical kinematic characteristics. No additional waves are visible when using the high-frequency approximation. Therefore, we claim that the low-frequency approximation works adequately to reveal the main features of seismic wave propagation in such kinds of models. Another important conclusion that follows from comparison of the modelling results in Figs 9 and 12 is that the bathymetry variations, though rather smooth, lead to considerable scattering and strong decay of the multiple waves. This should be taken into account when performing experiments in flat of tectonically disturbed areas. We expect that the results of our modelling will be important during the analysis of seismic data recorded by seismic stations deployed on ice floes, as was done by Schlindwein et al. (2007, 2015). In particular, our results will help to identify correct arrival times of the P and especially S waves from local seismicity and avoid misinterpretations related to other phases. SUPPORTING INFORMATION wave.avi Please note: Oxford University Press are 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. ACKNOWLEDGEMENTS We are grateful for the comments by Dr Vera Schlindwein and two other anonymous reviewers that improved the manuscript. AS and AY are supported by the Russian Foundation for Basic Research Grant # 18-35-20030. IK is supported by the Russian Science Foundation Grant # 18-17-00095. REFERENCES Aki K. , Richards P.G. , 1980 . Quantitative Seismology: Theory and Methods , Vol. 1 , WH Freeman & Co . Ben-Menahem A. , Singh S.J. , 2012 . Seismic Waves and Sources , Springer Science & Business Media . Bohlen T. , 2002 . Parallel 3-D viscoelastic finite difference seismic modelling , Comput. Geosci. , 28 ( 8 ), 887 – 899 . Google Scholar Crossref Search ADS Coutant O. , Virieux J. , Zollo A. , 1995 . Numerical source implementation in a 2D finite difference scheme for wave propagation , Bull. seism. Soc. Am. , 85 ( 5 ), 1507 – 1512 . Favorskaya A. , Petrov I. , Khokhlov N. , 2016 . Numerical modeling of wave processes during shelf seismic exploration , Procedia Comput. Sci. , 96 , 920 – 929 .. https://doi.org/10.1016/j.procs.2016.08.271 Google Scholar Crossref Search ADS Favorskaya A.V. , Petrov I.B. , 2016 . Wave responses from oil reservoirs in the Arctic shelf zone , Dokl. Earth Sci. , 466 ( 2 ), 214 – 217 .. https://doi.org/10.1134/S1028334X16020185 Google Scholar Crossref Search ADS Gaina C. , Medvedev S. , Torsvik T.H. , Koulakov I. , Werner S.C. , 2014 . 4D Arctic: a glimpse into the structure and evolution of the Arctic in the light of new geophysical maps, plate tectonics and tomographic models , Surv. Geophys. , 35 ( 5 ), 1095 – 1122 .. https://doi.org/10.1007/s10712-013-9254-y Google Scholar Crossref Search ADS PubMed Golubev V.I. , Petrov I.B. , Khokhlov N.I. , 2013 . Numerical simulation of seismic activity by the grid-characteristic method , Comput. Math. Math. Phys. , 53 ( 10 ), 1523 – 1533 .. https://doi.org/10.1134/S0965542513100060 Google Scholar Crossref Search ADS Graves R.W. , 1996 . Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences , Bull. seism. Soc. Am. , 86 ( 4 ), 1091 – 1106 . Hanyga A. , ed, 2016 . Seismic Wave Propagation in the Earth , Elsevier . Jakovlev A. , Bushenkova N.A. , Koulakov I.Yu. , Dobretsov N.L. , 2012 . Structure of the upper mantle in the circum-Arctic region from regional seismic tomography , Russ. Geol. Geophys. , 53 ( 10 ), 963 – 971 .. https://doi.org/10.1016/j.rgg.2012.08.001 Google Scholar Crossref Search ADS Korger E.I.M. , Schlindwein V. , 2013 . Seismicity and structure of the 85 E volcanic complex at the ultraslow spreading Gakkel Ridge from local earthquake tomography , Geophys. J. Int. , 196 ( 1 ), 539 – 551 .. https://doi.org/10.1093/gji/ggt390 Google Scholar Crossref Search ADS Koulakov I. , Gaina C. , Dobretsov N.L. , Vasilevsky A.N. , Bushenkova N.A. , 2013 . Plate reconstructions in the Arctic region based on joint analysis of gravity, magnetic, and seismic anomalies , Russ. Geol. Geophys. , 54 ( 8 ), 743 – 757 .. https://doi.org/10.1016/j.rgg.2013.07.007 Google Scholar Crossref Search ADS Kristek J. , Moczo P. , Galis M. , 2009 . A brief summary of some PML formulations and discretizations for the velocity-stress equation of seismic motion , Stud. Geophys. Geod. , 53 ( 4 ), 459 – 474 .. https://doi.org/10.1007/s11200-009-0034-6 Google Scholar Crossref Search ADS Laderach C. , Schlindwein V. , 2011 . Seismic arrays on drifting ice floes: experiences from four deployments in the Arctic Ocean , Seismol. Res. Lett. , 82 ( 4 ), 494 – 503 .. https://doi.org/10.1785/gssrl.82.4.494 Google Scholar Crossref Search ADS Maeda T. , Furumura T. , 2013 . FDM simulation of seismic waves, ocean acoustic waves, and tsunamis based on tsunami-coupled equations of motion , Pure appl. Geophys. , 170 ( 1–2 ), 109 – 127 .. https://doi.org/10.1007/s00024-011-0430-z Google Scholar Crossref Search ADS Muller C. , Jokat W. , 2000 . Seismic evidence for volcanic activity discovered in central Arctic , EOS, Trans. Am. geophys. Un. , 81 ( 24 ), 265 – 269 .. https://doi.org/10.1029/00EO00186 Google Scholar Crossref Search ADS Nakamura T. , Takenaka H. , Okamoto T. , Kaneda Y. , 2011 . A study of the finite difference solution for 3D seismic wavefields near a fluid-solid interface , Zisin, J. Seismol. Soc. Japan , 63 ( 3 ), 187 – 194 .. https://doi.org/10.4294/zisin.63.189 Nakamura T. , Takenaka H. , Okamoto T. , Kaneda Y. , 2012 . FDM simulation of seismic‐wave propagation for an aftershock of the 2009 Suruga Bay earthquake: effects of ocean‐bottom topography and seawater layer , Bull. seism. Soc. Am. , 102 ( 6 ), 2420 – 2435 .. https://doi.org/10.1785/0120110356 Google Scholar Crossref Search ADS Olsen K.B. , et al. . , 2006 . Strong shaking in Los Angeles expected from southern San Andreas earthquake , Geophys. Res. Lett. , 33 ( 7 ), doi: https://doi.org/10.1029/2005GL025472 . Pitarka A. , Irikura K. , Iwata T. , Sekiguchi H. , 1998 . Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogo-ken Nanbu (Kobe), Japan, earthquake , Bull. seism. Soc. Am. , 88 ( 2 ), 428 – 440 . Schaefer A.J. , Lebedev S. , 2013 . Global shear speed structure of the upper mantle and transition zone , Geophysical Journal International , 194 ( 1 ), 417 – 449 . Google Scholar Crossref Search ADS Schlindwein V. , Demuth A. , Korger E. , Laderach C. , Schmid F. , 2015 . Seismicity of the Arctic mid-ocean ridge system , Polar Sci. , 9 ( 1 ), 146 – 157 .. https://doi.org/10.1016/j.polar.2014.10.001 Google Scholar Crossref Search ADS Schlindwein V. , Muller C. , Jokat W. , 2005 . Seismoacoustic evidence for volcanic activity on the ultraslow-spreading Gakkel Ridge, Arctic Ocean , Geophys. Res. Lett. , 32 ( 18 ), doi: https://doi.org/10.1029/2005GL023767 . Schlindwein V. , Müller C. , Jokat W. , 2007 . Microseismicity of the ultraslow-spreading Gakkel ridge, Arctic Ocean: a pilot study , Geophys. J. Int. , 169 ( 1 ), 100 – 112 .. https://doi.org/10.1111/j.1365-246X.2006.03308.x Google Scholar Crossref Search ADS Schlindwein V. , Riedel C. , 2010 . Location and source mechanism of sound signals at Gakkel Ridge, Arctic Ocean: submarine Strombolian activity in the 1999–2001 volcanic episode , Geochem. Geophys. Geosyst. , 11 ( 1 ), doi: https://doi.org/10.1029/2009gc002706 . Van Vossen R , Robertsson J.O.A. , Chapman C.H. , 2002 . Finite-difference modeling of wave propagation in a fluid-solid configuration , Geophysics , 67 ( 2 ), 618 – 624 .. https://doi.org/10.1190/1.1468623 Google Scholar Crossref Search ADS Virieux J. , 1986 . P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method , Geophysics , 51 ( 4 ), 889 – 901 .. https://doi.org/10.1190/1.1442147 Google Scholar Crossref Search ADS This content is only available as a PDF. © The Author(s) 2019. 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 - Numerical modelling of seismic waves from earthquakes recorded by a network on ice floes
JO - Geophysical Journal International
DO - 10.1093/gji/ggz148
DA - 2019-07-01
UR - https://www.deepdyve.com/lp/oxford-university-press/numerical-modelling-of-seismic-waves-from-earthquakes-recorded-by-a-2s9kimWguB
SP - 74
VL - 218
IS - 1
DP - DeepDyve
ER -