1 Introduction

“We leave as we came and, God willing, as we shall return with peace and hope for all mankind.” These are the words of the Apollo 17 astronaut Gene Cernan, the last human to walk on the Moon, moments before he climbed aboard the lunar module for the return trip to Earth in 1972. With six manned landings from 1969 to 1972 and numerous unmanned visits, Moon is one of the most visited celestial bodies. However, for about the past 5 decades, no one has walked on the Moon. This may soon change as the National Aeronautics and Space Administration (NASA) recently unveiled a plan to not only land humankind on the Moon over the next few years, but also to construct a permanent lunar base (NASA 2019). Also, China successfully landed an unmanned space probe on the “far side” of the Moon making the Chinese spacecraft the first to ever land on this unexplored area of the Moon.

All these developments have placed lunar explorations at the center of attention once again. As part of the lunar exploration program, active and passive seismic data were acquired during various Apollo missions in the past. The objective of the Active Seismic Experiment (ASE) (Apollo 14 and 16) and Lunar Seismic Profiling Experiment (LSPE) (Apollo 17) was to produce and collect seismic waves to study the internal structure of the Moon. At the time, seismic refraction was the prominent seismic geophysical method for subsurface explorations. Therefore, it is not surprising that even up to now the majority of the lunar studies have analyzed lunar seismic data using the first arrivals of primary (P-) waves (e.g., Kovach and Watkins 1973a, b; Watkins and Kovach 1973; Cooper and Kovach 1974; Heffels et al. 2017). Generally, these exploration efforts have highlighted that the lunar regolith is composed of a very soft powder-like material interspersed with breccias and rocks of different sizes as well as massive boulders (e.g., Dainty et al. 1974; Dal Moro 2015b). The soft regolith material is a layer of unconsolidated debris of fine soil characterized by Vs values in the range of 30–60 m/s and VP of 104–114 m/s (AA. VV. 1972; Dal Moro 2015a, b). The properties of the soft regolith have been found to be consistent across different sites regardless of the location (Watkins and Kovach 1973).

As stated earlier, the general lunar regolith structure have been estimated using P-waves, and the contribution from shear (S-) waves and their interactions with P-waves (i.e., surface waves such as Rayleigh waves) have not been employed as much as they have deserved given their advantages (Sollberger et al. 2016). The propagation of surface waves contains useful information about the near-surface structure of a given medium. The dispersive nature of surface waves in a vertically-heterogeneous domain is the basis for surface wave testing methods. In such domains, different frequency components of a waveform travel at different phase velocities depending on the medium properties (e.g., stratification, density, and stiffness of each layer). This occurs because the different frequency components correspond to different wavelengths that sample various depths (Park et al. 1999). Longer wavelengths penetrate deeper into the subsurface and short wavelengths sample the shallower depths. The dispersion pattern of surface waves at a particular site therefore contains useful information about the properties of its subsurface. Establishing a phase velocity–frequency relationship (i.e., dispersion curve) allows an estimation of the underlying velocity profile. This occurs through solution of an iterative inversion problem where measured dispersion curves are matched to theoretical dispersion curves derived from wave propagation forward modeling through an idealized subsurface domain. The use of surface waves offers a number of advantages over seismic refraction including, but not limited to, detecting velocity reversals and providing higher signal-to-noise ratios (Foti et al. 2015). Moreover, surface wave methods measure the shear wave velocity (Vs) that is controlled by the soil skeleton and has a stronger correlation with geotechnical parameters (Foti et al. 2015).

Despite the earliest applications of surface waves in the 1950s (Van der Pol 1951; Jones 1955, 1958), surface waves methods were not widely used until the 1980s when the Spectral Analysis of Surface Waves (SASW) method was developed to process the surface wave information of seismic records (Nazarian et al. 1983; Stokoe et al. 1994). Eventually, the Multichannel Analysis of Surface Waves (MASW) method was developed in the late 1990s by a team at the Kansas Geological Survey (KGS) to address some of the limitations of the SASW and to provide practitioners with a robust surface wave testing method (Park et al. 1999; Xia et al. 1999). MASW uses several receivers along the ground surface to measure surface waves from either active sources or background seismic noise. The dispersion information derived from the multichannel recording is used to deduce a Vs profile, which is a direct measurement of the small-strain shear modulus and a proxy for subsurface stiffness at a given site. Thus far, MASW has been used in different engineering applications including seismic site classification (e.g., Kanli et al. 2006; Anbazhagan and Sitharam 2010; Coe et al. 2016), evaluation of liquefaction potential (e.g., Lin et al. 2004), bedrock profiling (e.g., Miller et al. 1999; Casto et al. 2010), assessment of ground improvement (e.g., Burke and Schofield 2008; Waddell et al. 2010; Finno et al. 2015; Mahvelati et al. 2016), and detection of underground anomalies (e.g., Miller et al. 2000; Ivanov et al. 2003).

Very few studies have employed active and/or passive surface wave measurements to produce lunar subsurface structures. Early applications of surface waves employed the horizontal-to-vertical spectral ratio (HVSR) method to estimate regolith thickness (e.g., Mark and Sutton 1975; Nakamura et al. 1975). It was only a few decades later that major technical advancement in surface wave analysis caused a new wave of studies in lunar near-surface characterization. Larose et al. (2005) used Apollo 17 passive measurements to generate Rayleigh wave dispersion images. Tanimoto et al. (2008) developed near-surface (< 15 m) velocity models for the Apollo 17 landing site by applying noise cross-correlation. Yeluru et al. (2008) proposed a method to perform MASW surveys with random receiver arrays for future space explorations. Sens-Schönfelder and Larose (2010) obtained a near-surface velocity structure for the Apollo 17 landing site using a single-trace Rayleigh-wave group velocity analysis. Dal Moro (2013, 2015a, b) performed joint analysis of single-trace Rayleigh-wave group velocity and HVSR on Apollo 14 and 16 data to retrieve Vs information for the uppermost 200 m of the Moon. In a recent study, Tsuji et al. (2018) performed seismic refraction, reflection, and surface wave testing on a laboratory model filled with lunar regolith simulant. Despite these efforts, a multichannel approach that describes the propagation of surface waves in a phase velocity–frequency domain has not been applied to any of the lunar ASE data. A lunar MASW study could provide additional insight into the fundamental nature of the lunar regolith and potentially address heterogeneity issues related to the previous body wave efforts and single-station surface wave approaches.

One of the challenges that has likely prevented such analysis is the lack of true multichannel data from previous exploration efforts. A Multichannel Simulation with One Receiver (MSOR) approach could address this limitation and generate an MASW record, which would allow a multichannel analysis to be performed. MSOR generates a multichannel record by only using a single receiver attached to a fixed location and multiple impacts delivered at successively increasing distance from the receiver (Ryden et al. 2001). Combining all the individual traces would then produce a multichannel record. MSOR assumes that there is source-receiver reciprocity whereby the recorded waveforms would be the same if source and receiver were interchanged (Rayleigh 1873; Knopoff and Gangi 1959).

Previous numerical simulations have confirmed that MASW and MSOR do indeed yield similar dispersion curves in one-dimensional layered profiles, soil profiles with a dipping interface, or even in profiles with a vertical fault (Lin and Ashlock 2016). Moreover, MSOR has been successfully applied in several real-world field experiments as a viable alternative to MASW to characterize soil sites (e.g., Tokeshi et al. 2013; Lin and Ashlock 2016; Dal Moro et al. 2018; Kumar and Mahajan 2020) and pavements (e.g., Ryden et al. 2001, 2004; Yuan et al. 2014; Lin and Ashlock 2015). Direct comparisons made between MASW and MSOR at soil sites that exhibit natural spatial variability comparable to the lunar subsurface demonstrated that MASW and MSOR are equivalent for practical purposes (Lin and Ashlock 2016). Similar comparisons demonstrated that for testing on pavements, MASW provides surface wave energy up to higher frequencies than MSOR (Lin and Ashlock 2015). However, in the frequency range covered by both, MSOR and MASW exhibited similar phase velocities on average. Although MSOR alleviates the problem of having only a few receivers, it necessitates the use of a repeatable impact source that generates waves with consistent timing (Park et al. 2002). Moreover, it is still somewhat unclear whether MSOR maintains reciprocity with a multichannel receiver approach given the nature of the significant variability present in the lunar near surface (i.e., interspersed soft lunar regolith and boulders).

Given the preceding discussion, the objective of this study was to re-analyze the active lunar seismic records in the context of a multichannel surface wave approach. In doing so, first, we used numerical simulation of seismic wave propagation to evaluate the differences between the results of a true MASW record and a multichannel record compiled with the MSOR approach in two highly heterogeneous domains with different scales of fluctuations. Then, we produced multichannel records from lunar active data acquired during Apollo 16 and processed such records for their surface wave information. For the purpose of this study, the focus was on the active seismic records from the Apollo 16 mission.

2 Methods

2.1 Apollo 16 Active Seismic Experiment and Waveforms

The signals used in this study were obtained during the Apollo 16 Active Seismic Experiment. Figure 1a shows the approximate location of landing sites of different Apollo missions, all of which were located on the near-side of the Moon. The Descartes Highlands is an area of lunar highlands located on the near-side that served as the landing site of Apollo 16 (Fig. 2a). It is located at a latitude of 8° 59′ 29″ S and a longitude of 15° 30′ 52″ E (AA.VV. 1972). The active survey line of Apollo 16 consisted of 3 geophones deployed at 45.72 m intervals, and two forms of sources were used to produce active seismic energy (Fig. 1b). The first source, called the “thumper”, was a hand-held source that contained 19 small explosives (Fig. 1c). Firing the explosives would drive the base of the thumper to collide with the regolith surface and produce input energy. The thumper source was used at locations within the array spread. Examining the frequency content of thumper firings has shown that the source has a predominant frequency of 22 Hz (AA. VV. 1972). The second source was a mortar launcher that launched four grenades to design impact locations at 150, 300, 900, and 1500 m away from the receiver array (AA.VV. 1972). The focus of this study was to evaluate the lunar near-surface structure, so we only used the waveforms from the thumper source. This source was less energetic and was fired at much closer offsets from the receivers. Therefore, it would generate surface waves with more high frequency energy. Figure 2b provides a schematic representation of the linear receiver array and thumper shot locations.

Fig. 1
figure 1

a Landing sites of Apollo mission (NASA/GSFC/Arizona State University); b Apollo 16 geophone/thumper line (NASA.gov; AS16-113-18345); c John Young practices with the thumper on Earth (NASA.gov; 72-H-226)

Full size image
Fig. 2
figure 2

a Aerial image of the Apollo 16 landing site (NASA/GSFC/Arizona State University); b schematic of the receiver layout and thumper source locations

Full size image

Figure 3 plots a sample trace of a thumper firing performed during the Apollo 16 exploration. There are two distinctive features about many of the lunar active records. First, many suffer from clipped (saturated) signals, especially at close source offset locations due to the limited dynamic range of the equipment used at the time. Second, the seismic traces exhibit a significant amount of ringy behavior. The ringy and scattered nature of traces has been attributed to low attenuation and the highly heterogeneous nature of the lunar regolith (Dainty et al. 1974; Dal Moro 2015b). Multiple studies have suggested that the low attenuation of the regolith is most likely due to multiple factors, including the absence of fluids in a highly-fractured material, the vacuum conditions present in the atmosphere, and the high temperature oscillations (Tittmann 1972, 1977; Dainty et al. 1974). These ringy signals were used to generate MSOR dispersion images and mimic multichannel analysis for the acquired surface waves.

Fig. 3
figure 3

An example of an active trace from the thumper source

Full size image

2.2 Wave Propagation Simulation in a Highly-Heterogeneous Domain

Prior to MSOR analysis, we conducted a numerical study that explored whether the MSOR assumption of source-receiver reciprocity was valid to model MASW in a highly heterogeneous domain similar to the one present in the near surface of the Moon. Previous studies have identified this reciprocity when characterizing terrestrial soil sites with spatial variability (Lin and Ashlock 2016), but the distinct nature of the lunar subsurface with its interspersed soft regolith and significant presence of boulders warranted additional study. We therefore employed the Spectral Element Method (SEM) research code SPECFEM2D (Tromp et al. 2008) to simulate wave propagation throughout a lunar regolith model.

The domain of interest was a single layer 10.0 m in depth and 55.0 m in length and was spatially discretized with a mesh size of 0.2 m × 0.2 m. As shown in Fig. 4, the bi-material medium modeled a soft powder stratum that contains rock-type boulders of different sizes. Vs values of the softer material ranged from 30 to 70 m/s centered around a mean of 50 m/s (Dal Moro 2015a). Within this range, 65% of elements have a Vs of 40–60 m/s meaning that the elements are not distributed uniformly. The stiff material corresponding to the boulders interspersed within the soft lunar regolith was assigned a velocity of 1000 m/s. To produce this model, an algorithm in MATLAB® was initially used to generate a spatially-correlated random field of Vs. The MATLAB® script applies a Gaussian correlation function to generate the random Vs field through lower-upper (LU) decomposition of the covariance matrix (Constantine and Wang 2012). The initial scale of fluctuation along the vertical and horizontal directions was set to θz = 1.0 m and θx = 4.0 m, which is consistent with terrestrial geologic materials (Phoon and Kulhawy 1999). Then, Vs values greater than 70 m/s were increased to 1000 m/s to represent a rock-type material. Following this model preparation step, most of the anomalous high-velocity bodies ended up having comparable dimensions in both directions (~ 1–2 m), with some more angular than the others. These dimensions agree well with recent efforts that indicate boulder sizes of 1–4 m are quite common across the Apollo 16 landing site (Watkins et al. 2019). Both the material categories were elastic with large quality factors (Q) to mimic the low amounts of attenuation in the lunar regolith. As evident in Fig. 4, the domain was subsequently dominated by the soft material.

Fig. 4
figure 4

Numerical setups in highly heterogeneous bi-material model: a MASW; and b MSOR

Full size image

During forward modeling, we specified a stress-free boundary condition to the top surface to mimic the ground surface and used Stacey absorbing boundaries along the vertical and bottom horizontal boundaries to prevent undesired wave reflections (Clayton and Engquist 1977; Stacey 1988). All the input sources were 20 Hz Ricker wavelets. A marching time step (dt) of 8.0 × 10− 6 s guaranteed the safety of forward modeling against Courant instability (Courant et al. 1928). We then ran each forward model for 625,000 time-steps, which resulted in a total duration of 5.0 s.

We used the synthetic model to compare the dispersion images collected from a true multichannel record and those collected using the MSOR approach. For the MASW record, the 20 Hz Ricker wavelet was placed at x = 2.86 m and the subsequent waveforms were collected by a series of 10 receivers equally positioned from x = 12.00 to 53.13 m (Fig. 4a). We compiled the 10-channel MSOR record using a single receiver at x = 53.13 m while placing multiple sources at 10 equidistant offset locations between x = 44.00 to x = 2.86 m (Fig. 4b).

3 Results and Discussion

3.1 Comparison of MSOR and MASW in a Highly Heterogeneous Domain

Figure 5 plots a sample of the waveform acquired by a receiver when simulating MASW with the source and receiver 22.85 m apart. The seismic trace appears to exhibit a ringy quality similar to the experimental seismogram acquired during the Apollo 16 ASE (Fig. 3). As previously noted, this attribute has been partially attributed to scattering from the heterogeneities present in the lunar near surface (Dainty et al. 1974; Dal Moro 2015a, b). The numerical modeling reinforces this theory that the complex nature of the lunar seismic traces is in part because of the heterogeneity in the shallow subsurface and scattering of the waves.

Fig. 5
figure 5

Sample seismogram from the wave propagation modeling in the highly-heterogeneous model

Full size image

Given the complexities introduced to the seismic traces by wave scattering, we evaluated the MSOR assumption of source-receiver reciprocity for the particular structure present in the lunar heterogeneities. This was accomplished by comparing dispersion images acquired from MSOR and MASW on the same domain but with different scatterer sizes. The scattering of seismic wave propagation in a medium is controlled by the relative size of heterogeneities (θ) compared to the dominant wavelength (λ) of the seismic wave. Different scattering regimes are expected based on the ratio of θ/λ. For instance, when the heterogeneous scale lengths are small compared to the seismic wavelength (θ/λ < 0.002) scattering effects are negligibly small (Yoon 2005). However, that is not the case in the lunar regolith model. For a dominant background Vs of 50 m/s, the dominant wavelength from a 20 Hz Ricker signal is 2.5 m which is on the same order as the size of the heterogeneous features. More specifically, θ/λ ratios between 0.02 and 1.59 (0.04 m < θ < 4 m for a λ of 2.5 m) fall into the Mie scattering regime (i.e., resonant scattering) where the scale of anomalies is roughly in the order of the seismic wavelength. In such cases, the incident waves are scattered with large angles relative to the incident direction and a reflection coda is expected in the wavefield.

Figure 6 presents the numerical dispersion images of MASW and MSOR records from the model in Fig. 4. The dispersion images were generated using the SeisImager/SW software package that uses the phase shift method to transfer the waveform data from the space-time domain into the phase velocity-frequency domain (Park et al. 1998). First, despite the scattering observed in the records, the dispersion images are both of promising quality. Note that the surface wave energy below 3 Hz cannot be relied on as such long wavelength components (> 30 m) cannot be produced with a domain of this size. Second, in the frequency range of 3.5 Hz–12.5 Hz, both dispersion images capture the dominant velocity of the domain. Extraction of a fundamental mode dispersion curve from these images would result in almost identical selections in that range as highlighted in Fig. 6 by the automatic selection algorithm implemented by the SeisImager/SW® software. However, subtle differences can be observed between the two images at higher frequencies (> 12.5 Hz). As the dominant frequency component of the waveform approaches 12.5 Hz and higher, the dominant wavelength starts to approach the size of the heterogeneities. The MASW dispersion image subsequently trends towards higher phase velocities. This likely indicates some sort of composite velocity that accounts for the higher velocity scatters within the domain. However, the MSOR image continues to show evidence of phase velocities of approximately 50 m/s (i.e., the mean domain Vs). This observation can be explained by referring to the main difference between MASW and MSOR. MSOR creates a multichannel record from a single receiver positioned at one location while MASW averages the subsurface across the spread length. As a result, it would not be surprising to observe that MASW is more influenced by higher velocity scatters especially at higher frequencies (> 12.5 Hz).

Fig. 6
figure 6

Numerical dispersion images and corresponding fundamental modes of: a MASW; and b MSOR

Full size image

We further studied the effect of heterogeneity size on the dispersion images using a second domain with much larger (i.e., approximately triple) nominal boulder sizes (Fig. 7a). By adopting a similar approach, we produced MASW and MSOR multichannel records and their associated dispersion images (Fig. 7b, c). A comparison between the dispersion images from this model with those of the previous model demonstrates that larger heterogeneities have had a more adverse influence on the quality of the dispersion images. Still, the dispersion images acquired with MASW and MSOR approach are similar in pattern in the frequency range of 6–14 Hz despite their subtle differences. A final remark here is that although a much larger portion of the second model was filled with boulders, the MASW and MSOR dispersion images are still dominated by the velocity of the softer material and are insensitive to that of boulders. This is again because dispersion analysis of surface waves presents an average stiffness of vertically-stacked layers along a spread, and any anomalous body that is smaller than half of the penetrating wavelength (λ/2) or one-fifth of spread length (L/5) in the lateral direction will not show up as a localized entity (i.e., layer) on the dispersion image (Mi et al. 2017).

Fig. 7
figure 7

a Large scatterer model; b MASW dispersion image; and c MSOR dispersion image

Full size image

The results from both models suggest that in the absence of MASW records, an MSOR approach can be a candidate for the analysis of surface waves in the large-contrast highly-heterogeneous domain of the lunar subsurface. Consequently, this method will be used in the following section to process actual lunar seismic records.

3.2 Multichannel Apollo 16 Lunar Dispersion Images Using MSOR

Dispersion images during MASW processing are generated by applying wavefield transforms to convert multichannel recordings from the space-time domain to the phase velocity-frequency domain. An ideal MASW survey therefore requires a multichannel (e.g., 24 channels) data acquisition system to be deployed. However, as described earlier, the testing setup from the Apollo 16 mission only acquired measurements from three receivers. To mimic a multichannel survey from this dataset, we adopted the MSOR approach and generated two different multichannel records. The MSOR requirement of a consistent source was satisfied in the Apollo 16 exploration efforts given that the small explosives of the thumper were of the same charge. Figure 8 shows the first compiled MSOR record. The second multichannel record involves stacking traces of the same source-offset distance to reduce the noise and to improve the low-frequency response of the averaged record. Stacking is a common practice in seismic geophysics to overcome the adverse influences of ambient and random experimental noise (Foti et al. 2015). Multiple individual records collected with the same source-receiver offset are stacked to suppress the incoherent ambient noise and improve the low-frequency response of dispersion image. Again, there is an assumption of a layered profile with negligible lateral variation when averaging traces acquired with the same source offset, but from different shot/receiver combinations. Figure 9 presents a hypothetical example where 5 individual records collected from the same source-receiver configuration were stacked. Table 1 lists the details of the (shot, receiver) pairs used to compile the untreated and treated records in this study.

Fig. 8
figure 8

Untreated MSOR record

Full size image
Fig. 9
figure 9

Example of stacking of seismic records: [left] 5 individual shot gathers; [right] stacked signal

Full size image
Table 1 Shot and receiver numbers used to compile the multi-receiver records (dx = 4.57 m)
Full size table

After generating the multichannel records, we generated their corresponding dispersion images. Figure 10a shows the dispersion image from the untreated MSOR record (Fig. 8). Although the quality of these images is less than ideal, they are of sufficient quality to extract useful information regarding the domain and represent an improvement from the previous notion that multichannel records are unusable since they provide completely blurred overtones (e.g., Dal Moro 2015a, b). In a similar manner, the dispersion image of the stacked waveforms was produced and is shown in Fig. 10b. A comparison of Fig. 10a, b shows that while the untreated dispersion image contains promising high-frequency energy (> 9 Hz), stacking different traces have improved the low-frequency energy (~ 3 to 7 Hz). Moreover, it should be noted that similar to the numerical dispersion images, the experimental dispersion images are insensitive to interspersed rocks and boulders and only capture an average stiffness of subsurface that is dominated by the soft material. Finally, the untreated and treated dispersion images were stacked to produce a combined image with broader spectral coverage (Fig. 10c). Following the generation of the stacked dispersion image, we extracted fundamental mode dispersion curve by examining the peaks in the surface wave energy in the phase velocity–frequency domain of the combined dispersion image. Figure 10c also includes the final selection for the fundamental mode dispersion curve based on the Apollo 16 MSOR records developed in this study.

Fig. 10
figure 10

Dispersion images of: a untreated MSOR record; b treated MSOR record; and c stacked images of a and b and picked fundamental mode

Full size image

3.3 Inversion of Multichannel Apollo 16 Records

The selected dispersion curve in Fig. 10c was the reference curve for inversion. We used the stochastic neighborhood algorithm (Wathlelet et al. 2004) implemented in the Geopsy software package to run the inversion on the reference dispersion curve. The inversion algorithm module searches within a constrained parameter space for 5-layer velocity models that minimize the misfit function. During the search, the Vs of each layer was constrained to 30–500 m/s and the VP profile was linked to Vs by a Poisson’s ratio of 0.33. The layer thicknesses were free to vary from 0 to 7 m. A root mean-square error (RMSE) of less than 5% between theoretical and measured dispersion curves was considered as the threshold for acceptable models in this study.

Figure 11 plots the one-dimensional Vs profile developed for the test site. The color-coded lines represent velocity models with misfit values below 5%. Consequently, the color-coded area provides a range of acceptable velocity models all of which practically describe the observed dispersion pattern well. This accounts for the non-uniqueness of surface wave dispersion inversion and is a measure of uncertainty often employed in other surface wave studies (e.g., Griffiths et al. 2016; Teague et al. 2018; Gouveia et al. 2019). Additionally, the horizontal dashed lines represent the minimum and maximum wavelength limitations. These limits are based on the one-third wavelength approximation for penetration depth of Rayleigh waves (Hayashi 2008). Also depicted on Fig. 11 are the inverted results from a recent surface wave study (Dal Moro 2015b). Multiple velocity models were acquired in Dal Moro (2015b) for the area bounded by the two gray lines. The dashed gray line is an estimated Vs profile from the Apollo 16 Preliminary Science Report.

Fig. 11
figure 11

Inversion results of this study against those of Dal Moro (2015b). Color-coded lines represent the range of velocity models with misfits below 0.05. Dashed gray line is the estimated shear-wave velocity variation with depth from the Apollo 16 report (AA. VV. 1972). The horizontal dashed lines represent the minimum and maximum wavelength limitations

Full size image

Unlike Dal Moro (2015b) that estimates a relatively constant velocity from the surface to a depth of 7 m, acceptable models from this study exhibit a gradually increasing Vs from the surface to a depth of 7 m with the best-fit velocity model suggesting two distinctive jumps in the Vs at depths of about 2–3 m and 7 m. An increasing velocity model within this range, as suggested by the results, follows the commonplace soil mechanics knowledge that the shear modulus (G) of soils is proportional to the square root of the effective stress (σ′): \({V}_{s}\propto \sqrt{G}\) and \(G\propto \sqrt{\sigma'}\) thus \({V}_{s}\propto {\sigma' }^{0.25}\) (Seed and Idriss 1970). Although the gravity on the Moon is about 1/6 of that on Earth, it is strong enough to slightly self-compress the lunar regolith (Ogino et al. 2016). For example, if one measures a Vs of 50 m/s at a depth of 1 m, then, the increased Vs due to the weight of overburden is expected to be about 75 m/s at a depth of 5 m [i.e., 50 × (5/1)0.25]. In fact, this point was already considered by Mitchell et al. (1972) in Chapter 8 of the Apollo 16 Preliminary Science Report (AA. VV. 1972), and they estimated a gradually-increasing Vs trend for lunar regolith as shown by the dashed gray line on Fig. 11. However, the trend suggested by Mitchell et al. is merely based on an engineering judgment of lunar regolith properties rather than actual seismic measurements.

The Vs estimates from this study are also higher than those of Dal Moro (2015b) in the 2.5–7 m depth range. In addition to the discussion presented above, the discrepancies between the two inversion results can be further ascribed to two factors: (1) the use of different analysis methods; and (2) different aerial coverage. With respect to the first factor, Dal Moro (2015b) extracted the dispersive properties of the medium (in a group velocity-frequency domain) by considering a single seismic trace. The group velocity represents the velocity of a packet of waves. In this study, a multi-trace approach was adopted that characterizes the dispersive properties in the phase velocity–frequency domain. In a normally dispersive medium, it can be mathematically proven that the group velocity is lower than the phase velocity (Foti et al. 2015). This could lead to differences in inversion results. For the second factor, a single-trace inversion focuses on one-point measurement and is highly influenced by the anomalies along the source-to-receiver path. On the other hand, in this study, receiver information from different source-offset combinations were stacked and compiled together to create a multichannel record.

A final point should be made about the thickness of regolith on the velocity model in Fig. 11a. The thickness of regolith varies depending on location. Table 2 lists the some of the previous estimates of the thickness of regolith from seismic experiments as well as crater morphologies. According to this table, conflicting estimates exist regarding the regolith thickness at the Apollo 16 landing site. While seismic experiments suggest an average thickness of 12.2 m, crater morphology suggests a shallower depth in the range of 3.1–7 m for lunar regolith at this location. According to the models with the least misfit values, Vs increases from an average of 80–120 m/s at depths ranging from 6 to 7 m. Given that previous seismic studies have suggested a VP of about 250 m/s for the mega-regolith layer (AA.VV. 1972), this depth could potentially represent the average thickness of regolith at the Apollo 16 site (i.e., Vs/VP = 0.48). A Vs of 120 m/s also agrees well with the lower bound estimate for Vs of the mega-regolith layer in Dal Moro (2015b). However, a definitive conclusion can only be made by inverting a dispersion curve that contains more low-frequency energy. To achieve this, a few updates might be helpful to consider in future space explorations. First, a true multichannel spread with at least 12 receivers and shot locations on both ends provides an ideal source-receiver configuration for surface wave testing. Second, using a source with a lower predominant frequency content than the thumper, combined with major technological advancements in instrumentation would guarantee that high-resolution imaging of the lunar near-surface becomes possible in future active experiments.

Table 2 Regolith thickness from previous studies for select Apollo mission landing sites
Full size table

4 Conclusions

Traditionally, the active seismic data from the Apollo missions were analyzed using the P-wave refraction method. However, recent advancements in surface wave testing and analysis have made it possible to revisit the archived Apollo seismic data. In this regard, very few studies have focused on the use of surface wave analysis of lunar data, especially with a multichannel approach. This study summarizes the inversion of the Rayleigh surface waves from the Apollo 16 active seismic records. We extracted a fundamental mode dispersion curve from a compiled multichannel record covered the 4–15 Hz frequency range. Inversion results provided an insight into the top 8 m of lunar subsurface; the acceptable velocity models at the landing site suggest soft material with a Vs of 40–50 m/s over the top 2 m that increases in stiffness until reaches a Vs of 95–145 m/s with an average of 120 m/s at a depth of 6–7 m. Also included in this study was a numerical simulation of wave propagation in a highly heterogeneous domain that only consisted of a soft and a rock-type material. Results of the numerical modeling supported the notion that the ringy behavior of lunar seismic traces, whether active or passive, is in part due to the complex heterogeneous formation of the regolith. Future lunar active experiments can be planned with more receivers, source offset locations on both ends of receiver array, and a source with a lower frequency content to improve our understanding of the near-surface structure of the Moon.