Structure of the Martian Crust Below InSight From Surface Waves and Body Waves Generated by Nearby Meteoroid Impacts

We measure group velocity dispersion of surface waves generated by two meteoroid impacts on Mars close to the lander of the InSight mission. This allows us to probe the crustal structure in the first few kilometers beneath the InSight lander. In combination with body wave arrival times from five impact events, we obtain direct seismic constraints on the seismic velocity of the crust in the vicinity of the InSight landing site. We confirm the existence of a uppermost low‐velocity layer with a mean thickness of ∼1.2 km, interpreted as layered volcanic materials, possibly interstratified with sedimentary and altered materials. Our joint inversion of surface and body waves shows a four‐layer model for the Martian crust, compatible with high‐ and low‐frequency P‐to‐S receiver functions estimated in previous studies.

• Group velocity dispersion generated by two meteoroid impacts (S0986c and S1034a) near the InSight landing site is measured • A shallow low-velocity crustal layer of ∼1.2 km thickness explains both surface wave and body wave observations • A four-layer crustal model is obtained, compatible with surface waves, body waves, and receiver function measurements

Supporting Information:
Supporting Information may be found in the online version of this article.
10.1029/2023GL104601 2 of 10 Recently, using high-frequency (0.25-2.0 Hz) RFs from three marsquakes, Shi et al. (2023) observed a converted S-wave at approximately 1 s after the direct P-wave, suggesting an additional discontinuity in the Martian shallow crust, interpreted as the transition between highly fractured and more coherent crustal materials.On a global scale, using the seismological constraints of Knapmeyer-Endrun et al. (2021) for a three-layer model beneath the InSight lander, in combination with gravity and topography data, Wieczorek et al. (2022) estimated that the average thickness of the Martian crust lies between 30 and 72 km.This estimate has recently been refined to a range of 42-56 km following the observation of multi-orbit surface wave signals from the largest seismic event recorded on Mars (Kim et al., 2023).At regional scale, further investigations of the Martian crustal properties have been conducted by Beghein et al. (2022), Kim et al. (2022), and Li et al. (2023).
In the second Martian year of the mission, eight impact events were seismically identified on Mars.Four impact events within a ∼300 km radius from of the lander were recorded (Garcia et al., 2022).They were first located with SEIS (Seismic Experiment for Internal Structure) data.Two large distant impacts were also detected (Posiolova et al., 2022) but will not be used in this study.Recently, Daubar et al. (2023) reported two additional verified impacts close to the InSight landing site.Thanks to orbital imaging, the knowledge of the crater locations, in combination with seismological analysis, is valuable information to constrain the crustal structure at the InSight landing site.
In this study, we present how we measure, for the first time, group velocity dispersion curves of surface waves generated by two meteoroid impacts close to the InSight lander.We perform a joint inversion of these dispersion curves, with P-and S-wave arrival times identified for five meteoroid impacts located at less than ∼250 km from InSight, to derive a distribution of crustal seismic models in the vicinity of the InSight lander.We then test the compatibility of our models with RFs estimated in previous studies.

Data Processing
The locations of the five impacts considered in this study are shown in Figure 1k.The body wave arrival times from small impacts have been estimated in the previous studies of Garcia et al. (2022) and Daubar et al. (2023).The P-wave was picked on the vertical component, and the S-wave on the radial component in direction of the dominant acoustic wave polarization (Garcia et al., 2022).The uncertainties of the picks strongly depend on the signal-to-noise ratio of these records, which is larger at night than during the day (Clinton et al., 2021).The body wave data set used in this study is shown in Table S1.1 in Supporting Information S1.
In the SEIS records (InSight Mars SEIS Data Service, 2019a) of S1034a and S0986c events (named according to the solar day on Mars during which the data were recorded, counted since InSight landed), we observe a significant coherence between ground velocities on the vertical and the radial components pointing in the direction of the impact, in the 0.7-1.1 Hz frequency range, presenting a 90° phase shift between these components.As shown in Figures 1c-1h (black ellipses), these features, characteristic of surface waves, are observed just after the S-wave arrival (Garcia et al., 2022).Such coherencies are often observed on Earth for Rayleigh surface waves (e.g., Panning et al., 2015) which are very likely generated by a superficial source like an impact.This surface wave signal is very clear for the S1034a event, and a similar signal is observed with a lower signal-to-noise ratio for the S0986c event.
To investigate the existence of surface waves, both signals are first bandpass filtered into narrow spectral bands (Figure S3.1 in Supporting Information S1).Various different parameterizations have been tested, and the 0.68-1.4Hz and 0.6-1.25 Hz frequency bands contain the most consistent energies for S1034a and S0986c, respectively.The envelope of each narrow-band filtered signal is then computed and the time is converted into group velocity using the epicentral distances provided in Table S1.1 in Supporting Information S1.After interpolation, for each narrow Gaussian filter frequency, the envelope gives directly the probability density function (pdf) of the seismic energy (Gaudot et al., 2020).The two resulting probability plots are shown in Figures 1i  and 1j.The probability diagram for S1034a (Figure 1i) exhibits a prominent pattern over the whole period range (0.95-1.35 s).In contrast, the diagram for S0986c is more patchy than for S1034a.This effect may be due to a scattering effect of surface waves in a diffuse crust (Menina et al., 2023), since S0986c is almost twice as far away as S1034a.It could also be caused by a less energetic impact which is consistent with the crater diameter estimations of respectively 6.1 and 9.2 m for S0986c and S1034a (Daubar et al., 2023).The region encircled by the black ellipse between 0.95 and 1.05 s in Figures 1i and 1j corresponds to large probabilities shared by the two

Inversion Method
The inverse problem consists in retrieving the seismic velocities as a function of depth in the vicinity of the InSight landing site, from surface waves and body waves measurements.This problem is, however, ill-posed: Different seismic velocity models can give similar dispersion curves and body wave arrival times.To overcome this difficulty, we solve the inverse problem with a Bayesian approach based on a Monte Carlo Markov chain method (e.g., Mosegaard & Tarantola, 1995;Tarantola, 2005).To estimate the posterior distribution of the parameters, we use the Metropolis algorithm (Hastings, 1970;Metropolis et al., 1953).
The inversion output consists of an ensemble of crustal structure models that fit the data set.For each sampled model, the synthetic group velocity dispersion curves are computed using the Computer Programs in Seismology software package (Herrmann, 2013), and we rely on the TauP toolkit (Crotwell et al., 1999) to calculate the Pand S-wave arrival times.Assuming the observations to be independent and that modeling uncertainties may be neglected, the cost function   is then defined as follows: Since the origin time of the seismic events remains unknown,   corresponds for body waves to the misfit between the observed and computed differential arrival times t S − t P , taking into account the error bars σ S−P , for the five different impact events (see values in Table S1.1 in Supporting Information S1).Superscripts throughout refer to observations (obs) and computed data (calc).For surface waves, we followed the approach initially developed by Panning et al. (2015) and developed further by Drilleau et al. (2020) and Gaudot et al. (2020).Here, the calculated group velocity dispersion curve estimated for a given velocity model is compared to the data represented by a pdf (Figures 1i and 1j).For each period i corresponding to a narrow band-pass filter B, the intersection between the calculated group velocity dispersion curve and the pdf gives a probability ρ i .The sum of 1 − ρ i over the total number of band pass filters B (5 for S0986c and 11 for S1034a) determines the misfit value for surface waves.

Prior Information
Previous studies of RFs favor a three-layer crustal model below the InSight landing site (Joshi et al., 2023;Kim et al., 2021;Knapmeyer-Endrun et al., 2021).Recently, using high-frequency RFs (0.25-2.0 Hz), Shi et al. (2023) discovered an additional discontinuity at a depth of approximately 2 km beneath the InSight lander.This finding has been confirmed by Carrasco et al. (2023) using Rayleigh wave ellipticity of large seismic events.To make our results comparable to these previous studies (see Section 4), we fixed the number of layers in the crust to four, over a half-space.The base of the fourth layer is considered to be the discontinuity between the crust and mantle (the Moho), and is allowed to vary between 31 and 47 km depth, as suggested by Knapmeyer-Endrun et al. (2021).The depth of the first three interfaces can vary freely between the surface of Mars and the Moho.V S in the four intra-crustal layers is randomly sampled within a relatively broad velocity range, between [0.5 and 4.5] km/s.Below, V S in the mantle can take values between [4.25 and 4.55] km/s, as constrained by Drilleau, Samuel, Garcia, Wieczorek, et al. (2022) and Drilleau, Samuel, Garcia, Rivoldini, et al. (2022).The range of V P /V S values allowed is 1.7-1.9.Density is deduced from a scaling relationship with V P , using Birch's law (Birch, 1961).The a priori distribution of V S is shown in Figure 2a.Because we impose that the seismic velocities within the crust increase with depth, the pdf favors larger velocity values with depth.

Results
The a posteriori distribution of V S is shown in Figure 2b.V P and density are, however, not constrained and are therefore not shown.The marginal distributions of V S and of the bottom depth of the four crustal layers are displayed in Figures 2c-2j.The prior and posterior marginal distributions are shown in blue and gray, respectively.Only posterior distribution values larger than the prior ones are considered to indicate significant information contained in the data.The crustal structure in the vicinity of the InSight landing site is characterized by a uppermost layer of a mean thickness of 1.2 ± 0.2 km and a mean V S of 1.5 ± 0.1 km/s.A second layer with a bottom depth located at 11.5 ± 3.1 km and V S = 2.4 ± 0.2 km/s is observed.The third layer shows a mean V S value of 3.8 ± 0.4 km/s, but the bottom depth is not constrained by our data set, as well as the depth and V S of the fourth layer.
Because of the small period range considered, the surface wave sensitivity kernels (Figure S4.1f in Supporting Information S1) show that the crustal structure can be probed down to ∼2.5 km depth at most with surface waves only.Except for the farthest impact (S0981c), the turning points of the body wave ray paths are located at depths shallower than 20 km (Figures S4.1a-S4.1e in Supporting Information S1), explaining the difficulty to constrain the seismic structure in the deep crust.By performing inversions of surface waves and body waves independently (Texts S.5 and S.6 in Supporting Information S1), it is obvious that the input group velocity dispersion diagrams provide significant constraints on the properties of the shallow layer, contrary to the t S − t P data.To illustrate this, we investigated the sensitivity of the group velocity dispersion curve to the shallow structure in Text S.7 in Supporting Information S1.
The fit to the data is displayed in Figures 3a, 3b and 3c-3g for surface waves and body waves, respectively.First of all, it should be noted that the inversion scheme is able to reproduce the observed distributions of group velocity with period, for both S0986c and S1034a, and the t S − t P arrival times.The contours of the dispersion curves associated with the whole set of accepted models (shown in black in Figures 3a and 3b) all match with the area of high probabilities.The marginal distributions of t S − t P (Figure 3c-3g) fit the observed differential arrival times within uncertainties (red lines), except for S0981c.This impact event is the farthest, located at 244 km from InSight.A hypothesis could be that the body waves of S0981c traveled through some heterogeneities with smaller velocity values along their ray path.We verified that performing the same inversion by removing the t S − t P data

Constraints From Receiver Functions
To assess the robustness of the results, we tested the compatibility of the crustal velocity models derived in the present study with high-frequency (HF) and low-frequency (LF) RFs measured by Shi et al. (2023) and Joshi et al. (2023), respectively.Synthetic RFs were computed for a subset of 1,010 models randomly sampled in the whole set of accepted models.Details on the measured and synthetic RFs are provided in Text S.9 in Supporting Information S1.
Comparison of the results obtained by the surface waves and body waves inversion, with results from RFs, is shown in Figure 4.For the subset of 1,010 models, the misfit value   (Equation 1) estimated by the surface waves and body waves inversion is represented in Figure 4a.The cross-correlation coefficients between the observed and synthetic HF and LF RFs are shown in Figures 4b and 4c.Considering the best models (in red), a good agreement is found for the first two layers.The arrival times and amplitudes of the main phases near 1 and 2.5 s for both HF and LF data sets, are reasonably well explained by the best models (Figures 4g and 4h).
For synthetic HF RFs whose cross-correlation coefficients with the observed RF are greater than 0.7 (in red in Figure 4g), the arrival time difference for the P-to-S wave from the ∼1.2 km discontinuity relative to the direct P-wave is ∼0.2 s, resulting in a broad direct P-wave that deviates slightly from zero lag-time in the synthetic HF RFs.A bulge at around 1 s in the synthetic HF RFs is observed, which is the multiple (PpPs) from the  1i and 1j.The output dispersion diagram is represented through contour lines enclosing 30%, 60%, and 90% of the data.Panels (c)-(g) show the marginal distributions of t S − t P .The red line and dashed lines are the observed data and uncertainties detailed in Table S1.1 in Supporting Information S1.
∼1.2 km discontinuity and is roughly consistent with the ∼1 s signal in the observed HF RFs.This case was not discussed by Shi et al. (2023), who only tested two end-member models.The slight deviation of the ∼1 s signal in the synthetic HF RFs from that of the observed HF RFs is probably due to heterogeneities in the shallow crust beneath the InSight lander (Shi et al., 2023), and to the difference in the back azimuth of the marsquakes used for the velocity model inversion and the RFs calculation.
The comparison with the LF RFs data allows us to extract models with a deep crustal structure compatible with RFs data, below the sensitivity of the surface wave and P-and S-wave data.It is evident that some of the models derived from surface wave and body wave inversion of close impacts can explain the completely independent LF RFs waveforms up to 10 s (Figure 4h), including all of the main arrivals first identified in Knapmeyer-Endrun et al. (2021).The LF RFs are compatible with discontinuities of the third and fourth layers located at ∼25 and ∼38 km depth.The correlations between the parameters of each layer are provided in Text S.10 in Supporting Information S1. Figure 4i summarizes the mean values and uncertainties of the crustal seismic velocity model estimated in our study, from the joint inversion of surface wave and P-and S-body waves, and using constraints from RFs.

Geological Interpretation
Our crustal structure corresponds fairly well with the geological knowledge around InSight, which landed on a plain of the Western Elysium Planitia, north of the Martian dichotomy boundary (Figure 1k).A simplified geological model is proposed in Figures 4d-4f.Within ∼100 km of InSight, the plain is mantled by poorly consolidated sand with limited dust cover (M.Golombek et al., 2020).Below, the plain is composed of lava flows (e.g., M. Golombek et al., 2018;Warner et al., 2022), dated to Early Hesperian by Tanaka et al. (2014) from >2 km crater count, supported by mafic mineral near-infrared spectra on rims of small impact craters (Pan et al., 2020).A regional volcanic resurfacing event occurred at ∼2.6-1.7 Ga (Early Amazonian), covering most craters below ∼2 km-in-diameter, with a thickness ranging from ∼60-140 m (Warner et al., 2022) near the landing site and up to 300 m at regional scale (M.Golombek et al., 2018).
The plain shows sinuous wrinkle ridges (e.g., M. P. Golombek et al., 2001;Maxwell et al., 1975;Watters, 1988), interpreted as compressive tectonic structures (e.g., M. Golombek & Phillips, 2009;Mangold et al., 1998;Mueller & Golombek, 2004), created during the cooling and contraction of the magmatic plume below the volcanic plain, or the whole planet.Although the formation of wrinkle ridges requires the presence of any layered upper crust, regardless of the nature of the layering, they are commonly interpreted to form in volcanic terrains composed of stacks of lava flows, based on lunar analogy.Moreover, for a few large impact craters, their internal wall shows dark layered material beyond a depth of 500 m, interpreted to be composed of damaged lava layers (Ansan, 2023).The maximum thickness of the Early/Late Hesperian volcanic material is not well constrained in the InSight landing site because of the lack of deep outcrops.
Early Amazonian lava flows embayed remnant buttes such as those located a few kilometers southwest of the S1034a impact crater.They are probably composed of the same materials as those of Nepenthes Mensae, interpreted as volcanic deposits, impact breccias, and sediments of Noachian age with aprons formed by Hesperian mass-wasted deposits (Tanaka et al., 2014).This area experienced many geologic processes associated to hydrous minerals (e.g., Carter et al., 2013), which are observed in different large impact craters such as Gale and Robert Sharp (e.g., Brossier et al., 2021;Ehlmann & Buz, 2015).We therefore cannot exclude the possibility of sedimentary rock below the landing site without being able to formally demonstrate it (Pan et al., 2020).At the interface of the Amazonian lavas on which the lander touched down and the Noachian-Hesperian Nepenthes Mensae materials, deposits of Early Amazonian-Hesperian Medusae Fossae Formation could be interlayered, due to their volcanic ash origin reworked by fluvial processes (e.g., Lefort et al., 2012;Scott & Tanaka, 1986;Zimbelman & Griffin, 2010).
From these observations in the close vicinity of landing site, at least the first 500 m of the subsurface consists of a stack of lava flows, cratered and tectonically deformed.This seems to be in agreement with the shallow 1.2 km thick layer found in this study.Our analysis cannot resolve the likely significant variation of seismic velocities with depth due to compaction and pore closure.However, our results are compatible within error bars with the compaction subsurface model proposed by Daubar et al. (2020) and Larmat et al. (2020), based from an extrapolation or the Earth lava flow model of Lesage et al. (2018).For the second deeper layer with a bottom depth at ∼11.5 km, the geologic outcrops are missing.However, magmatic materials could be good end-members as basalts, gabbros, and more felsic material, as observed in rocks scattered along the Curiosity rover in Gale Crater (Cousin et al., 2017;Sautter et al., 2015).The seismic velocities in the lower crust are compatible with basalts and gabbros, in good agreement with Drilleau, Samuel, Garcia, Wieczorek, et al. (2022) and Drilleau, Samuel, Garcia, Rivoldini, et al. (2022), who constrained the average crustal structure at larger scale from marsquake body wave measurements.

Conclusion
For the first time, group velocity dispersion diagrams generated by two meteoroid impacts located in the vicinity of the InSight landing site were measured.Considering a four-layer model, we performed a Bayesian joint inversion of the two dispersion diagrams and the differential arrival times t S − t P of five impact events, to derive the crustal structure.Our results suggest that the seismic data can be explained by a shallow low-velocity layer of 1.2 ± 0.2 km thick, with V S = 1.5 ± 0.1 km/s, consolidating the finding of Shi et al. (2023) and Carrasco et al. (2023), using a completely independent data set.We propose that this layer might be made of layered volcanic materials, possibly interstratified with sedimentary and altered materials.A second layer of more cohesive material, possibly made of basaltic materials mixed with felsic rocks, is constrained with a bottom depth located at 11.5 ± 3.1 km, and V S = 2.4 ± 0.2 km/s.To assess the robustness of our results, we verified that our four-layer models are compatible with RFs calculated in previous studies.This comparison allows us to select models compatible with all observations and to define a seismic velocity model of the whole crust below InSight (Figure 4i).Our model will help to define a robust anchoring point to derive three-dimensional crustal models of Mars, in combination with gravity and topography data (Drilleau, Samuel, Garcia, Wieczorek, et al., 2022;Drilleau, Samuel, Garcia, Rivoldini, et al., 2022;Kim et al., 2023;Wieczorek et al., 2022).

Figure 1 .
Figure 1.Recordings of S1034a (left) and S0986c (right) events.VEL Z and VEL BAZ are the vertical and radial ground velocity components, along the impact direction.Time is relative to the Mars Quake Service event start time.(a, b) Band pass filtered records between 0.9 and 1.1 Hz.Vertical dashed lines indicates arrival times of P-(gray) and S-(green) waves.(c, d) Spectrogram of VEL Z component.(e, f) Coherence between VEL Z and VEL BAZ .(g, h) Phase relation between VEL Z and VEL BAZ at maximum coherence value.Black ellipses present regions of high coherence and 90° phase shift between VEL Z and VEL BAZ , a characteristic feature of surface waves.(i, j) Group velocity dispersion probabilities of VEL Z in panels (a, b).The black ellipse indicates large probabilities shared by the two diagrams.The black dashed lines delineate the period ranges used to perform the structure inversion (see Figures 3a and 3b).The pdf is normalized by the maximum probability in each frequency band.(k) Location of the five seismic events considered in this study.Impact craters are indicated with an arrow on each HiRISE image (McEwen et al., 2007) (see Text S.2 in Supporting Information S1).The source for the S1160a is a cluster, and the two largest craters are shown.Central part shows the MOLA shade topographic map and Tanaka et al. (2014)'s geologic units in the vicinity of InSight (red square) and the meteoroid impacts (yellow points).

Figure 2 .
Figure 2. Inversion results.Panels (a) and (b) show the a priori and a posteriori pdfs of the 1-D V S profiles, respectively.Blue and red colors show small and large probabilities.Panels (c, e, g, i) are the marginal prior (in blue) and posterior (in gray) distributions of V S in the layers.Panels (d, f, h, j) are the marginal prior (in blue) and posterior (in gray) depth distributions of the discontinuities.The mean value and uncertainty are indicated in red for those parameters that are constrained by the data.

Figure 3 .
Figure 3. Fit of data for surface waves and body waves.Panels (a) and (b) show the input dispersion diagram for S1034a and S0986c, respectively.They correspond to the areas delimited by the black dashed lines in Figures 1i and 1j.The output dispersion diagram is represented through contour lines enclosing 30%, 60%, and 90% of the data.Panels (c)-(g) show the marginal distributions of t S − t P .The red line and dashed lines are the observed data and uncertainties detailed in TableS1.1 in Supporting Information S1.

Figure 4 .
Figure 4. Comparison with results from receiver function (RFs).Panels (a)-(c) show 1,010 models randomly sampled in the a posteriori distribution (Figure 2b).Panel (a) is the result from the inversion of surface waves and body waves measurements.Panels (b) and (c) display the results from high-frequency (HF) and low-frequency (LF) RFs, respectively.Panels (d)-(f) show our simplified geological model beneath the InSight lander.Panel (d) describes the crustal structure deduced from our study.Panel (e) is a zoom on the upper crustal layers.Panel (f) shows a synthetic geologic crustal model for Mars based on Knapmeyer-Endrun et al. (2022).Panels (g) and (h) present the fit to RFs waveforms for HF and LF RFs, respectively.Solid and dashed black lines show the observed data with standard deviation.The vertical shaded boxes mark the arrival times of the converted phases detected by the previous studies ofShi et al. (2023) andJoshi et al. (2023).Blue and red colors correspond to the misfit value in panel (a), and to the correlation coefficient in panels (b, c, g, h).Panel (i) summarizes the mean seismic velocity model estimated in our study.