High‐Frequency Receiver Functions With Event S1222a Reveal a Discontinuity in the Martian Shallow Crust

The shallow crustal structure of Mars records the evolutionary history of the planet, which is crucial for understanding the early Martian geological environment. Until now, seismic constraints on the Martian crust have come primarily from the receiver functions (RFs). However, analysis of the Mars RFs did not focus on the shallow structure (1–5 km) so far due to the limitation of the signal‐to‐noise ratio at high frequencies for most events. Here, we take advantage of the S1222a and six other marsquakes, which exhibit high signal‐to‐noise ratios, to probe the shallow structure of Mars. We observe a converted S‐wave at approximately 1 s after the direct P‐wave in the high‐frequency P‐wave RFs. This suggests a discontinutity at 2‐km depth between highly fractured and more coherent crustal materials.

Furthermore, under the constraints of such a three-layer seismic model beneath the lander, Wieczorek et al. (2022) estimated the average Martian crustal thickness to be between 30 and 72 km using gravity and topography data. However, previous seismic studies provided little information about the shallow structures, between 1 and 5 km, which is also important for constraining Martian geological history.
The largest marsquake since the beginning of the mission occurred on Martian sol 1222 (∼4.7 Mw), which contains broadband energy and exhibits high a signal-to-noise ratio (SNR) even at high frequencies . The location quality of S1222a is evaluated as quality A according to the analysis from the Marsquake Service (MQS) (InSight Marsquake Service, 2022), meaning that both epicentral distance and back azimuth of the event can be reliably determined. The epicenter of S1222a is approximately 37° from the InSight lander, and the back azimuth is approximately 101°   (Figure 1). We conduct P-wave RF analysis of this event because the RF method is very sensitive to velocity discontinuities. On Earth, RFs have been routinely applied to study structures from near-surface sedimentary layering (e.g., Cunningham & Lekic, 2020;Licciardi & Agostinetti, 2014;Victor et al., 2020) to the structure of the mantle transition zone (e.g., Xu et al., 2018). On Mars and the Moon, RFs mostly focus on the crustal scale so far (Durán et al., 2022;Knapmeyer-Endrun et al., 2021;Vinnik et al., 2001). The high-frequency energy of S1222a allows us to investigate the shallow crustal structures of Mars with high-frequency RFs.
In this paper, we compute the high-frequency RFs of S1222a and six other high-SNR marsquakes. We describe the data and our computational methods in Section 2. Then, in Section 3, we analyze the low-frequency (0.25-1.0 Hz) and high-frequency (0.25-2.0 Hz) RFs and propose a possible geological model beneath the InSight lander. We conclude that a newly observed converted S-wave at approximately 1 s differential travel time, extracted from the high-frequency RFs, can be related to a shallow interface at a depth of approximately 2 km.

Data and Methods
We used the 20 samples-per-second seismic data of the three very broad band components of SEIS (InSight Mars SEIS Data Service, 2019) and removed glitches (Scholz et al., 2020). After applying a symmetric taper with a Hanning window and detrending, we deconvolved the instrumental responses from the original records as the instrumental responses of the three components (UVW) are slightly different. We transformed the original records into velocity records with a high pass filter of 0.01 Hz. Subsequently, seismic waveforms of the UVW system are rotated to the vertical-north-east system (ZNE) according to the dip and azimuth information of the UVW system. Next, based on the back azimuth of S1222a (∼101°), we rotated the seismograms from the ZNE system to the vertical-radial-tangential system (ZRT) and calculated the radial P-wave RF of S1222a with the time-domain iterative deconvolution method (Ligorria & Ammon, 1999). The length of the source time function might be various for different marsquake events (InSight Marsquake Service, 2022), and there are scattered waves in the marsquake record , both of which are likely to make the RF calculated with only one time window unreliable; thus, it is necessary to test different time windows (Durán et al., 2022;Knapmeyer-Endrun et al., 2021). To determine the robustness of the crustal converted waves (0-15 s) in the RF, we computed 85 RFs for 85 different time windows (starting at the fixed 10 s before the direct P-wave and ending between 15 and 100 s, with a 1-s interval, after the direct P-wave). We finally obtained the mean RF and the corresponding ±2 standard deviations (2 ) by the bootstrap method (Efron & Tibshirani, 1993).
To test the signal's frequency content, we calculated P-wave RFs of S1222a in different frequency bands (0.25-2.0 Hz, 2.0-4.0 Hz, 4.0-6.0 Hz, and 6.0-8.0 Hz), and we found that RFs computed with frequency bands above 2.0 Hz have only one peak at zero lag-time ( Figure S1 in Supporting Information S1), which means that the horizontal and vertical components oscillate almost synchronously. This is likely caused by marsquake-induced linearly polarized vibrations of the lander or the common resonances above 2.0 Hz of either the InSight seismometer or the lander Dahmen et al., 2021;Hurst et al., 2021); thus, we chose to only use frequencies below 2.0 Hz. We used two different frequency bands. One is in the low-frequency band, where waveforms are filtered between 0.25 and 1.0 Hz, and the Gaussian width factor used to compute the RFs is set to 2. The other has a higher frequency limit, set to 2.0 Hz, and with a Gaussian width factor of 4. The Gaussian width factor acts as a low-pass filter, and half of its value approximately corresponds to the upper limit of the bandpass filtering (Ligorria & Ammon, 1999). Moreover, to confirm our observations for S1222a, we also computed RFs of other marsquakes detected before with high SNRs in the same way as for S1222a. In order to select high-SNR events, we first calculated the SNRs of seismic waveforms, filtered between 0.25-1.0 Hz, in two components (Z and R) for all events of quality A in the MQS catalog (InSight Marsquake Service, 2022) using Equation 1: where and are the root mean squares of noise and signal waveforms in a 100-s time window before and after the incident P arrival, respectively. All information on the arrival time of events is provided by the MQS (InSight Marsquake Service, 2022). We only selected high-SNR events (SNR>10.0 dB for both Z and R components) because the SNRs of waveforms might significantly influence the extraction of converted waves. Finally, we retained a total of seven events, including S1222a ( Figure 1 and Table S1 in Supporting Information S1). Five out of the selected events (S0173a, S0235b, S0809a, S1094b, and S1133c) have also already been included in previous papers (Durán et al., 2022;Joshi et al., 2022;Kim et al., 2022;Knapmeyer-Endrun et al., 2021), allowing us to confirm P-wave RF results from different studies.
Among the selected events, S0173a, S0235b, S0809a, and S1133c are close to each other, with epicentral distances between 28° and 30° and back azimuths between 74° and 91°. S1000a and S1094b are the two impact events on Mars detected recently (Posiolova et al., 2022). The epicentral distance and back azimuth of S1094b are approximately 60 and 40°, respectively, whereas S1000a is very far from the InSight lander, with an epicentral distance of approximately 125° and a back azimuth of approximately 34°   (Figure 1). For event S1000a, because of its large epicentral distance and the possible shadow zone in the Martian mantle, resulting in only weak energy before the PP wave Horleston et al., 2022), we used the PP wave to calculate RFs.
We divided events into two groups based on the seismic waveforms' SNRs at high frequencies. In Figure 2 and Figure S2 in Supporting Information S1, we computed the amplitude spectra up to 2.0 Hz of the noise and signal waveforms for the selected events. The signal energy of events S1000a, S1094b, S1133c, and S1222a is greater than their noise energy at high frequencies (>1.0 Hz). However, S0173a, S0235b, and S0809a have comparable noise and signal energy at high frequencies, resulting in a low SNR (<10.0 dB) in the high-frequency band (1.0-2.0 Hz) (Table S1 in Supporting Information S1). Therefore, based on the SNR of high-frequency signals (>1.0 Hz), we classified events S0173a, S0235b, and S0809a as type 1 events, and events S1000a, S1094b, S1133c, and S1222a as type 2 events.  (Smith et al., 2001) and is projected with the Cylindrical Equidistant method. The InSight seismometer is marked with a black triangle, and the pink stars represent the selected events. The locations of the selected events are provided by the MQS (InSight Marsquake Service, 2022).

Results and Discussion
In this section, we first computed low-frequency RFs (0.25-1.0 Hz) for the seven events and validated our approach by retrieving the three positive converted P-to-S waves detected in previous RF studies (Durán et al., 2022;Joshi et al., 2022;Kim et al., 2022;Knapmeyer-Endrun et al., 2021). Then, we computed high-frequency RFs (0.25-2.0 Hz), observed a converted wave at approximately 1 s differential travel time, and investigated which velocity models fit this ∼1-s signal. Finally, we proposed a possible model of the shallow crust beneath the InSight lander.

Analysis of Receiver Functions
We first computed P-wave RFs in the low-frequency band (0.25-1.0 Hz), following the processing described in Section 2. The 2 uncertainty of all the RFs are narrow, indicating that the RFs of all the selected events are reliable over different time windows (Figure 3a and Figure S3 in Supporting Information S1). The RFs clearly show the three positive converted P-to-S waves after the incident P-wave (between 2 and 3 s, 4 and 5 s, and 7 and 8 s, vertical orange shade boxes in Figure 3a), already highlighted in previous studies (Durán et al., 2022;Joshi et al., 2022;Kim et al., 2022;Knapmeyer-Endrun et al., 2021). The direct P-wave at the zero lag-time is weak for event S1000a, which might be attributed to the complex ray path of the PP wave of S1000a. Type 2 event: S1222a. In each panel, the blue and orange solid lines represent the amplitude spectra of the signal windows of the radial (R) and vertical (Z) components, respectively. The solid black lines denote the amplitude spectra of the noise windows of the two components. The noise and signal windows are the 100-s windows before and after the direct P-wave, respectively.
Then, we computed P-wave RFs in the high-frequency band (0.25-2.0 Hz) (Figure 3b and Figure S4 in Supporting Information S1). For three out of four type 2 events (S1000a, S1094b, and S1222a), their high-frequency P-wave RFs present a converted P-to-S wave at 1.0 ± 0.3 s after the incident P-wave, and the corresponding 2 respectively. The blue and cyan lines represent RFs of type 1 and 2 events, respectively. The gray shade in each RF stands for the 2 uncertainty. The vertical orange shaded boxes mark the arrival times of the converted phases detected before (Durán et al., 2022;Joshi et al., 2022;Kim et al., 2022;Knapmeyer-Endrun et al., 2021). The pink shade box marks the newly detected converted wave at approximately 1 s in the high-frequency RFs of type 2 events. (c) 1-D velocity profile for the lava-flow model (model 1). (d) Synthetic RFs of model 1. The blue and black solid lines are synthetic low-frequency (lower) and high-frequency (upper) RFs, respectively. (e) and (f) are the same as (c) and (d), but for the highly fractured model (model 2). The ray paths corresponding to the 1-s converted S-wave are inserted in the (c) and (e) diagrams. The blue solid and dashed lines represent the direct P-and converted S-waves, respectively. The pink arrow denotes the 1-s converted S-wave that can be recovered in the synthetic P-wave RFs. uncertainty demonstrate that this 1-s converted S-wave is reliable. Although the amplitude of the 1-s converted S-wave for S1222a is weak, that for S1000a and S1094b is strong; thus, it is not a spurious signal. Additionally, the RFs calculated with the same frequency band and Gaussian width factor by the water-level deconvolution method (Langston, 1979) confirm that the converted S-wave at approximately 1 s is not an artifact of the calculation method ( Figure S5 in Supporting Information S1). In contrast, for the high-frequency P-wave RFs of all three type 1 events (in blue in Figure 3b), there is no additional signal within 8 s other than the three converted S-waves that are already seen in the low-frequency RFs.
We also noticed that the 1-s converted S-wave is absent in the high-frequency RF of S1133c (type 2 event). Event S1133c is close to S0173a, S0235b, and S0809a, whose back azimuths are approximately 80 ± 10°, while those for S1000a, S1094b, and S1222a are approximately 34, 40, and 101°, respectively (Figure 1 and Table S1 in Supporting Information S1). There is a negative signal at ∼1 s in all of the high-frequency RFs from S0235b, S0809a, S0173a, and S1133c. Therefore, there may be a back azimuth-dependent heterogeneity that contributes to this discrepancy. However, this hypothesis needs to be investigated in the future due to the limited amount of current high-SNR data; here, we consider that the 1-s converted S-wave originates from a velocity discontinuity in the shallow crust.

RF Modeling and Geological Interpretation
In this section, we combined our results with some geological observations to explain the ∼1-s converted S-wave. InSight was deployed in the Elysium Planitia, a region covered with a wide range of basaltic lava flows. The upper 170 m of the stratigraphy at the landing site itself consists of sand-dominated (0-3 m depth) and decameter-scale (10-30 m depth) basaltic regolith underlain by ∼140 m of Hesperian-Early Amazonian basaltic lavas (hereafter: lava flow) that include an internal sedimentary unit between approximately 30 and 80 m depth (Hobiger et al., 2021). Although basaltic lavas at the landing site are <170 m thick, in the vicinity around the landing site, the thickness of the Hesperian-early Amazonian basaltic lava flows can be greater and reach thicknesses up to approximately 0.5 km or more Pan et al., 2020;Warner et al., 2017Warner et al., , 2022. Orbital imaging and infrared spectroscopy of excavated materials in nearby impact craters indicate the presence of both layered rocks and Fe-Mg phyllosilicates and suggest that layered sedimentary rocks likely extend from 170 m to depths of >1 km (Pan et al., 2020). The knowledge, like the thickness and velocity, of both the deep lava flows and sediments is not sufficient, so we explored two different models by incorporating an additional interface on top of the three crustal discontinuities of Mars from Knapmeyer-Endrun et al. (2021).
In model 1 (Figure 3c), the first interface represents a boundary between a pure lava flow and a primordial crustal material. The P-wave and S-wave velocities (Vp and Vs) are assigned to be 5 and 3 km/s, respectively, which are based on terrestrial data for basalts from Iceland that have very little crack damage and were measured and reported by Vinciguerra et al. (2005). In model 2 (Figure 3e), the first interface corresponds to a discontinuity between a highly fractured crustal material and a more coherent material, and Vp and Vs are set to 2 km/s and 1 km/s, respectively, which are well within the velocity ranges of fractured basalts and sediments at the landing site based on previous estimations (Hobiger et al., 2021;Knapmeyer-Endrun et al., 2018). We placed the two full crustal models utilized for modeling in Table S2 in Supporting Information S1.
We calculated synthetic seismic waveforms with the two models by employing the Raysum procedure (Frederiksen & Bostock, 2000) and computed low-frequency and high-frequency RFs (Figures 3d and 3f) by the time-domain iterative deconvolution method (Ligorria & Ammon, 1999). The ray parameter of synthetic seismic waveforms is set to 0.11 s/km, which is the mean of ray parameters of type 2 events in this study, and other parameters for the calculation of RFs of the synthetic seismograms are the same as those used with observed seismic waveforms. In order to fit the 1-s converted phase, the first interface of both models is set to be 2 km depth. But it should be noted that this depth may have a certain trade-off with the velocity model since we cannot well constrain the absolute velocity.
Both models can generate the ∼1-s converted wave. In model 1, the 1-s converted S-wave in the synthetic high-frequency RF corresponds to the multiple PsPs + PpSs of the first interface (Figures 3c and 3d). This is because the Vs of the lava-flow material is larger than that of the underlying crust, resulting in a negative Ps wave, while the amplitude polarity of the multiple PsPs + PpSs of the first interface is positive. In comparison, in model 2, the 1-s converted wave is the P-to-S wave of the first interface (inserted figure in Figure 3e) and can 7 of 9 be extracted by the synthetic high-frequency P-wave RF. However, the differential arrival times of the multiples of the first interface are approximately 3 s (PpPs) and 4 s (PsPs + PpSs), respectively, coinciding with that of the two crustal converted S-waves of the underlying 8-km (∼2.8 s) and 20-km (∼4.5 s) discontinuities; thus, the two multiples of the first interface cannot be distinguished. Meanwhile, for the synthetic low-frequency RF, there is a bulge following the direct P-wave, which can also be found in the observed low-frequency RFs of S1094b and S1222a (Figures 3a and 3f).
We prefer model 2 for the following three reasons. First, according to model 2 (Figure 3e), the two-way travel time between the ground surface and the first interface for the S-wave at zero offset is approximately 4 s, which falls in the first broad arrival (3-8 s) of horizontal-component seismic autocorrelation functions of InSight . In contrast, model 1 predicts a two-way travel time of approximately 1.3 s for the S-wave, which is hard to correlate with any arrivals in the horizontal autocorrelations. Second, the 1 km/s of Vs for the 2-km highly fractured crustal material is very close to the mean Vs (∼1.3 km/s) for the first 2-km depth of the velocity model obtained by Daubar et al. (2020) from the self-compaction lava flow models of volcanic material (Lesage et al., 2018). It is most likely that within the first interface of model 2, the velocity increases with increasing depth due to increased pressure and closure of fractures, which is compatible with the reduction and compression of the lava flow. And this is furthermore reinforced by the fact that the amplitude of the ∼1-s converted wave seems weaker than those of the synthetics. Last but not least, the highly fractured shallow crust of model 2, containing both fractured basalts and sediments, can closely reflect the upper crustal materials actually observed at the InSight landing site Pan et al., 2020;Warner et al., 2022). However, the 2-km thick lava flow of model 1 is significantly greater than that observed at the landing site itself. As a result, model 2 provides a better model of the shallow crust beneath the InSight lander, where there is a 2-km depth discontinuity which corresponds to a physical discontinuity instead of a chemical discontinuity (Figure 4). And the shallower crustal structure (<500 m) can not be distinguished by the P-wave RF ( Figure S6 in Supporting Information S1).
The 2-km depth discontinuity found in this study is probably heterogeneous and only limited to the InSight landing site. First, the 1-s signal is absent in RFs of events with back azimuth from 74-91°, likely hinting a heterogenous shallow crust beneath the InSight. Next, the remarkable variation in the global crustal thickness of Mars  suggests that the Martian crustal structure varies by region. Furthermore, the Martian geological map shows eight different geological units within 10 degrees of the landing site (Pan et al., 2020;Tanaka et al., 2014), indicating that the area of the Elysium Planitia near the InSight experienced a different evolutionary history. Finally, the high-resolution orbital images (meter scale)  reveal that the surface topography around the landing site varies distinctly. Therefore, it is reasonable to speculate that the 2 km thick highly fractured crustal material beneath the InSight is a heterogeneous and very local result.

Conclusions
By using seismic waveforms of event S1222a and six other high-SNR marsquakes, we obtained reliable low-frequency and high-frequency RFs beneath the InSight lander. The low-frequency RFs are consistent with previous works, while for events whose incident P-waves have high SNRs above 1.0 Hz, we observed a positive converted S-wave at approximately 1 s after the incident P-wave in the high-frequency RFs. Combined with the geological characteristics observed adjacent to the lander, this 1-s positive converted S-wave indicates a discontinuity at approximately 2 km depth, implying a highly fractured crust located above the stronger underlying crust. This discontinuity is likely heterogeneous and only limited to the InSight landing site. There is a 2-km depth discontinuity between a highly fractured crustal material and a more coherent crustal material. It should be noted that the distribution of these fractures is only a schematic and the scales of fractures are not necessarily correct.  (Hunter, 2007).