Constraints on Midmantle Discontinuities in Mars from Analysis of InSight Seismic Data Open Access

Throughout the 1440 Sols (a Sol is a Martian day and corresponds to ∼24 hr 40 min) that the InSight mission was active on the surface of Mars, more than 2700 seismic events, including marsquakes and meteorite impacts, were detected and cataloged (Ceylan et al., 2022; InSight Marsquake Service, 2023). Despite significant improvement in our knowledge of Mars’ interior structure, including observations of intracrustal layering, Moho, a thick lithosphere, a mantle that appears to be equivalent to Earth’s upper mantle from a mineralogical point of view, and an entirely liquid core that makes up about 50% of its radius, there are fewer constraints on the structure of the mid‐to‐lower mantle (depth range ∼1000–1700 km) (e.g., Drilleau et al., 2022; Durán, Khan, Ceylan, Charalambous, et al., 2022; Durán, Khan, Ceylan, Zenhäusern, et al., 2022; Lognonné et al., 2023). This is in part due to the geographical distribution of the low‐magnitude seismic events (⁠$Mw∼2.5–4.6$⁠) with generally poor signal‐to‐noise ratio (Ceylan et al., 2023; Kawamura et al., 2023) that have been recorded and in part because of intense near‐receiver scattering presenting observational challenges to detecting seismic phases that have traversed Mars’ mid‐to‐lower mantle (Durán et al., 2022).

Of importance about mantle structure is the depth to and sharpness of Mars’ olivine transition seismic discontinuities (MOTDs). On Earth, these result from the transformation of olivine to its higher‐pressure polymorphs, wadsleyite, and ringwoodite, at depths of about 410 and 660 km, respectively, that bound the mantle transition zone (e.g., Helffrich, 2000). Because of the Mars’ smaller size and its relatively large core, however, only the discontinuity associated with the transformation of olivine, which occurs somewhere between 12 and 14 GPa, corresponding to depths of 1000–1200 km, is expected to be observed, as illustrated in Figure S1, available in the supplemental material to this article. Given the pressure, temperature, and compositional dependence of this mineral transformation, imaging an MOTD holds the potential of providing insights into the thermochemical structure of the mantle. The depth, and hence pressure, at which the discontinuity could be observed would give a temperature constraint at that depth because it is primarily influenced by mantle temperature (Munch et al., 2018), although composition may have a greater effect on the depth of this transition at higher temperatures (Khan, Liebske, et al., 2021). Projecting that temperature upward and downward to the boundaries of the mantle, the base of the lithosphere and the core‐mantle boundary, would provide important geophysical constraints on for example, the style of convection in Mars (Cheng et al., 2024; Murphy and King, 2024), the absence of a present‐day magnetic field (Mittelholz et al., 2020), and the flexure of Mars’s lithosphere (Broquet et al., 2020).

Seismic phases sensitive to midmantle structure include body‐wave conversions at, triplications from, and underside reflections off the MOTDs. Huang et al. (2022) recently reported on the seismic detection of a Martian MOTD at a depth of about 1000 km through analysis of triplicated waveforms. Moreover, Deng and Levander (2023), using single‐station autocorrelation analysis, reported observation of a reflection response that was interpreted to emanate from the Martian olivine–wadsleyite transition. However, where Huang et al. (2022) only had five teleseismic events at their disposal, there are now 13 events available at the appropriate epicentral distance range for detecting seismic phases that have interacted with MOTDs. Building on the earlier results and given the larger number of events at our disposal, we report on a complete analysis of the Insight event data that have P‐ and S‐wave turning points between ∼1000 and 1200 km depth and, as a consequence, have spent part of their trajectory in the MOTD. Following this, we realign, relocate, and subsequently determine more precise event back azimuths in an attempt to improve the detection of MOTD‐interacting phases. Finally, we compare the observed data with predictions based on the most up‐to‐date models of Mars’ interior structure and attempt to infer the characteristics of the Martian MOTD based on the available evidence presented here.

As a result of the increase in seismic wave velocities associated with a Martian MOTD, three arrivals corresponding to direct, reflected, and refracted P and S waves are predicted to be observed within a short time window at distances between ∼65° and 85°. This phenomenon is referred to as “triplication” and is illustrated in Figure 1 based on synthetic P‐ and S‐waveform sections and ray paths for a representative Martian seismic velocity model. The Mars models used herein are based on a joint geophysical–mineral physics parameterization (e.g., Khan, Ceylan, et al., 2021; Durán et al., 2022) that incorporates the olivine transition and are obtained from the most up‐to‐date differential body‐wave travel time and astronomic–geodetic data set (Khan et al., 2023) (see Data and Resources for model availability). The high‐frequency (HF) (1 Hz) waveforms were computed using AxiSEM by Nissen‐Meyer et al. (2014) for a source located at 40 km depth and dominated by normal faulting (Jacob et al., 2022). Although this source may not represent the focal mechanism of events outside of Cerberus Fossae, which is the most seismically active region seen by InSight (Stähler et al., 2022), it serves to capture the main characteristics of Martian triplicated waves.

More generally, detection of separate triplicated waveform branches, as illustrated in Figure 1, depends not only on the velocity structure but also on the duration of the source time function and the frequency range under consideration (e.g., Bissig et al., 2022). Triplicated S waves exhibit longer differential travel times than P waves because S waves travel at slower velocities than P waves. Consequently, where triplicated P waves are expected to arrive within ∼10 s of the first arrival, the time window for triplicated S waves increases to ∼20 s after the main S‐wave arrival. Figure 1 also shows that where P waves exhibit a well‐defined triplicated waveform pattern, amplitudes of triplicated S waves potentially suffer from interference of other seismic phases, including the depth phases and associated triplications. Finally for the Martian model considered here, interference arises from negative upper mantle S‐wave velocity gradients, which act to produce an additional late triplicated branch (yellow ray paths in Fig. 1d,e), corresponding to shallow propagation.

The Marsquake Service (MQS) has classified all seismic events recorded by InSight according to frequency content as either low‐frequency (LF), HF, or wideband (WB) events (Ceylan et al., 2022). The LF family comprises events with energy mainly below 1 Hz and mostly clear P‐ and S‐wave arrivals with an overall duration of ∼10–30 min. The HF family consists of events with predominant energy above 1 Hz (up to ∼12 Hz) and with a duration ranging from 5 to 20 min. The WB events, on the other hand, include events with energy spanning from well below 1 Hz to well above 4 Hz. The WB class also includes some events previously classified as LF, for example, the imaged meteoroid impacts (Posiolova et al., 2022) and the largest event detected during the mission (Kawamura et al., 2023), but also marsquakes previously grouped into the HF family that were considered to be trapped waves propagating in the Martian crust (van Driel et al., 2021).

Of importance in the analysis of Martian seismic data is removal of nonseismic signals related to atmospheric perturbations, transient signals (e.g., glitches), lander reverberations, subsurface resonances, and instrument artifacts (Ceylan et al., 2021). To limit their impact, MQS have produced a complete catalog of denoised event waveforms through application of deep learning models (see Dahmen et al., 2024 for details).

Events are labeled by the InSight team by mission Sol of occurrence and sublabeled alphabetically for Sols with more than one event. The events considered here and by Huang et al. (2022), including MQS epicentral distance estimations and quality labels (InSight Marsquake Service, 2023) are summarized in Table 1. Of these, only two events are quality A, characterized by high signal‐to‐noise ratio (SNR), clearly distinguishable P‐ and S‐wave arrivals, and determinable back azimuths. The remaining events are quality B and C, with lower SNR and partly indistinct phase arrivals. Of the five events analyzed by Huang et al. (2022), we excluded S0345a (quality C), partly because of unclear phase arrivals and partly because its waveform envelopes are similar to those of more distant events (Ceylan et al., 2023). Significant deviations also exist for S0234c (quality C). Where we find an S−P differential travel time of 329.6 s based on the denoised data, Huang et al. (2022) uses 395 s; yet, the P wave is unobservable in the “original” data. In addition, for impact event S1094b, Huang et al. (2022) used a back azimuth of 35.9° ± 6.5° in place of 51.4° derived from orbital images of the impact (Posiolova et al., 2022). The WB class events, meanwhile, are characterized by relatively large uncertainties in arrival times and therefore in epicentral distances. Thus, initial refinement of the main P‐ and S‐wave arrivals is required prior to searching for triplicated phases.

For this, we follow Durán et al. (2022) and use raw and polarized waveforms and equivalent time‐domain envelopes, in addition to three‐component scalograms to enhance body‐wave arrivals in the low magnitude, noisy marsquake signals. The procedure is illustrated for WB event S1153a in Figure 2. Equivalent figures of denoised and nondenoised waveforms for all 13 events considered here are shown in Figures S2–S13. Although denoised waveforms seemingly provide clearer phase arrivals, the removal of important signal has been noticed in some cases. Consequently, we simultaneously analyze both denoised and nondenoised waveforms throughout. The scalogram indicates the arrival of energy from well below 0.1 to 1 Hz and above.

A zoom into the band‐pass‐filtered denoised waveforms and their time‐domain envelopes are shown in Figure 2 (Panel I) for time windows containing the P‐ and S‐wave arrivals, respectively. To increase the SNR of body waves, we follow Durán et al. (2022) and apply a time‐domain polarization filter that exploits the linear polarization characteristic of body waves. Through principal‐component analysis, this method enhances the linearly polarized components in the signal as indicated in the filtered polarized vertical (Z) and horizontal (N) component waveforms for event S1153a. Because P‐ and S‐wave arrivals are distinguished not only by their strong linear polarization, but also by their near‐vertical and near‐horizontal inclinations (for a close‐to‐vertical incidence angle), respectively, we also perform a three‐component polarization analysis and use TwistPy (see Data and Resources) to filter in the time–frequency domain. This filtering workflow effectively identifies and suppresses (or preserves) features in the waveforms by making use of the polarization attributes of the signal and has been applied to remove glitches (Durán et al., 2024) and enhance body waves (Brinkman et al., 2023). To further strengthen body‐wave arrivals, we apply a filter mask that preserves signals with ellipticity values below 0.5 and inclinations between 0° and 30° for P waves and 60° and 90° for S waves, measured from vertical. The polarized and polarization‐filtered (after application of TwistPy) waveforms exhibit a series of high‐amplitude packages of energy associated with body waves.

A complementary method for identifying seismic phases involves the use of narrow‐band‐filtered time‐domain waveforms and their envelopes, labeled filter banks (e.g., Khan, Ceylan, et al., 2021; Durán et al., 2022). To this end, we filter the velocity traces in frequency bands half an octave wide, centered on frequencies ranging from 1/4 to 1 Hz, which captures the signal’s LF energy. Waveform filter banks and time‐domain envelopes, along with raw, polarized, and polarization‐filtered waveforms (only envelopes are displayed for the latter) are shown in Figure 2 (Panel II) for the denoised data set. The onset of energy associated with both P‐ and S‐wave arrivals is apparent in the envelopes across most frequencies and in the raw, polarized, and polarization‐filtered waveforms.

For a phase to be identified as a body‐wave arrival, it must conform to the following set of criteria: (1) be present across different frequency bands; (2) be present in both, raw and polarized traces/envelopes; and (3) differential travel time between phase pairs follows a sequence for example, PPP−PP < PP−P. In addition to these criteria, we also require that body‐wave arrivals must be consistent in both denoised and nondenoised waveforms, as well as in the polarization‐filtered waveforms. Following these criteria, we selected the P‐ and S‐wave arrivals indicated by vertical green lines in Figure 2 (Panel II), whereas arrival‐time uncertainties are represented by green bars. A similar analysis has been applied to all events (see Figs. S2–S13).

Triplicated S‐wave identification typically relies on transverse (T)‐component waveforms, which require estimating the back azimuth. For this, we use particle‐motion estimates from time–frequency polarization analysis and hodograms based on the work of Zenhäusern et al. (2022). The time–frequency polarization analysis relies on transforming the data into the time–frequency domain using continuous wavelet transformation and then performing eigen‐analysis on the resulting spectral matrix. The instantaneous back azimuth is determined by projecting the semimajor axis of the best‐fitting polarization ellipse onto the horizontal plane. Because this axis is more stable for rectilinear signals, it provides a reliable constraint for P‐ and S‐wave arrivals.

The estimated azimuth as a function of time and frequency for the entire waveform is shown in Figure 2c (Panel III). To determine azimuth based on P and S waves, we computed azimuthal probability densities for the time windows centered on the P‐ and S‐wave arrival picks, widths for which are determined by the assigned arrival‐time uncertainty, for frequencies ranging between 0.3 and 0.7 Hz and 0.5 and 0.9 Hz, respectively. The maximum of the P‐wave azimuthal probability density curve (not shown for brevity) is considered representative of the back azimuth of the event, whereas the S wave is used for corroboration (see Zenhäusern et al., 2022 for details). The so‐obtained back azimuths are indicated in Figure 2m,o (Panel III, orange solid and black line).

To validate the azimuths, we compute P‐ and S‐wave horizontal particle motion hodograms, using the same time and frequency range as in the polarization analysis. The filtered waveforms, centered on the main body‐wave arrivals, along with the corresponding hodograms, are displayed in Figure 2m–o (Panel III) and are seen to align with the direction determined from the time–frequency polarization analysis. Similar polarization analysis is conducted for all events (see Figs. S2–S13). Table 1 summarizes estimated back azimuths for all events. To each event, we have also assigned a quality label based on the reliability of the estimated back azimuth. For comparison, MQS only determined the back azimuth for S1102a.

Our final set of body‐wave picks and corresponding uncertainties of P‐ and S‐wave arrivals for the 13 events considered in this study are listed in Table 1, whereas the realigned envelopes, along with their updated differential travel‐time uncertainties, are summarized in Figure 3. Although uncertainty in differential S−P time, which is a proxy for epicentral distance, is significantly reduced for most body waves compared to those of MQS, some, particularly the weaker low‐frequency P waves (e.g., S1237a and S1157b), do not meet our criteria and are, therefore, discarded from further analysis.

Polarized and polarization-filtered P‐ and S-waveform sections rotated into the ZRT system are shown in Figure 4, with only the clearest of the two shown to avoid clutter. For comparison, predicted P‐ and S‐wave travel times based on the most up‐to‐date Mars models from Khan et al. (2023) are also shown (see Data and Resources). The comparison allows several observations to be made. First, phase refinement in the absence of the denoised data would not have been possible. Second, of the 13 events considered, we focus on the 7 events that cover the appropriate epicentral distance range (S−P differential travel time of ∼360–490 s), in which triplicated phases are expected. Third, in spite of careful data selection and processing, waveforms remain noisy. This is not unexpected in view of the fact that most events are either quality B or C. Fourth, with a noncontinuous epicentral distance coverage, observing the characteristic triplicated waveform pattern (compare with Fig. 1) is challenging. Fifth, uncertainties in S−P differential travel time (see Table 1) average about ±5 s, implying that individual events can be considered reasonably well located in distance, but not necessarily in direction (uncertainties associated with back‐azimuth determination of quality B and C events is significant). This implies that the rotation into the ZRT system is possibly only approximate, as a result of which the sought‐after seismic phases are potentially present on multiple components. Finally, because individual arrivals occur within ∼10 and 20 s of the main P and S waves, respectively, and most mantle‐traversing body waves are detected at frequencies between 0.2 and 0.8 Hz, resolution is generally insufficient, particularly in the case of P waves, to distinguish first from secondary arrivals. In the case of event S1102a (quality A), for example, the main P wave, the depth phase, and the triplicated P waves are expected to arrive within just a few seconds of each other, making identification of either of these challenging. Nonetheless, a number of repeating wave packets are present in the transverse‐component waveforms. For the three events that group together toward the lower triplication branch, S1102a–S0421a–S0167b, two recurring wave packets, indicated by arrows in Figure 4, are seen in the waveforms and at about the same relative location. For S1415a–S1153a, toward the middle of the triplication, a single recurring wave packet is observed. On account of the recurrence, the indicated wavepackets are unlikely to be related to random noise and can therefore be employed for model down‐selection. This is illustrated in Figure 4, in which the predicted travel time ranges based on the models of Khan et al. (2023) (gray lines) have been sorted according to overlap with the two wavepackets seen in the seismograms of events S1102a–S0421a–S0167b. The corresponding models are also shown in Figure 4 and indicate MOTD depths in the ranges 987–1052 km (red models matching travel time of first wavepacket) and 1075–1122 km (yellow models matching travel time of second wavepacket) (see inset in Fig. 4), respectively. Events S1415a–S1153a provide little additional information as these appear to be sitting in the center of the triplication, which is compatible with the observation of a single large wavepacket.

Incidentally, our shallower depth range is compatible with the results from Huang et al. (2022), who found a midmantle discontinuity at a depth of 1006 ± 40 km using a waveform matching approach. Although the latter potentially represents a more quantitative means of assessing the depth of the midmantle discontinuity, the study of Huang et al. (2022) nevertheless relied on several assumptions, including known source mechanisms and separate fitting of P and S waveforms using distinct thermochemical models, which showed that a singular model among the structures that fit both types of data could not be isolated. Deng and Levander (2023) interpreted a peak in the noise autocorrelation at a lag of ∼280–285 s to be associated with the olivine transition and, using premission seismic velocity models of Mars, estimated its depth to range between 1110 and 1170 km. Converting their time lag using an average velocity model created from the Khan et al. (2023) model database (see Data and Resources) corresponds to a depth range of 1071 ± 5 km (280 s) and 1092 ± 5 km (285 s), which is compatible with the deeper range found here but incompatible with that of Huang et al. (2022).

Thus, although evidence for a midmantle seismic discontinuity continues to be strengthened, the uncertainties on the location of the Martian MOTD at present are such that any thermochemical inferences made from the mapped MOTD depths must be considered provisional. A joint approach incorporating noise autocorrelation analysis, however, may be a viable way forward to improve upon current estimates of MOTD depth.

We reanalyzed the entire InSight seismic data set for marsquakes and impacts that have interacted with the Martian midmantle seismic discontinuity associated with the olivine transition, which is predicted to occur in the midmantle of Mars. The depth of this discontinuity is desirable because it anchors the thermochemical structure of Mars’s mantle. Our search resulted in seven events at the required epicentral distance from the InSight station. Based on matching differential travel times of triplicated body‐wave phases, we find tentative evidence for a seismic discontinuity in either of the depth ranges 987–1052 km or 1075–1122 km, in line with results obtained from two earlier studies. The uncertainty reflects our current inability to properly separate the different body‐wave phases in the observed waveforms. This is related to the limited number of useful events and their generally poor SNR (most of the available events are quality B and C). With only half of the considered events located at epicentral distances where interaction with a midmantle discontinuity is expected, the poor epicentral distance distribution of the events makes identification of the waveform pattern characteristic of triplications challenging. A combination of different data sets may help resolve this impasse. To seismically image the Martian midmantle structure in more detail in future missions will require a significantly improved SNR, which can be achieved through proper shielding of the seismometer(s) and solid grounding on bedrock. On account of the low seismicity of the planet, with only a single magnitude‐4.6 event recorded in 4 yr, mission longevity will be another important factor.

The InSight event catalog V14 (comprising all events until the end of the mission) and waveform data are available from the Incorporated Research Institutions for Seismology Data Management Center (IRIS‐DMC), National Aeronautics and Space Administration Planetary Data System (NASA PDS), SEIS‐InSight data portal, Institut de Physique du Globe de Paris (IPGP) data center (https://www.insight.ethz.ch/seismicity/catalog/v14, last accessed March 2025) and MarsQuake Service catalog by InSight Marsquake Service (2023). The data were processed with ObsPy, NumPy, SciPy, and TwistPy, and visualizations were created with Matplotlib. Synthetic waveforms were computed employing AxiSEM (Nissen‐Meyer et al., 2014), and travel‐time curves and ray paths employing TauP (Crotwell et al., 1999). The interior Martian structure models from Khan et al. (2023) that are used in the present analysis are available in digital format from doi: 10.18715/IPGP.2023.llxn7e6d. Raw and denoised waveform data for the events considered here are available from doi: 10.5281/zenodo.14228799. The TwistPy data are available at https://twistpy.org/ (last accessed November 2024). The supplemental material includes text and 13 figures, which describe the olivine phase diagram under Martian conditions and data processing, phase refinement and azimuth analysis for all events considered in this study.

The authors acknowledge that there are no conflicts of interest recorded.

The authors acknowledge National Aeronautics and Space Administration (NASA), CNES, partner agencies and institutions (UKSA, SSO, DLR, JPL, IPGP‐CNRS, ETHZ, ICL, MPS‐MPG), and the operators of Jet Propulsion Laboratory (JPL), SISMOC, MSDS, Incorporated Research Institutions for Seismology Data Management Center (IRIS‐DMC) and PDS for providing SEED SEIS data. C. D., A. K., and D. G. would like to acknowledge support from ETH through the ETH+ funding scheme (ETH+02 19‐1: “Planet Mars”). Part of this research was also supported by the Swiss National Science Foundation through Grant 197369. The authors are grateful to three referees for constructive comments that helped improve the clarity of the article. Finally, the authors would like to thank Keith Koper for efficient editorial handling. This is InSight contribution 350.