# Global SH-wave propagation in a 2D whole Moon model using the parallel hybrid PSM/FDM method

Xianghua Jiang, Yanbin Wang, Yanfang Qin, Hiroshi Takenaka

► To cite this version:

Xianghua Jiang, Yanbin Wang, Yanfang Qin, Hiroshi Takenaka. Global SH-wave propagation in a 2D whole Moon model using the parallel hybrid PSM/FDM method. Earthquake Science, 2015, 28, pp.163-174. 10.1007/s11589-015-0121-4. insu-03579971

HAL Id: insu-03579971

https://insu.hal.science/insu-03579971

Submitted on 18 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH PAPER

# Global SH-wave propagation in a 2D whole Moon model using the parallel hybrid PSM/FDM method

Xianghua Jiang · Yanbin Wang · Yanfang Qin · Hiroshi Takenaka

Received: 24 March 2015 / Accepted: 5 May 2015 / Published online: 6 June 2015
© The Author(s) 2015. This article is published with open access at Springerlink.com

**Abstract** We present numerical modeling of SH-wave propagation for the recently proposed whole Moon model and try to improve our understanding of lunar seismic wave propagation. We use a hybrid PSM/FDM method on staggered grids to solve the wave equations and implement the calculation on a parallel PC cluster to improve the computing efficiency. Features of global SH-wave propagation are firstly discussed for a 100-km shallow and 900-km deep moonquakes, respectively. Effects of frequency range and lateral variation of crust thickness are then investigated with various models. Our synthetic waveforms are finally compared with observed Apollo data to show the features of wave propagation that were produced by our model and those not reproduced by our models. Our numerical modeling show that the low-velocity upper crust plays significant role in the development of reverberating wave trains. Increasing frequency enhances the strength and duration of the reverberations. Surface multiples dominate wavefields for shallow event. Core–mantle reflections can be clearly identified for deep event at low frequency. The layered whole Moon model and the low-velocity upper crust produce the reverberating wave trains following each phases consistent with

observation. However, more realistic Moon model should be considered in order to explain the strong and slow decay scattering between various phases shown on observation data.

**Keywords** Whole Moon model · Seismic wavefield · SH-wave propagation · Hybrid method · Parallel computing

## 1 Introduction

The Apollo Passive Seismic Experiment consisted of a network of four seismometers deployed on the near side of the Moon during the Apollo missions between 1969 and 1972, which operated continuously until 1977 and recorded more than 12,000 moonquakes (Lognonné 2005). These data provided the basis for seismic studies of the lunar interior. Analyses of the lunar seismic data revealed the one-dimensional (e.g., Nakamura 1983; Khan and Mosegaard 2002; Lognonné et al. 2003; Gagnepain-Beyneix et al. 2006), lateral heterogeneous (e.g., Chenet et al. 2006; Zhao et al. 2008) lunar seismic velocity structure in the crust and mantle, the possible mechanism of moonquakes (e.g., Nakamura 1978; Nakamura et al. 1979; Koyama and Nakamura 1980) and their implications on the composition and physical properties of the lunar interior (e.g., Khan et al. 2006, 2007). Recent reanalysis of lunar seismograms (Weber et al. 2011; Garcia et al. 2011) determined the velocity structure of the lunar core and core–mantle transition zone. Based on both seismological and geodesic data, Garcia et al. (2011) constructed a very preliminary reference Moon model (VPREMOON).

Lunar seismograms differ greatly from the typical terrestrial seismic signals. They are characterized with reverberating wave trains of very long duration and a slow

<footnote>X. Jiang · Y. Wang (✉)
Department of Geophysics, School of Earth and Space Sciences,
Peking University, Beijing 100871, China
e-mail: ybwang@pku.edu.cn</footnote>

<footnote>Y. Qin
Equipe de Géosciences Marines, Institut de Physique du Globe
de Paris, 4 Place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote>H. Takenaka
Department of Earth Sciences, Faculty of Science, Okayama
University, 3-1-1 Tsushima-Naka, Kita-ku, Okayama 700-8530,
Japan</footnote>
decay of amplitude. Except for the direct P wave and S wave for some events, various secondary seismic phases and surface waves are difficult to identify from the intense scattering and reverberating waveforms. Therefore, the process of seismic wave propagation inside the whole Moon is still not well understood from the observed lunar seismograms. However, numerical modeling of seismic wave propagation in the whole Moon model had been conducted to improve our understanding of lunar seismic wavefield. Dainty et al. (1974) and Nakamura (1977) applied the diffusion theory to a Moon model with randomly distributed scatters near the surface and closely reproduced the envelope of typical lunar seismograms from impacts and moonquakes. Their results explained some features of the scattering and reverberating lunar seismic signals of long duration and suggested the existence of a shallow scattering zone near the surface. Lawrence and Johnson (2010) and Blanchette-Guertin et al. (2012) applied a seismic phonon method to calculate global synthetic lunar seismograms and investigated the effects of seismic scattering in a highly heterogeneous regolith layer on high-frequency lunar signals. Wang et al. (2013) recently performed P-SV seismic wave propagation modeling for a 2D lateral heterogeneous whole Moon model with a pseudospectral (PSM) and finite difference (FDM) hybrid method. They calculated P-SV body wave propagation in the whole Moon radiated from both shallow and deep moonquakes, compared synthetics with observed Apollo seismic data, and found that the near-surface low-velocity layer contributes significantly to the development of waveform reverberation.

In an isotropic 2D whole Moon model, the horizontally polarized SH-waves are recoded in the transverse component. SH-waves are decoupled from P and SV waves; therefore, the SH-wavefield contains only the direct SH-waves and their interactions with interfaces are much simpler than the P and SV wavefields recorded in the radial and vertical components. For a similar 2D whole-Earth model, several numerical methods had been applied to simulate global SH-wave propagation. Cummins et al. (1994) proposed a direct solution method to calculate global SH-wavefield for a spherically symmetric model. Igel and Weber (1995) calculated 2D global SH-wave propagation with a high-order FDM for a 2D axisymmetric whole-Earth model. Igel and Gudmundsson (1997) investigated the frequency-dependent effects of S and SS waves with the similar FDM method. Thorne et al. (2007) studied lateral heterogeneity in the lower mantle beneath the Cocos Plate with a 2D global SH-wavefield simulation based on FDM. Jahnke et al. (2008) performed 2D global SH-wave modeling on a parallel computer. Toyokuni et al. (2005) presented a new quasi-spherical FDM approach for 2D global wavefield simulation. Wang and Takenaka (2011) applied a PSM method to calculate 2D global SH-wave propagation

with arbitrary lateral heterogeneities. Wang et al. (2014) developed a parallel PSM/FDM hybrid method for calculating global SH-wave propagation on PC cluster. The SH-wave component observation data had often been used to study the localized shear-velocity heterogeneities in the Earth’s interior for its simple waveforms (e.g., Wen 2002; Wang and Chen 2009).

In this paper, we attempted to simulate global SH-wave propagation in a 2D whole Moon model proposed by Weber et al. (2011) and Garcia et al. (2011). We apply a hybrid PSM/FDM method on staggered grids to solve the 2D SH-wave equations. We implemented parallel computing on a PC cluster for the hybrid method to improve computing efficiency. The 2D whole Moon model is decomposed into sub-domains in the radial direction and assigned to different processors of a PC cluster. We applied the parallel hybrid method to a lateral homogeneous whole Moon model for both deep and shallow moonquakes to show the process of SH-wave propagation in the whole Moon. Comparisons between the synthetics and observed Apollo seismograms from typical moonquake events were performed to illustrate wavefield features that our simulations reproduced and those that were not achieved. We also investigated the effects of lateral varying crustal thickness and frequency content on synthetic waveforms with various numerical models.

## 2 SH-wave equations and the parallel PSM/FDM scheme

### 2.1 Elastodynamic wave equations

Considering a cylindrical coordinate system $(r, \theta, z)$ and assuming invariance for all field variables in $z$, we get a 2D cylindrical coordinate system $(r, \theta)$. For an isotropic and elastic medium, SH-wave propagation in such 2D coordinates can be formulated by the following equations of momentum conservation in velocity-stress form:

$$ \rho \frac{\partial v_z}{\partial t} = \frac{\partial \sigma_{zr}}{\partial r} + \frac{1}{r} \left( \frac{\partial \sigma_{z\theta}}{\partial \theta} + \sigma_{zr} \right) + f_z, \quad (1) $$

where $v_z$ is the velocity component in anti-plane direction, $\rho$ is the mass density, $\sigma_{zr}$, $\sigma_{z\theta}$ are components of stress tensor, and $f_z$ is the body force component in $z$ direction. The constitutive relations between stress and velocity are given by:

$$ \begin{aligned} \frac{\partial \sigma_{zr}}{\partial t} &= \mu \frac{\partial v_z}{\partial r} \\ \frac{\partial \sigma_{z\theta}}{\partial t} &= \frac{\mu}{r} \frac{\partial v_z}{\partial \theta}, \end{aligned} \quad (2) $$

where $\mu$ is the shear modulus of the medium.
## 2.2 The parallel hybrid PSM/FDM scheme

We calculate SH-wave propagation in a 2D whole Moon model, which corresponds to a cross section of the 3D Moon cutting through a great circle. The 2D model is defined from the Moon’s surface down to a depth of 1485.0 km in the liquid outer core in the radial direction and 0°–360° in the lateral direction. In order to solve Eqs. (1) and (2) for the 2D model, we discretize the whole Moon long r and θ directions with staggered grids as shown in Fig. 1. The two stress components are discretized into half grid spacing between the anti-plane velocities in both directions. In radial direction, the grids are distributed uniformly with constant grid spacing Δr. In lateral direction, the number of grids is the same for all radii and hence the grid spacing reduces gradually with increasing depth.

We apply a PSM/FDM hybrid scheme to solve Eqs. (1) and (2) for the discretized model on the staggered grids (Wang et al. 2014). In the lateral direction, derivatives with respect to θ are calculated with the PSM method. Let $f(j\Delta\theta)(j = 0, 1, 2, ..., N-1)$ represent discretized field variable in θ direction, we need to evaluate its derivative half-way between grids at $(j \pm 1/2)\Delta\theta$. The derivative is calculated in the wavenumber domain by multiplication and inverse Fourier transformation:

$$ \frac{d}{d\theta} f\left[\left(j \pm \frac{1}{2}\right) \Delta\theta\right] = \frac{\Delta k}{2\pi} \sum_{l=0}^{N-1} i(l\Delta k) e^{\pm i l \Delta k \Delta\theta/2} F(l\Delta k) e^{i2\pi j l/N}, \\ j = 0, 1, 2, ..., N-1, \qquad (3) $$

where $F(l\Delta k)(l = 0, 1, 2, ..., N-1)$ is the Fourier transformation of $f(j\Delta\theta)$. In the radial direction, derivatives with respect to r are calculated with a fourth-

<figure><img src="image_3.png" /><figcaption>Fig. 1 Discretization of stress and velocity components for the 2D global model on staggered grids</figcaption></figure>

order accuracy FDM on staggered grids. For field variable $f(i\Delta r)(i = 0, 1, 2, ..., M-1)$, its derivative half-way between grids at $(i \pm 1/2)\Delta r$ is calculated with:

$$ \frac{d}{dr} f(i\Delta r) = \frac{1}{\Delta r} \left\{ c_1 \left[ f\left(\left(i + \frac{1}{2}\right)\Delta r\right) - f\left(\left(i - \frac{1}{2}\right)\Delta r\right) \right] + c_2 \left[ f\left(\left(i - \frac{3}{2}\right)\Delta r\right) - f\left(\left(i + \frac{3}{2}\right)\Delta r\right) \right] \right\}, \qquad (4) $$

where $c_1 = 9/8$, $c_2 = 1/24$. The time domain calculation of $v_z$ is performed with a second-order FDM scheme (Furumura et al. 1998). Time interval $\Delta t$ is determined from the minimum grid spacing and the maximum seismic velocity in the model by satisfying stability condition for both PSM and FDM schemes (Wang et al. 2014).

The artificial boundary condition at the bottom of the model located in the outer core is treated by introducing a damping zone of 20 gird points (Cerjan et al. 1985). The free surface boundary condition at the Moon’s surface is introduced into the calculation by satisfying a traction-free condition. Body force component $f_z$ in Eq. (1) is calculated by a combination of corresponding moment tensor components as given by Wang et al. (2014). Distribution of body force component in space and source-time function is approximated with a bell-shaped Herrmann’s function of unit area (Herrmann 1979).

Lunar seismic observation and numerical modeling for P-SV wavefield (Wang et al. 2013) showed that the lunar seismograms are dominated with high-frequency waves. Hence, it is important to model high-frequency wavefield in order to enhance our understanding of lunar seismic wave propagation. Comparing with the common hybrid method used by Wang et al. (2013), the parallel PSM/FDM algorithm used in this study enables us to improve calculating efficiency and hence calculate for higher frequency wavefield. FDM used in the radial direction is a localized derivative operator and hence allows us to divide the whole Moon model into sub-domains along the radial direction. In each sub-domain, the number of grids in both radial and lateral directions is the same. Each sub-domain is assigned to a processor of a PC cluster. In each time step, spatial derivatives are calculated in each processor simultaneously. Since the fourth-order staggered grid FDM is used in the radial direction, the values of the top and bottom two layers in each sub-domain are exchanged between neighboring processors. Overlap regions of two layers of grids are added to both the top and bottom of each sub-domain to store the exchanged data. The data exchange and interprocess communication are implemented by the Message Passing Interface. The good accuracy and efficiency enhancement of the parallel PSM/FDM scheme were discussed for 3D regional model by Qin et al. (2012) and SH-
wave propagation modeling for 2D whole-Earth model by Wang et al. (2014). The following is the procedure to solve Eqs. (1) and (2) with the parallel PSM/FDM scheme at each time step: Firstly, the whole Moon model is divided into sub-domains along the radial direction and assigned to processors; In each sub-domain, the spatial derivatives of velocity in Eq. (2) are then evaluated with PSM or FDM and data in the overlap regions are exchanged; Stress components are calculated in Eq. (2) with time differencing and then their spatial derivatives in Eq. (1) are evaluated with PSM or FDM and data in the overlap regions are exchanged; Finally the velocity component in Eq. (1) are calculated with time domain differencing.

## 3 SH-wave propagation in the whole moon

We first show characteristics of global propagation of SH-wave in the whole Moon model. Recent reanalysis of Apollo seismic waveforms suggested the existence of a liquid outer core and a solid inner core. Our 1D whole Moon model is derived from the recently proposed velocity and density model by Weber et al. (2011) and the $Q_s$ value given in VPREMOON model by Garcia et al. (2011) as shown in Fig. 2. The Moon’s crust in this model is 39 km thick that includes the upper and lower layers. The upper crust is thinner but has a very low S-wave velocity (1.8 km/s) as compared to that in the lower crust (3.2 km/s). The Moon’s mantle extends to a depth of 1257.1 km with a discontinuity at 738.0 km deep. Variation of S-wave velocity in the mantle is smoother than that for the P-wave

velocity. A transition zone of 150.0 km thick was proposed by Weber et al. (2011) between the mantle and outer core. It has S-wave velocity as low as that in the lower crust. The 1-km layer just below the free surface is considered as a fractured and fragmented zone that has very low S-wave velocity.

The whole Moon is discretized into 8192 grids in the lateral direction that leads to a maximum lateral grid spacing of 1.33 km at the surface and a minimum grid spacing of 0.19 km at the base of the model. In the radial direction, the model is discretized into 1650 grids with uniform grid spacing of 0.9 km. Width of the bell-shaped source-time function is 3.0 s. The minimum wavelength in the calculation is 5.4 km considering the minimum S-wave velocity in the upper crust. The number of grid points per minimum wavelength equals to or great than 4.0 in the lateral direction and equals to 6.0 in the radial direction that is sufficient for PSM and fourth-order staggered grid FDM, respectively. The time interval determined from the minimum grid spacing, the maximum wave velocity, and the stability condition is 0.015 s in the calculation. We calculated 200,000 time steps for wave propagation of 3000 s in the whole Moon. The whole model is divided into 30 sub-domains along the radial direction and assigned to 30 processors on a PC cluster. Comparing with the modeling for global P-SV wave propagation in the whole Moon (Wang et al. 2013), the parallel modeling for SH-wave in this paper requires less parameters and less computer memory and computation time for the same number of grid points. Therefore, it is possible to discretize the whole Moon with finer grid spacing and extend the higher

<figure><img src="image_4.png" /><figcaption>**Fig. 2** The whole Moon model used in the modeling. a S-wave velocity and density distribution along depth (Weber et al. 2011), b anelastic attenuation factor for S wave (Garcia et al. 2011)</figcaption></figure>
frequency range in the modeling. Characteristics of global SH-wave propagation can be investigated at higher frequency range than that for global P-SV wave propagation.

## 3.1 The shallow moonquake

Except for the thermal moonquakes and meteoroid impacts, most identified moonquakes occurred inside the Moon are classified into shallow moonquakes and deep moonquakes. Focal depth of the 28 detected shallow moonquakes is between 50 and 220 km in the upper mantle. Number of shallow moonquakes is much less than the deep moonquakes, but their magnitudes are larger. Although the focal mechanisms of shallow events are not well understood from the Apollo observation, it is inferred that they are comparable to intraplate earthquakes and are quite likely caused by tectonic reason (Nakamura et al. 1979). Wang et al. (2013) applied a double-couple line source corresponding to a 45° dipping fault for P-SV wave modeling. For SH-wave modeling in this study, a seismic source with moment tensor component $M_{rz} = 1.0$ is used that is commonly adopted in whole-Earth seismic wave propagation modeling (Wang et al. 2014). The focal depth of the simulated moonquake is 100 km as shown in Fig. 3.

Figure 4a shows synthetic seismograms at the Moon’s surface for the shallow moonquake together with theoretical travel time curves for major phases. Direct S wave with large amplitude can be seen up to epicentral distance of 145°. Surface reflections and multiple reflections such as sS, SS, sSS, and SSS appear at epicentral distance greater than 45°. Interactions between the direct S and these reflections cause large amplitude between 45° and 150°. Core–mantle reflection ScS is strong and clear at all distances followed by sScS. Reflection from the top of the core–mantle transition zone is also visible before ScS. The core–mantle reflections (ScS)₂ and (sScS)₂ are still clearly identified at about 1300 and 1400 s. The most important features of wave propagation at later time are the development of wave trains of surface multiples such as SS, SSS, and SSSS. Comparing with SH-wave propagation in the whole Earth (Wang et al. 2014), waveforms of all phases for the Moon show strong reverberations following the first arrivals. The duration of the reverberating wave train increases with the propagation distance that is very obvious for SS, SSS, and SSSS. Such phenomenon is observed for global P-SV wave propagation in the Moon as well (Wang et al. 2013). This can be explained as the trapping wave effect in the low-velocity upper crust. The S-wave velocity in the upper crust is much lower than that in the lower crust and mantle. Transmission of energy from lower crust into this layer is rather efficient than that propagates back into the lower crust. Multiple reflections occurred in this layer causes reverberations that are observed

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Configuration of the whole Moon model used in the modeling. Shallow and deep circles show crust–mantle and core–mantle boundaries. Dashed circles are mantle discontinuity, top boundary of core–mantle transition zone, and inner–outer core boundary. Stars are locations of the shallow and deep moonquakes used in the modeling</figcaption></figure>

on the synthetic waveforms. With increasing epicentral distance, constructive interference between the multiple reflections and later incident waves enhance the reverberating wave trains. This suggests that the trapping effect and the reverberations occurred within the low-velocity upper crust which contributes partly to the long wave trains observed for the moonquake.

## 3.2 The deep moonquake

Most of the identified moonquakes are deep events occurred repeatedly at depths of 700–1200 km within about 300 distinct foci. Magnitudes of the deep events are smaller than that for the shallow events. Their occurrence is supposed to be related to the tidal force raised on the Moon by the Earth and could be slipping along a nearly horizontal plane (Nakamura 1978; Koyama and Nakamura 1980). Zhao et al. (2012) presented a new tomographic model for the Moon. They found a correlation between the distribution of deep moonquakes and seismic wave velocity variations in the deep lunar mantle that is similar to earthquakes. They suggested that the occurrence of deep moonquakes is affected by the lunar structural heterogeneity in addition to the tidal stresses that is quite similar to the shallow and deep earthquakes occurred in the crust and upper mantle of the Earth. In global P-SV wave propagation modeling for a deep moonquake, Wang et al. (2013) applied a double-couple source corresponding to a horizontal plane. For SH-wave propagation in this paper, we adopt a seismic source with moment tensor component
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Synthetic seismograms at the surface for the whole Moon model. **a** is for a 100-km and **b** is for a 900-km deep moonquake. The ray-theoretical arrival times for major phases are shown in red curves. The amplitude of all seismograms are multiplied by a factor of 1.0 × 10<sup>6</sup></figcaption></figure>

$M_{0z}$ = 1.0 that describes a horizontal slipping. The focal depth of the simulated deep moonquake is 900 km (Fig. 3).

Figure 4b shows synthetic seismograms for the deep moonquake with theoretical arrival times for major phases. We see very clear direct S waves with strong energy. Comparing with shallow moonquake, surface reflections and multiples do not interfere with direct S wave that allows S waveforms to be identified clearly at all distances. The core–mantle reflection ScS is rather obvious for deep event as compared to shallow one. Reflection from the top of core–mantle transition zone can be clearly discerned before ScS. Later reflections from core–mantle boundaries such as sScS and (ScS)<sub>2</sub> are also clearly identified at later time. Comparing with the waveforms for shallow events, reverberations following direct S and core–mantle reflections are weaker. The most significant differences between two events are the very weakly developed surface reflections and multiples for the deep one. This can be explained by the very small incidence angle for all epicentral distance for deep event for which the waves are efficiently transmitted into and out of the low-velocity upper crust layer. Hence the wave trapping effect and multiple reflections within the low-velocity layer is much less developed for deep event than that for shallow one. This is observed in case of P-SV wave propagation as well. The most important feature for the deep event is the very clearly shown reflections from the core–mantle boundary. Comparing with the shallow event (Fig. 4a) and the P-SV wave cases (Wang et al. 2013), the ScS for SH-wave is the most clear

and identifiable for a large range of epicentral distance. Garcia et al. (2011) showed that the amplitude of the core–mantle reflection for vertically polarized S wave (ScSV) is smaller than that for horizontally polarized S wave (ScSH). They suggested that ScSH can be detected from the Apollo data and can be used to support the existence of the liquid core and constrain the radius of the core rather than ScSV. Our modeling for the deep moonquake shows that the ScSH is the most clearly identified core–mantle reflection and thus supports their suggestions. Hence ScSH observed from deep moonquake should be the most efficient seismic phase for the study of the Moon’s core.

## 4 Frequency effects on synthetic waveforms

The observed moonquake waveforms are dominated with high-frequency reverberation and coda. For global P-SV wave propagation in the Moon, numerical modeling showed that amplitude and duration of reverberations increase with dominated frequency (Wang et al. 2013). In this section, we study frequency effects on synthetic waveforms for SH-wave propagation. In order to model higher frequency wavefield, we reduced the range of the model used in previous modeling from 360° to 180° in the lateral direction. The number of grids in the lateral direction is still 8192 and lateral grid spacing ranges from 0.67 km at the surface to 0.09 at the bottom of the model. The number of radial grids is 3000 with uniform grid
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Comparing of synthetic seismograms for different frequency for a shallow moonquake with focal depth of 100 km. The duration of source-time function for **a** and **b** is 3.0 and 1.67 s, respectively. The amplitude of all seismograms are multiplied by a factor of 5.0 × 10⁵</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6** Comparing of synthetic seismograms for different frequency for a deep moonquake with focal depth of 900 km. The duration of source-time function for **a** and **b** is 3.0 and 1.67 s, respectively. The amplitude of all seismograms are multiplied by a factor of 5.0 × 10⁵</figcaption></figure>

spacing of 0.5 km. With the reduced grid spacing, the width of the bell-shaped source-time function reduces to 1.67 s. The minimum wavelength in the calculation is 3.0 km. Time interval determined from stability condition is 0.01 s. The calculation is carried out for 150,000 time steps to show wave propagation of a duration of 1500 s. The whole model is divided into 30 sub-domains along radial direction and assigned to 30 processors. In order to compare with waveforms of the previous models, the same focal mechanisms are used in the modeling.

Figure 5 shows comparing of synthetic waveforms up to epicentral distance of 150°. Figure 5a shows results calculated for the 360° model with source-time function of 3.0 s duration in previous section. Figure 5b is the result for the 180° model with source-time function of 1.67 s duration calculated in this section. We see that the amplitude of the direct wave and the following reverberations increases with increasing frequency. For higher frequency model, the reverberations decay slower and their duration becomes longer. Hence the reverberations enhance with increasing frequency. Similar phenomenon is observed for P-SV wave propagation in the Moon as well (Wang et al. 2013). The reason is that for shorter wavelength component, the trapping effect and multiple reflections occurred within the low-velocity layer is stronger. The comparison indicates that the seismic phases could be identified more clearly on lower frequency waveforms.

Comparing of synthetic waveforms of different frequency for deep moonquake is shown in Fig. 6. The strength and duration of reverberations following each phase increase with frequency. However, identification of

direct wave and reflections from core–mantle boundary become less distinct with increasing frequency. Comparing with shallow moonquake, the enhancement of reverberation with frequency is weaker. This is because waves from deep moonquake propagate into the low-velocity layer with very small incidence angle at all epicentral distances. The trapping effect and multiple reflections occurred within the low-velocity layer is weaker for deep event than that for shallow one. Garcia et al. (2011) compared observed waveforms at two different frequency ranges and found that the energy of ScSH is lower at high frequencies than at low frequencies and suggested that the search for this seismic phase should be performed at frequencies as low as possible. Results of our comparison show similar phenomenon and support their suggestion.

## 5 Effects of crustal thickness on synthetic waveforms

For seismic wave propagation on the Earth, numerical modeling showed that the low-velocity near-surface layer has significant effects on seismic body wave and generates seismic coda (e.g., Robertsson and Holliger 1997). Analysis of observations and numerical modeling suggested that the lateral variation of crust thickness and velocity has pronounced effects on regional seismic wave propagation and the thick zone of lower seismic velocities in the continental crust acts to trap S-wave energy and generates crustally guided phases (Furumura and Kennett 1997; Kennett and Furumura 2001). The most recent research shows that the lateral variation of crust thickness on the
Moon is much larger than that on the Earth, especially across some young impact basins where crust varies from near zero to about 60 km drastically (Hikida and Wieczorek 2007; Wieczorek et al. 2013). Hence, it is very interesting to investigate the effects of lateral variation of crust thickness on wave propagation on the Moon.

Similar to the modeling performed in Sect. 3, we calculated SH-wave propagation for a 100-km shallow moonquake and a 900-km deep moonquake. The whole Moon model, the number of grids, grid spacing, time interval, width of the source-time function, and focal mechanism are the same as used in Sect. 3 except for that the crust thickness varies laterally. Considering the average crust thickness of 40 km on the Moon, the crust incorporated in the modeling varies from 20 to 60 km laterally. The top 8 km and 20 km layers are set to be the upper crust for the 20 and 60 km crust, respectively. For the shallow moonquake modeling, the crust is 20 km thick between 0° and 180° and 60 km thick from 180° to 360°. For the deep moonquake modeling, thickness of the crust is 60 km between 0° and 180° and 20 km between 180° and 360°. The seismic source is located at 180° for both models. With such models, it is possible to show difference of wave propagation caused by varying crust thickness.

Synthetic seismograms for both the shallow and deep moonquakes are shown in Fig. 7. For shallow moonquake, the basic features of synthetic waveforms for both the thin and thick crust parts are similar to that for the whole model shown in Fig. 3. We see clear direct wave and core–mantle reflections. Surface reflections and multiples dominate the latter of the wavefield with increasing strength. Comparing between the thin and thick crust regions shows that the reverberations following each phases become stronger in the thick crust region. This suggests that wave trapping effects and multiple reflections are more efficient in the thick crust region. Such phenomenon is similar to the S-wave trapping effects in thick zone of lower velocity continental crust on the Earth (Kennett and Furumura 2001). One interesting phenomenon that is not observed in Fig. 3 is the very efficient transmission of waves from mantle into the low-velocity thick crust region occurred at the boundaries between two regions at 180° and 360°. The transmitted waves are very strong and propagate as trapped waves and dominate the whole wavefield above the 60-km thick crust region. Hence the identification of each phase is not so easy at this region as compared to the 20-km crust region. For deep moonquake, the basic features of wave propagation are similar to the model shown in Fig. 4 as well. We still observe the enhancement of reverberations following each phase at the thick crust region. Transmission of waves at the boundaries at 180° and 360° from mantle into the thick crust region are efficient and cause multiple reflections that contaminate the wavefield.

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Synthetic seismograms at the surface for the whole Moon model with lateral varying crust thickness. The *upper* and *lower* panels show results for 900 km deep and 100 km shallow moonquakes with crust thickness distribution, respectively. “H” means the thickness of the crust thickness. The seismic source is located at 180°. The amplitude of all seismograms are multiplied by a factor of $1.0 \times 10^6$.</figcaption></figure>

However, the effects of crust thickness on wave propagation are not so obvious for deep moonquake than that for shallow one and suggest that each phase could be more easily identified on deep moonquake waveforms even considering the lateral variation of crust thickness.

## 6 Comparison with the observed waveforms

In this section, we try to make comparison between our synthetics with Apollo observation data tentatively to show features of wave propagation that exhibited in our results and those not reproduced by our modeling. We compared with two shallow (SH1, SH2) and one deep (A1) events typically recorded on three Apollo stations (S12, S14, S16).
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 8** Comparison between observed transverse component Apollo seismograms (black) with synthetic SH-waveforms (red). Blue lines show theoretical arrival times for major annotated phases. The maximum amplitude of observation (unit in DU) and synthetics (unit in cm/s) are annotated for each waveform. The event-station pair for each figure is: a SH1-S12, b SH1-S16, c SH2-S12, d SH2-S14, e A1-S12, f A1-S16. Epicentral distance for each pair: SH1-S12: 89.13°, SH1-S16: 85.73°, SH2-S12: 69.73°, SH2-S14: 67.43°, A1-S12: 18.15°, A1-S16: 51.22°</figcaption></figure>

These events were used for comparisons with global P-SV wave propagation modeling as well (Wang et al. 2013). The original waveform data are processed with the same procedure as used by Wang et al. (2013). We first band-pass filter the data between 0.1 Hz and the Nyquist frequency (3.3127 Hz), remove large amplitude noise by a despiking algorithm, and then correct the gain difference between the two horizontal components. The two corrected horizontal components are finally rotated into the radial and transverse components for comparison. Wang et al. (2013) compared the radial and vertical components for P-SV waves. In this paper, we compare the obtained transverse component with our synthetic SH-waveform. Since our calculation is performed for a 2D model, we firstly corrected our synthetic waveforms from 2D to 3D for geometric spreading and source radiation and then band-pass filtered them with the same frequency band as used for observations before comparing.

Figure 8 shows comparison of waveforms at each event-station pair. For shallow events, the calculated direct

S wave corresponds to a sharp increase of amplitude in observation that was identified as S wave. Because of the trapping effects and multiple reflections within the low-velocity surface layer, the direct S wave, the following sS and core-mantle reflections form reverberating wave trains lasting for about 400 s and they are not clearly separated. The wave train following the direct S wave is consistent with the observations and suggests that it contributes mainly to this portion of the observations. The second wave train lasting for about 1000 s is formed by the reverberating surface multiples SS, SSS, and SSSS. They contribute mainly to the middle portion of the observations and are consistent with data. At later time of the synthetics, reverberations caused by trapped waves in the surface layer agree well with observations and suggest their contribution at later portion. For deep event, the most interesting feature is the clearly identified and separated core-mantle reflections on synthetics. However, the corresponding arrivals cannot be identified on observations that are similar to the case for shallow events. Each phase shows reverberating
waveform caused by surface low-velocity layer that is
consistent with observation. With increasing epicentral
distance, strength and duration of the reverberating wave
train and their contribution to the waveform increase.

Although our modeling produced the reverberating
features of each wave train that are consistent with the
observations, the decay of our synthetic wave train is not as
slow as that shown in the data. The strong scattering effects
between seismic phases presented on observation are not
reproduced in our synthetics. A similar phenomenon is also
shown for global P-SV wave propagation modeling (Wang
et al. 2013). This suggests that the low-velocity surface
layer presented in our model contributes significantly to the
development of reverberating waveforms, but it is not
enough to produce the strong scattering and slow decay
features of the observations with the present layered model.
Similar phenomenon can be observed for earthquake as
well. The spindle-shaped waveforms frequently appeared
in the crustal earthquakes due to wave trapping effect in
low-velocity layer (e.g., Furumura and Kennett 1997;
Kennett and Furumura 2001). These features also appeared
due to strong seismic wave scattering from small-scale
heterogeneity with characteristic scale of several kilome-
ters (e.g., Nielsen et al. 2003; Furumura and Kennett 2005;
Takahashi et al. 2007; Kennett and Furumura 2008; Sato
et al. 2012; Takemura and Yoshimoto 2014). This implies
that small-scale heterogeneities should be required for
modeling of realistic Moon structure. For earthquake ob-
servation, there is no wave before the direct SH-wave on
transverse component for isotropic medium. However,
transverse components for moonquakes in Fig. 7 show very
strong energy before direct SH-wave and after direct P
wave that do not appear on our synthetics. This discrepancy
again suggests that the present isotropic 1D layered model
is not enough to produce the characteristics of observations
and more realistic models should be considered in order to
well understand seismic wave propagation in the Moon.

7 Conclusions

We applied a hybrid PSM/FDM method to model global
SH-wave propagation in the whole Moon and try to im-
prove our understanding of lunar seismic wave propaga-
tion. In order to model high-frequency wavefield as
presented on lunar seismic observations, we implemented
parallel modeling on a PC cluster to improve the comput-
ing efficiency. The modeling is performed for shallow and
deep moonquakes, respectively, to show SH-wave
propagation in the whole Moon. Effects of frequency range
and lateral variation of crust thickness on wave propagation
are investigated with several models. Finally, we compared
our synthetics with Apollo observation data and showed

major features of wave propagation that are produced by
our modeling and the limitations of our present models.

For shallow moonquake, the direct wave and core–mantle
reflections can be identified at a limited range of epicentral
distance because of interactions between them and other
phases. Waveform of each phase shows strong reverberations
caused by wave trapping effects and multiple reflections oc-
curred in the surface low-velocity upper crust. Strong rever-
berating surface multiples dominate later in the wave
propagation. For deep moonquake, the direct wave and core–
mantle reflections can be clearly identified over wide range of
epicentral distance. Reverberations in waveform and surface
multiples are weaker than those for shallow event.

Increasing frequency enhances the strength and duration
of waveform reverberation and reduces the possibility of
phase identification. This suggests that body wave phases
such as core–mantle reflections can be more clearly iden-
tified for deep moonquake at low-frequency range.
Development and deployment of very broad band seis-
mometers in the future lunar mission could provide further
data to constrain the inner structure of the Moon. Lateral
variation of crust thickness significantly influences lunar
seismic wave propagation. Wavefield above the thick crust
region tends to be more complicated because of transmitted
wave from mantle into the crust.

Comparing with Apollo observation data that show that,
our modeling produces the reverberating features of each
phase but cannot reproduce the strong and slow decay
scattering between each phase. For earthquake observation,
similar phenomenon can be explained with wave trapping
effect in low-velocity layer and scattering from small-scale
heterogeneity. This suggests that more realistic Moon
model should be considered in the future numerical mod-
eling in order to improve our understanding of lunar seis-
mic wave propagation.

**Acknowledgments** This research was supported by the National Natural Science Foundation of China (Grants 41374046 and 41174034).

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

References

Blanchette-Guertin J-F, Johnson CL, Lawrence JF (2012) Modeling seismic waveforms in a highly scattering Moon 43rd lunar and planetary science conference Abstract 1473 Lunar & Planetary Institute Huston

<figure><img src="image_16.png" /></figure>
Cerjan C, Kosloff D, Kosloff R, Reshef M (1985) A nonreflecting boundary condition for discrete acoustic and elastic wave equation. Geophysics 50:705–708

Chenet H, Lognonné P, Wieczorek M, Mizutani H (2006) Lateral variations of lunar crustal thickness from the Apollo seismic data set. Earth Planet Sci Lett 243:1–14

Cummins PR, Geller RJ, Hatori T, Takeuchi N (1994) DSM complete synthetic seismograms: SH, spherically symmetric, case. Geophys Res Lett 21:533–536

Dainty AM, Toksöz MN, Anderson KR, Pines PJ, Nakamura Y, Latham G (1974) Seismic scattering and shallow structure of the moon in Oceanus Procellarum. Moon 91:11–29

Furumura T, Kennett BLN, Furumura M (1998) Seismic wavefield calculation for laterally heterogeneous whole Earth models using the pseudospectral method. Geophys J Int 135:845–860

Furumura T, Kennett BLN (1997) On the nature of regional seismic phases—II. On the influence of structural barriers. Geophys J Int 129:221–234

Furumura T, Kennett BLN (2005) Subduction zone guided waves and the heterogeneity structure of the subducted plate: intensity anomalies in northern Japan. J Geophys Res 110:B10302. doi:10.1029/2004JB003486

Garcia RF, Gagnepain-Beyneix J, Chevrot S, Lognonné P (2011) Very preliminary reference Moon model. Phys Earth Planet Inter 188:96–113. doi:10.1016/j.pepi.2011.06.015

Gagnepain-Beyneix J, Lognonné P, Chenet H, Lombardi D, Spohn T (2006) A seismic model of the lunar mantle and constraints on temperature and mineralogy. Phys Earth Planet Inter 159:140–166

Herrmann RB (1979) SH-wave generation by dislocation source-a numerical study. Bull. Seismol Soc Am 69:1–15

Hikida H, Wieczorek MA (2007) Crustal thickness of the moon: new constraints from gravity inversions using polyhedral shape models. Icarus 192:150–166

Igel H, Weber M (1995) SH-wave propagation in the whole mantle using high-order finite differences. Geophys Res Lett 22:731–734

Igel H, Gudmundsson O (1997) Frequency-dependent effects on travel times and waveforms of long-period S and SS waves. Phys Earth Planet Inter 104:229–246

Jahnke G, Thorne MS, Cochard A, Igel H (2008) Global SH-wave propagation using a parallel axisymmetric spherical finite-difference scheme: application to whole mantle scattering. Geophys J Int 173:815–826. doi:10.1111/j.1365-246X.2008.03744.x

Kennett BLN, Furumura T (2001) Regional phases in continental and oceanic environments. Geophys J Int 146:562–568. doi:10.1046/j.1365-246x.2001.01467.x

Kennett BLN, Furumura T (2008) Stochastic waveguide in the lithosphere: Indonesian subduction zone to Australian craton. Geophys J Int 172:363–382. doi:10.1111/j.1365-246X.2007.03647.x

Khan A, Mosegaard K (2002) An inquiry into the lunar interior: a nonlinear inversion of the Apollo seismic data. J Geophy Res 107. doi:10.1029/2001JE001658

Khan A, Maclennan J, Taylor SR, Connolly JAD (2006) Are the Earth and the Moon compositionally alike? Inferences on lunar composition and implications for lunar origin and evolution from geophysical modeling. J Geophys Res 111:E05005. doi:10.1029/2005JE002608

Khan A, Connolly JAD, Maclennan J, Mosegaard K (2007) Joint inversion of seismic and gravity data for lunar composition and thermal state. Geophys J Int 168:243–258. doi:10.1111/j.1365-246X.2006.03200.x

Koyama J, Nakamura Y (1980) Focal mechanism of deep moon-quakes. In: Bedini (ed) Lunar and Planetary Science Conference Proceedings, vol. 11 of Lunar and Planetary Science Conference Proceedings, pp 1855–1865

Lawrence JF, Johnson CL (2010) Synthetic seismograms with high-frequency scattering for the Moon 41st Lunar and Planetary Science Conference Abstract 2701 Lunar & Planetary Institute Huston

Lognonné P, Gagnepain-Beyneix J, Chenet H (2003) A new seismic model of the moon: implications for structure, thermal evolution and formation of the moon. Earth Planet Sci Lett 211:27–44

Lognonné P (2005) Planetary seismology. Ann Rev Earth Planet Sci 33:191–1934

Nakamura Y (1977) Seismic energy transmission in an intensively scattering environment. J Geophys Res 43:389–399

Nakamura Y (1978) Al moonquakes—source distribution and mechanism In: Lunar and Planetary Science Conference Proceedings, vol. 9 of Lunar and Planetary Science Conference Proceedings, pp 3589–3607

Nakamura Y (1983) Seismic velocity structure of the lunar mantle. J Geophys Res 88:677–686

Nakamura Y, Latham GV, Dorman HJ, Ibrahim A, Koyama J, Horvath P (1979) Shallow moonquakes—depth, distribution and implications as to the present state of the lunar interior. In: Hinners (ed) Lunar and Planetary Science Conference Proceedings, vol. 10. Pergamon Press, New York, pp 2299–2309

Nielsen L, Thybo H, Levander A, Solodilov N (2003) Origin of upper-mantle seismic scattering—evidence from Russian peaceful nuclear explosion data. Geophys J Int 154:196–204

Qin Y, Wang Y, Takenaka H, Zhang X (2012) Seismic ground motion amplification in a 3D sedimentary basin: effect of the vertical velocity gradient. J Geophys Eng 9:761–772. doi:10.1088/1742-2132/9/6/761

Robertsson JOA, Holliger K (1997) Modeling of seismic wave propagation near the Earth’s surface. Phys Earth Planet Inter 104:193–211

Sato H, Fehler M, Maeda T (2012) seismic wave propagation and scattering in the heterogeneous Earth structure, 2nd edn. Springer-Verlag, New York

Takahashi T, Sato H, Nishimura T, Obara K (2007) Strong inhomogeneity beneath quaternary volcanoes revealed from the peak delay analysis of S-wave seismograms of microearthquakes in northeastern Japan. Geophys J Int 168:90–99. doi:10.1111/j.1365-246X.2006.03197.x

Takemura S, Yoshimoto K (2014) Strong seismic wave scattering in the low-velocity anomaly associated with subduction of oceanic plate. Geophys J Int 197:1016–1032. doi:10.1093/gji/ggu031

Thorne MS, Lay T, Garnero EJ, Jahnke G, Igel H (2007) Seismic imaging of the laterally varying D'' region beneath the cocos plate. Geophys J Int 170:635–648. doi:10.1111/j.1365-246X.2006.03279.x

Toyokuni G, Takenaka H, Wang Y, Kennett BLN (2005) Quasi-spherical approach for seismic wave modeling in a 2D slice of a global Earth model with lateral heterogeneity. Geophys Res. Lett 32:L09305. doi:10.1029/2004GL022180

Wang T, Chen L (2009) Distinct velocity variations around the base of the upper mantle beneath northeast Asia. Phys Earth Planet Inter. doi:10.1016/j.pepi.2008.09.021

Wang Y, Takenaka H (2011) SH-wavefield simulation for a laterally heterogeneous whole-Earth model using the pseudospectral method. Sci China Earth Sci 54:1940–1947. doi:10.1007/s11430-011-4244-8

Wang Y, Takenaka H, Jiang X, Lei J (2013) Modelling two-dimensional global seismic wave propagation in a laterally heterogeneous whole-Moon model. Geophys J Int. doi:10.1093/gji/ggs094

Wang Y, Luo L, Qin Y, Zhang X (2014) Global SH-wavefield calculation for a two dimensional whole-Earth model with the Proceedings, vol. 11 of Lunar and Planetary Science Conference Proceedings, pp 1855–1865
parallel hybrid PSM/FDM algorithm. Earthq Sci. doi:10.1007/s11589-014-0085-9

Weber RC, Lin P, Garnero EJ, Williams Q, Lagnonné P (2011) Seismic detection of the lunar core. Science 331:309–312

Wen L (2002) An SH hybrid method and shear velocity structures in the lowermost mantle beneath the central Pacific and South Atlantic Oceans. J Geophys Res 107:2055. doi:10.1029/2001JB000499

Wieczorek MA, Neumann GA, Nimmo F, Kiefer WS, Taylor GJ, Melosh HJ, Phillips RJ, Solomon SC, Andrews-Hanna JC,

Asmar SW et al (2013) The crust of the Moon as seen by GRAIL. Science 339:671–675

Zhao D, Lei J, Liu L (2008) Seismic tomography of the Moon. Chin Sci Bull 53:3897–3970

Zhao D, Arai T, Liu L, Ohtani E (2012) Seismic tomography and geochemical evidence for lunar mantle heterogeneity: comparing with Earth. Global Planet Change 90:29–36