Azimuthal Anisotropy of the Megaregolith at the Apollo 14 Landing Site

The characterization of the megaregolith on the Moon has been investigated in various ways including analysis of lunar meteorites, remote sensing of mineralogy and gravity, and deriving a seismic velocity profile. In this study, we propose a method for analyzing azimuthal anisotropy of the megaregolith. We call this method deep‐moonquake seismic interferometry applied to S‐wave coda (DMSI‐S). DMSI‐S can turn the records of deep moonquakes into recordings from virtual active sources. The retrieved virtual sources coincide with the station positions, and thus, we obtain virtual zero‐offset (pulse‐echo) measurements. DMSI‐S is applied to seven clusters of deep moonquakes recorded at the Apollo 14 landing site, resulting in virtual zero‐offset measurements at the Apollo station 14. We use the S‐wave recordings retrieved from DMSI‐S to analyze azimuthal anisotropy. We find weak anisotropy at the layer where the megaregolith is assumed to be present. We interpret our result to show that the megaregolith at this site is characterized by a layer (or layers) of impact material, following the Imbrium impact, with internal alignment of the crushed material.


Introduction
Seismically driven analyses have aided to broaden our knowledge of numerous aspects of the Moon's interior: wave propagation in the near surface (e.g., Dal Moro, 2015; Sens-Schönfelder & Larose, 2010), regional-scale Poisson's ratio (Zhao et al., 2008), mechanism of the moonquakes (e.g., Frohlich & Nakamura, 2009;Kawamura et al., 2017;Toksöz et al., 1977), and velocity profiles (e.g., Nakamura, 1983;Toksöz, 1974;Toksöz et al., 1974) including investigations which particularly focused on the crustal thickness (e.g., Chenet et al., 2006). These seismic analyses, which used moonquakes and ambient noise recorded during the National Aeronautics and Space Administration (NASA)'s Apollo missions from July 1969 to September 1977, were directly or indirectly based on the identifications of the moonquakes. The identification process was a difficult and tedious task due to the data quality of the 70s records (e.g., Nakamura, 2005). Still, a number of studies, including Nakamura et al. (1981) and Bulow et al. (2005Bulow et al. ( , 2007, made significant contributions to advance this topic. Event detection is still ongoing as can be found from a recent example by Dimech et al. (2017), who preliminary succeeded in identifying an additional~50,000 moonquakes (near-surface local events) using a pattern-recognition algorithm at the Apollo 17 site.
For all the analyses introduced above, which are all important for the eventual aim of unveiling how the Earth and the Moon were formed and are related to each other, there is still one seismological topic not yet applied. That is seismic anisotropy (Backus, 1965). While there are numerous studies of seismic anisotropy on the Earth (e.g., Bartalis et al., 2006;Crampin & Chastin, 2003;Dalton & Gaherty, 2013;Ju et al., 2019;Montagner & Tanimoto, 1991), to our knowledge, similar studies have not been conducted for the Moon. The reason behind this seemingly relates to reasons like the data quality, source-receiver configurations, and the fact that neither surface waves nor ballistic body waves are readily observed in single Apollo seismograms.
In this study, we take a step toward analyzing the Moon's seismic anisotropy, specifically in the near surface (shallower than a few kilometers) where the megaregolith is considered present. The megaregolith is a global layer of impact debris containing polymict ejecta and comminuted melt sheet generated by the numerous large impacts (Hartmann, 1973). Although the target depth is constrained by the data resource and the method we apply, this topic is still interesting because estimates of the anisotropy of the megaregolith might give additional information on whether directional alignment is present, which potentially would expand the current knowledge of the megaregolith and, thus, how the Moon has evolved.
As explained later in detail, our study location is technically fixed at the Apollo 14 landing site (3.7°S, 17.5°W). According to a review work by Hiesinger and Head (2006, and references therein), the structure of the Apollo 14 site was formed hypothetically, yet most plausibly explained to date, by the so-called Procellarum impact roughly 4.2 Ga ago, which had brought up residual material of the magma ocean (Longhi, 1977), known as KREEP (high potassium, rare earth elements, and phosphorous; e.g., Wieczorek & Phillips, 2000), from the lower crust and upper mantle. About 3.8 Ga ago, material, such as Imbrium ejecta (Fra Mauro Formation), was thrown to the Apollo 14 site by the Imbrium impact, which would explain why the Imbrium ejecta was found at this location (e.g., Jolliff et al., 2000). The thickness of the Imbrium ejecta, which can be seen as a part of the megaregolith (but not necessarily), is assumed to be not greater than 130 m (McGetchin et al., 1973). Underneath the Imbrium ejecta, the megaregolith is considered to overlay the crust (Hiesinger & Head, 2006;Hörz et al., 1991).
In order to analyze the seismic anisotropy on the Moon, we use a technique called seismic interferometry (SI; e.g., Claerbout, 1968;Schuster, 2009;Wapenaar, 2004). There are numerous ways that SI can be applied in various situations. Redatuming of seismic records is one of them (e.g., Alshuhail et al., 2012;Draganov et al., 2006), and this is what we use. It allows us to virtually change source positions to receiver (geophone) positions in completely data-driven way. Hence, it is possible to turn moonquakes (and their coda) records into virtual-shot gathers as if active sources were shot at the receiver locations. Nishitsuji et al. (2016) previously applied body-wave SI to deep moonquakes. The authors produced an estimate of the crustal thickness based on the acoustic reflection response retrieved from P-wave coda. In this study, we augment their idea by applying SI to S-wave coda. We use autocorrelation to retrieve zero-offset records (sources and receivers are co-located). For the sake of brevity, we term this technique deep-moonquakes seismic interferometry applied to S-wave coda (DMSI-S). Since we aim to calculate the azimuthal anisotropy of the subsurface, we then rotate the polarization of the retrieved reflection responses. The azimuthal analysis we perform is based on the phenomenon of S-wave splitting (Hess, 1964).
In the following, we show how DMSI-S can be applied and what we can interpret from the retrieved results.

Selection of Study Area and Data
NASA's Apollo Passive Seismic Experiment (widely known as PSE) deployed five seismometers. The seismometer at Station 11 stopped working after 21 Earth days (and we do not use these data), whereas the four remaining seismometers (Figure 1) transmitted signals to the Earth for several years, until 30 September 1977. The study regions we can choose are exactly at these stations-with DMSI-S reflection traces are estimated at station locations. From the four stations (Apollo stations 12, 14, 15, and 16), we select Apollo station 14 because of the following: (i) Seven moonquake clusters are available at horizontally close proximity to this station; (ii) the back-azimuth distribution of the clusters varies widely. A large number of events simply give more chances to increase the data quality after the stacking step of SI. The wide coverage of the back azimuth mitigates possible directional biases leading to artifacts or ghosts (e.g., Draganov et al., 2010), which are otherwise dominant especially when the ray parameters are not close to zero. According to the above two criteria, we could also select the Apollo 12 site. However, the S-phase onsets for A40 and A41 at Apollo 12 were not available in Nakamura (2005), which means that only five clusters would be available to us for utilization. Note that the Apollo 14 landing site was located at the Fra Mauro Formation (Swann et al., 1977) where relatively thin basalts are expected to be present at the surface (Head, 1975;Watkins & Kovach, 1972), but older buried layers, buried impact melt, or dike intrusions could also be present.
To visualize how different depths would lead to different steepness (or ray parameter) of the travel-time functions for different phases, we show travel times of the shallowest (depth of 747 km for A15) and the deepest (depth of 1037 km at A9) clusters as a function of epicentral distance in Figure 2. We use the velocity model proposed by Weber et al. (2011). The cyan rectangle in Figure 2 indicates a range of epicentral distances between 0 and 25°within which fall all of the clusters we use. The travel times for the body-wave arrivals can be estimated at 14.2°for A15 and 3.2°f or A9, respectively. Note that the PcP and PKiKP phases of the one-way travel-time profile shown in Nishitsuji et al. (2016) should be correctly read as PKiKP and SKiKS, respectively.
The lunar seismic data we use are fetched from the Moon Seismic Monitor where both Lunar Surface Gravimeter Experiment (LSG) and PSE data are available in ASCII format. The Moon Seismic Monitor is the interface of the space-science data archive called Data ARchives and Transmission System (DARTS), which is under operation of the Center for Science-satellite Operation and Data Archive (C-SODA) established by   Weber et al. (2011). Examples of shallower (747 km at A15) and deeper (1,037 km at A9) hypocentral depths are shown. The cyan rectangle shows the range of epicentral distances we used as shown in Table 1. Expected seismic phases are labeled along each travel-time curve.

Journal of Geophysical Research: Planets
the Japan Aerospace Exploration Agency (JAXA). The event catalogue of the deep moonquakes was summarized by Nakamura et al. (1981) with additional counts given by Bulow et al. (2005Bulow et al. ( , 2007. A summary of the cluster coordinates, together with the ray parameters, is shown in Table 1. We entirely remove short-periodic signals, oversaturated amplitude, and dead traces without compensating any of their values during our quality control before we move to the DMSI-S processing.

Deep-Moonquake Seismic Interferometry Applied to S-Wave Coda
SI allows one to redatum source positions to receiver positions. Speaking of the foundation of SI when a horizontally layered medium is considered, Claerbout (1968) showed that one can create zero-offset (pulse-echo) reflection measurement from virtual shots co-located with the receivers from the autocorrelation of the transmission response at the receivers; this was later proven by Wapenaar (2004) for any inhomogeneous elastic medium. From theirs and other theoretical foundations, many SI applications spawned for the Earth, but also for the Moon. Sens-Schönfelder and Larose (2010) applied SI to ambient noise recorded at the Apollo 17 station and successfully obtained S-wave velocity for a depth range down to 10 m. Thereafter, based on a study of Ruigrok and Wapenaar (2012), Nishitsuji et al. (2016) developed SI application to P-wave coda from the deep moonquakes.
With DMSI-S in this study, we use the S-wave coda rather than P-wave coda of the deep moonquakes for the purpose of retrieving S-wave properties. By modifying the work of Nishitsuji et al. (2016), whose theoretical foundation is based on Wapenaar (2003) and Ruigrok and Wapenaar (2012), DMSI-S can be written as is the S-wave coda associated with a deep moonquake that occurred at hypocentral location X i and was detected at receiver location X R , M i (t) is the source-time function of the ith deep moonquake, δ(t) is the Kronecker delta function, R S (X R , X R , t) is the zero-offset reflection response retrieved at location X R from a co-located virtual source emitting energy downward around the vertical, and M t ð Þ is the average of the autocorrelated M i (t) in the cluster. The symbol "*" stands for convolutional operation. The symbols t and -t denote causal and acausal time, respectively. As can be seen from equation (1), we use S-wave coda but exclude the direct S phases. This gives rise to the concern of encountering spurious events (nonphysical signals) in the retrieved reflection response due to the correlation of two primary reverberations, such as the first free-surface multiple reflections of the S phase after reverberating at some impedance contrasts. Still, we can feasibly assume that the summation in equation (1) over the different events of the cluster would lower the contribution to retrieval of such spurious responses because of the varying hypocentral depths and locations of the cluster (those ranges are shown in Figure 1 and Table 1  First, based on the event catalogue (e.g., Nakamura et al., 1981), we fetch lunar seismograms of the selected deep moonquakes for the vertical and the two horizontal components recorded by the Apollo LP seismometer. Note that an amplitude difference was found between the two horizontal components (Garcia et al., 2011). We proceed and apply the processing to both data without and with gain correction to study the effect of this correction. For the gain correction, we use a value of 0.748 taken from Garcia et al. (2011). We then apply to the obtained data a band-pass filter of 0.2-2.0 Hz (Nakamura, 2005), after removing the instrument response (Latham et al., 1969). Since the signals of the deep moonquakes are often weak, the above filtering is expected to increase the signal-to-noise ratio (Nakamura, 2005). Another point of using this bandwidth is to minimize the effect of surface waves whose dominant frequency is around 4-12 Hz (Larose et al., 2005). While Bulow et al. (2005) and Nakamura (2005) stacked individual traces to boost seismic phases for the identification of events, one of the reasons we stack traces is to suppress possible remnants of the surface waves in the DMSI-S operation (equation (1)). Further explanations can be found in Nishitsuji et al. (2016).
We show an example of a deep-moonquake seismogram for the three components (one vertical and two horizontal) in Figure 3. For the arrival times of the (direct) S phases, we refer to Nakamura (2005) who identified them by stacking individual traces. Since DMSI-S uses the S-wave coda part without the direct S phases, we add 5% of extra time to the estimated S-wave arrival times and extract a time window of 33 min (explained further below) starting after the added 5% extra time. Latham et al. (1970) and Nakamura (1977) suggested that the coda could be recorded for up to 60 min due to the low attenuation and highly scattering environment of the Moon, while Blanchette-Guertin et al. (2012) indicated that the lunar coda lasts for 30 to 60 min after the S-phase arrival.
In addition to the two horizontal components, we also make use of the vertical component but only for quality-control purpose of traces when we choose events from the catalogue. Note that NASA reported that the vertical component at the Apollo station 14 has shown instability. We, however, do not observe response  (Nakamura et al., 1981). Three components, one vertical and two horizontal ones, are shown. The light-yellow rectangle shows the range of S-wave coda used in this study, whereas the magenta rectangle-5% of additional time from the S-phase arrival-indicates the duration we exclude from our analysis. The arrival times of the P phase and the S phase are after Nakamura et al. (1981) and Nakamura (2005), respectively.

Journal of Geophysical Research: Planets
instability of the vertical component in the traces we use. As reported in Latham et al. (1971), this instability can be corrected for by removing the feedback filter.

DMSI-S Application
Having extracted the S-wave coda from each selected event, we subsequently normalize the individual amplitudes per cluster to equally treat their magnitudes; that is, in the summation process in equation (1) all moonquakes will contribute to the final result in a comparable way. Then, we apply DMSI-S to each cluster. Note that phase conversions from S to P waves are not considerable because we use for DMSI-S deep moonquakes with first arrivals characterized by small ray parameter. Moreover, the use of horizontal components further attenuates the presence of (scattered) subvertical P-waves. To see how DMSI-S results differ when we change the coda duration, we conduct a simple experiment as shown in Figure 4. Here, we compare the results by using a duration of 200, 500, 1,000, and 1,980 s (the one we eventually select). These durations correspond to the bars shown in Figure 3. Qualitatively speaking, the shorter the duration we use, the lower the signal-to-noise ratio in Figure 4, especially for times later than 7 s. In other words, the effects after DMSI-S at earlier times (e.g., <7 s) are considered to be rather marginal. According to Blanchette-Guertin et al. (2012), we might further increase the duration we use (e.g., 60 min). However, the number of events we can use in practice generally begins to decrease as the probability of encountering unwanted events (e.g., meteorite impacts) and other unwanted phenomena (e.g., amplitude saturation due to mechanical trouble) increases for those extended portions. Our choice of 33 min, which we obtain by taking the minimum S-wave coda duration (Blanchette-Guertin et al., 2012) and prolonging it by 10%, is dictated by our observation that within such time lengths we do not encounter such unwanted phenomena. Although  Table 1), the responses we actually obtain, for example, in Figure 4d, turn out to be visually similar. As we explain later, stacking of individual DMSI-S results computed from the separate clusters per station should lead to stabilization of the DMSI-S results around the stationary phase of the reflection response.

Azimuthal Rotation
Using the DMSI-S results, we proceed to calculate the azimuthal anisotropy. To do that, we rotate the moonquake records to obtain horizontal particle oscillation between 0 and 360°, with increments of 5°. Subsequently, we apply DMSI-S to estimate a reflection response for these different horizontal polarizations. We obtain seven azimuthal rotation DMSI-S results from the seven clusters (Table 1). Based on the event catalogue, none of the clusters we use should provide a ray parameter of exactly zero except after long scattering (Table 1 and Figure 2). To mitigate possible biases due to the horizontal offset between Apollo 14 and the separate clusters (i.e., cases that the retrieval process might result in erroneous reflections) through the DMSI-S processing, we stack autocorrelations from those clusters (hence events therein). Such directional biases (i.e., sources are not equally distributed in the space domain) would lead to retrieval of artifacts (nonphysical events) in the reflection responses (e.g., Melo et al., 2013). However, as Draganov et al. (2010) showed for instance, the stacking is an effective process for suppressing such artifacts. Here, making use of the advantageous wide back-azimuthal distribution of the seven clusters with respect to Apollo 14 (see Figure 1), we proceed to stack the retrieved results from all clusters to average out potential directional biases. The stacked result of the seven clusters is shown in Figure 5 after muting the retrieved band-limited delta pulse at t = 0 (equation (1)). Note that the traces for 180-360°, which are symmetric to the 0-180°traces, are shown for visualization purpose. As can be seen in Figure 5, the difference with and without the gain correction is small.

Azimuthal Anisotropy
The degree of seismic anisotropy in an elastic medium can be estimated through analyses of the seismic-velocity variation as a function of azimuth. Backus (1965), as a pioneering work, derived an equation of the azimuthal anisotropy considering the mantle of the Earth. In his paper, the solutions for both P and S waves are given. Later, Crampin (1977) reviewed a number of studies on seismic anisotropy including Backus (1965) and his other papers and further introduced the detailed solution for S waves where SH and SV waves were separately written.
In our study of the Moon, we do not know whether the propagation of the S phases is expected to be dominated by the SH phases, SV phases, or a combination of them. Thus, we assume that the DMSI-S results ( Figure 5) contain a combination of SH and SV phases. For this, we combine equations (13) and (14) from Crampin (1977) for the variation of the SH and SV phases, leading to the following equation: where u are unknowns to be fitted with the rotated DMSI-S results and θ is the polarization direction (azimuth). Using equation (2), we estimate the degree of the azimuthal anisotropy. Tonegawa et al. (2015) described S-wave anisotropy on the Earth using a processing flow similar to ours. Also, their synthetic test (Case 2 in their Discussion) demonstrates that the estimation of the anisotropy with this processing is robust.
Note that velocity information is not needed for estimating the anisotropy per se. However, we use the velocity model proposed by Weber et al. (2011) to obtain a rough estimation of the depth of the megaregolith. With their velocity model, which assumes the megaregolith present at the depth of 1 km (e.g., Cooper et al., 1974;Weber et al., 2011), the reflection from the base of the megaregolith should arrive at about 4 s in Figures 4 and 5, while the largest peak in Figure 5 around 1 s could be a part of the source-time function of the retrieved virtual source. Since we do not intend to evaluate the sources, hereafter we focus on the second largest peak in Figure 5 around 5 s, which we interpret as the primary reflection from the base of the megaregolith at the scale of our vertical resolution (about 0.8 km, that is, a quarter of the dominant wavelength at 1 Hz). We use this reflection for further anisotropy analysis. To visualize what we mean by the base of the megaregolith at the seismic scale of our resolution, we sketch a highly schematic diagram of the megaregolith and the surrounding subsurface layers in Figure 6. Note that for our vertical resolution, features 10.1029/2019JE006126

Journal of Geophysical Research: Planets
shorter than about 0.8 km will not be captured as separate structures but would rather appear as one reflection arrival due to the tuning effect (Kallweit & Wood, 1982). We discuss interpretation further below. (2) is used to fit the arrival time of the peak in Figure 5 around 5 s. The fitting appearance is quite good as witnessed by the small resulting residuals (Figures 7b  and 7d, respectively). To qualitatively evaluate the fitting, we use the sum of the squared errors. For  Figures 7b and 7d, the sum of the squared errors of the fitting is 3.4 * 10 −4 and 6.7 * 10 −4 , respectively. Using equation (2), we estimate the azimuthal anisotropy at Apollo station 14 from the data with and without the gain correction to be 0.93%. The value is expressed as a peak-to-peak percentage calculated by subtracting the fastest time from the slowest time and dividing the result by the average time. To evaluate the fitting robustness in Figure 7a, we alternatively perform a bootstrapping test. After randomly selecting four out of the seven clusters recorded at Apollo station 14, DMSI-S is repeatedly carried out. The four  Figure 8. In Table 2, we also summarize the fitting parameters for both the full and the bootstrapping test. Note that the fluctuations (standard deviation) in relation to the full stack are calculated based on the bootstrapping test.

Figures 7a and 7c show our results when equation
6. Discussion

Validation of the Azimuthal Anisotropy
Various studies have investigated azimuthal anisotropy on the Earth (e.g., Backus, 1965;Crampin, 1977;Wild & Crampin, 1991). The value ranges of the azimuthal anisotropy for the polarized S-wave velocity depend on the depth and location. For example, Yao et al. (2010) found the azimuthal anisotropy to be 0-3% at the depth of 0-25 km in southeast Tibet, Dalton & Gaherty (2013) found it to be 4-5% in the continental crust of northwestern Canada, and Tonegawa et al. (2015) estimated it to be 1-2% at the marine sediments of the northwestern Pacific. In contrast, we know little about it for the Moon. The reason for this is most likely the quality of the available data, which limits the application of the techniques used regularly for estimation of the azimuthal anisotropy of the Earth.
The azimuthal variation in two-way travel time of the retrieved reflection we use is about 0.05 s. To quantify the uncertainty in the retrieval of the reflection, we perform a bootstrapping test. We select 148 traces (equal to the total number of all traces from all clusters) randomly from a pool containing all traces collected from all clusters (allowing repetition of traces) and apply to them DMSI-S. We repeat this selection and retrieval 100 times. The bootstrapping is performed for one polarization-0°. The standard deviation in the two-way travel time we obtain from the test is 0.015 s. Thus, the azimuthal variation of about 0.05 s in the retrieved reflection is more than 3 times higher than the uncertainty in the retrieved reflection of 0.015 s.
With the DMSI-S results, we estimate that the azimuthal anisotropy at Apollo station 14 is 0.93% when the full-stack result is used ( Table 2). The uncertainty of the azimuthal anisotropy can be ±0.32, which we define to be equal to the standard deviation based on the bootstrapping results ( Figure 8). The fitting features like peaks around 15°and troughs around 135°can be found in both Figures 7 and 8, which suggests that the fitting is qualitatively robust.
However, the bootstrapping results in Table 2 indicate that the range of parameters using equation (2), where both SV and SH phases are assumed to be present as we explain earlier, cannot be estimated stably -the standard deviation for some of the fitting parameters in Table 2 is comparable to the corresponding full-stack parameter. This indicates that our assumption over equation (2) might not be reflecting the reality. Instead, the presence of the sagittal symmetry in Figures 7 and 8 could be indicating that the wavefield we

Journal of Geophysical Research: Planets
use is dominated by the SV phase with relatively small dependence on the SH phase. In this case, a more appropriate fitting, in terms of stability, is to use equation (17) in Crampin (1977): Note that u 3 is not in Crampin (1977), but we include it in our relation to account for the unknown symmetry axis for our data. Figure 9 shows the full-stack results when we fit the DMSI-S results ( Figure 5) using equation (3) with and without the gain correction. Figure 10 shows the bootstrapping tests corresponding to the results in Figure 9a. The bootstrapping method we use for the result in Figure 10 is identical to that we use for the results in Figure 8. Qualitatively, the fitting with equation (3) (Figures 9 and 10) does not improve on the fitting with equation (2) (Figures 7 and 8). However, the parameter variations per bootstrapping in Table 3 become less scattered, suggesting that fitting using equation (3) is more stable. The estimated anisotropy is 0.75% ± 0.32, which is lower than when using equation (2). Without using the fitting curves, we obtain identical value when equation (2) is used.

Interpretation
While our DMSI-S technique is able to generate pseudo seismic experiments by utilizing the natural moonquakes, locations are inherently limited by the available stations. In this study, the location is technically fixed at the Apollo 14 landing site. Applying the velocity model of Weber et al. (2011) to the DMSI-S  (2) is used for the second-largest peak in Figure 5. (b) Fitting residuals of Figure 7a. (c) As for (a) but for the data with gain correction. (d) As for (b) but for the data with gain correction. The fitting parameters are given in Table 2. 10.1029/2019JE006126

Journal of Geophysical Research: Planets
results shown in Figure 5, we should be looking at depths between 0 and 1.25 km where the megaregolith is present (Tong & Garcia, 2015). According to Hartmann (2003), the thickness of the megaregolith is around 3 km (not necessarily at the Apollo 14 site, but in general), and the crust, which could be highly fractured (Hiesinger & Head, 2006), is beneath it.
The retrieved reflection we interpret should be seen as representing a subsurface feature or subsurface features around the depth of 1.25 km. Due to our vertical resolution of 0.8 km, multiple features around this depth would not be seen as separate reflection arrivals, but as a single tuned arrival with longer apparent time function (Kallweit & Wood, 1982). Therefore, we are limited by this scale in our discussion of the characteristics of the reflection and azimuthal anisotropy. Considering the general thickness of the megaregolith, Figure 8. Results of bootstrapping tests for Figure 7. Four out of seven clusters are randomly selected; the selection is repeated four times. The fitting parameters are given in Table 2. Note. The bootstrapping is carried out four times. The ranges next to the full stacks stand for standard deviation estimated from the four bootstrappings.

10.1029/2019JE006126
Journal of Geophysical Research: Planets the reflection depth around 1.25 km should represent part of the megaregolith or be close around its base. Nonetheless, it is difficult to ascertain exactly at which part of the megaregolith we are looking especially due to the possible gradational change of feature as a function of depth. The distinct reflector we see can be explained by the following possibilities: (i) an impact melt layer, (ii) a buried lava flow, (iii) mafic intrusion, (iv) a deep-seated fault offset, and (v) rotated impact block following the collapse of an impact basin.
According to Williams and Zuber (1998) and Dibb and Kiefer (2015), the basalt's thickness at the center of the Imbrium Basin is around 5 km, which should be deep enough to affect the megaregolith at the Apollo 14 site. Assuming that the impact (or successive impacts) hit and crushed the ground, which resulted in delivering the ejecta in the radial direction away from the impact center in a wide range of azimuths (Schultz & Crawford, 2016), we expect to detect azimuthal anisotropy from the results. The estimated local azimuthal anisotropy could result from the change of the thickness of the layer of ejected and deposited crushed material in the radial direction. Our resolution, though, might not allow detecting such a change alone. Another reason for the azimuthal anisotropy could be alignment of the ejected crushed material in the radial direction inside the deposited ejecta layer (or layers). We would detect such material alignment along the line connecting Apollo 14 with the center of the Imbrium basin. We are aware that our results do not have sufficient coverage to generalize the characterization of the megaregolith around this region. The characteristics away from the Apollo 14 site might be quite different from our observations. However, we interpret that a high-energy provider like the Imbrium impact should profoundly affect the subsurface characteristics on regional scale. Hence, the weak anisotropy of the reflection we analyze might indicate that the layer or  Table 3.

10.1029/2019JE006126
Journal of Geophysical Research: Planets layers inside the impact block do show some weak alignment of their forming material (e.g., crushed pieces of original rock, as opposed to alignment of the building minerals of a rock type), which is a direct result of the Imbrium impact. The detected weak anisotropy might also be partly caused by the change of the thickness of the ejecta layer (or layers) in the radial direction.
Assuming that the layer(s) of the impact block are situated around the depth of 1.25 km of the megaregolith, such blocky material could lead to relatively lower density due to extra spaces between blocks. In this case, the scenario of the impact block also agrees with the low gravity value observed around the Apollo 14 site from the Gravity Recovery and Interior Laboratory (GRAIL) mission by NASA (Besserer et al., 2014;Wieczorek et al., 2013).
We cannot exclude other possibilities giving rise to the retrieved reflection response, such as a layer of melt or lava flow. Sheet-like features would mean more directional alignment of material, which we do not  (3). The fitting parameters are given in Table 3. observe in terms of the anisotropy. On the other hand, the presence of a magmatic intrusion or melt around a certain azimuthal direction might be responsible for the estimated weak anisotropy. Yet another, and even more probable, case is that the material characterizing the retrieved reflection to be a combination of the above-mentioned possibilities which we cannot simply measure separately.
Although our study does not provide spatial information away from the Apollo 14 location, knowing the characteristic of the megaregolith from various sites in future could give a useful clue for constraining lunar evolution models such as the asymmetric floatation crust (e.g., Arai et al., 2008), the serial magmatism (e.g., Dygert et al., 2017), and the impact modification (e.g., Piskorz & Stevenson, 2014). Also, our proposed method can be applied for other planetary missions, for instance, Mars Insight if detectable natural quakes (e.g., Giardini et al., 2020) are present underneath.

Conclusions
We proposed a method that gives a new perspective of seismic azimuthal anisotropy on the Moon. The method-DMSI-S-is applied to characterize the megaregolith at the Apollo 14 landing site. Based on the small value of the azimuthal anisotropy we estimated by DMSI-S, the presence of a horizontally layered impact block following the Imbrium impact is implied. Since DMSI-S can be applied to any natural quakes when a suitable source-receiver configuration is available, DMSI-S has a potential to contribute as a tool for characterization of planetary geology.