Seismic detection of a deep mantle discontinuity within Mars by InSight

Observations of Mantle Triplications on Mars.

Earth’s MTZ was first discovered by analyzing the behavior of P-wave travel time curves (35). Early studies on MTZ discontinuities utilized refracted waves bottoming above and below the 410- and 660-km depths, forming travel time triplications, and inverted the data from these triplications for one-dimensional (1D) velocity structures of Earth’s mantle (36, 37). Here we follow a similar method to identify the P and S triplications that resulted from the interactions with a midmantle discontinuity within Mars using the InSight data (Fig. 1 A and B). As an example, the predicted travel time curves from a prelanding model (38) exhibit two intersecting triplications (Fig. 1 C and D): the AB and BC branches are associated with a weak seismic discontinuity from the phase change from orthopyroxene (opx) to high-pressure clinopyroxene (HP-cpx) at ∼800 km depth (39), while the CD, DE, and EF branches are related to the seismic discontinuity from the postolivine transition at ∼1,000 km depth. Here we simply refer to these two discontinuities as 800 and 1,000. On Earth, the opx to HP-cpx transition provides one possible mechanism to explain the X discontinuity at ∼300 km depth, but the lower proportion of opx in Earth’s upper mantle at these depths (<10 vol %) compared to the expected Martian mantle (∼20 vol %) hinders this phase transition from producing a global discontinuity (39). Our study focuses on the analyses of CD, DE, and EF branches of the predicted triplications to constrain the depth and sharpness of the 1,000 discontinuity which is likely to be associated with the postolivine transitions and subsequently infer the thermal and compositional state of the Martian mantle.

The InSight Marsquake Service (MQS) classifies seismic events into low- and high-frequency families based on their dominant frequency content (9, 10). The low-frequency family, including both low-frequency (LF) and broadband (BB) events, is interpreted as a set of marsquakes that originate in the deep crust or mantle and radiate seismic waves deep into the mantle. In contrast, the high-frequency family is attributed to shallow events with most of their energy trapped in the crust (40). Here we focus on the LF and BB events to search for deep mantle triplications. From the MQS catalog, we down-selected five LF and BB events that have $t_S-t_P$ values (differential travel times between P and S waves) consistent with the 60° to 85° distance range expected for mantle triplications: S0185a, S0234c, S0345a, S1094b, and S1102a (SI Appendix, Table S1-1) (41, 42). Among these events, S1094b (∼Mw 4) and S1102a (∼Mw 3) are the largest events and have energy reaching above 5 Hz, while the other three events are dominated by energy between 0.1 and 1 Hz. Due to crustal scattering, it is not straightforward to pick P and S arrivals directly from the unfiltered seismograms, especially for S0234c and S0345a that have lower signal-to-noise ratios (SNRs). Therefore, we used a polarization filter (43), in addition to a Butterworth bandpass filter (0.3 to 0.9 or 0.5 to 0.9 Hz), to improve the SNR of the P- and S-wave signals, while suppressing signals that are not rectilinearly polarized, such as noise and scattered waves. We also carried out frequency-dependent polarization analysis (FDPA) to validate the vertical polarizations for P waves and horizontal polarizations for S waves (44). All of our P- and S-arrival picks must show consistent vertical or horizontal polarizations from the polarization filter and FDPA methods. Finally, to ensure our picks are not contaminated by wind or pressure changes, we performed a comodulation analysis between the wind/pressure data and seismic data (45). Our analyses only kept P- and S-arrival picks that had excess seismic energy above the environmental noise. The resulting dataset of P- and S-arrival picks was then used to determine the epicentral distances and to search for later-arriving triplicated phases.

To facilitate body wave identifications, we aligned the waveforms on their first P and S arrivals and ordered them by $t_S-t_P$ (Fig. 2). This alignment reveals secondary arrivals following the P and S waves whose travel times are consistent with the predicted move-outs of mantle triplications and, furthermore, have amplitudes and polarities similar to the direct P and S arrivals. One possibility is that these arrivals are depth phases (e.g., pP and sS for a range of event depths), but we rule out this interpretation based on the following lines of evidence: 1) The two event depths obtained from interpreting these arrivals as pP and sS phases from a single event are inconsistent with each other for the examined events. 2) The lack of surface waves and short codas of most events imply these LF and BB events occurred in the deep crust or uppermost mantle (10). If source depths are greater than 30 km, the pP and sS phases are then predicted to arrive later than 12 and 22 s after the first P and S arrivals, respectively (11–13); thus, they are unlikely to interfere with the mantle triplications. 3) Depth phases travel twice through the crust at the source side, thereby losing more energy due to crustal scattering and attenuation (4, 40), and would have overall lower amplitudes than the direct phase regardless of the radiation pattern for the event. Furthermore, for surface events (e.g., a meteoroid impact), we would not expect to observe depth phases at all. We also exclude the possibility that these phases arise from crustal conversions or reflections because these phases would only have amplitudes ∼20% of the first P and S waves (47) and would be further suppressed by the polarization filter due to noise contamination (SI Appendix, Fig. S4-11). We only observe crustal S reflections from the two large events, S1094b and S1102a, arriving ∼10 s after the S wave. These phases do not interfere with the candidate S triplications (SI Appendix, Figs. S1-1 and S1-3). With these considerations, we interpret these secondary arrivals as candidate P and S triplications caused by the interactions with the 1,000 discontinuity and investigate their characters further with seismic waveform modeling. We also demonstrate with synthetic modeling that the magnitude thresholds for detecting the P and S triplications are Mw 2.5 (SI Appendix, Fig. S4-9) and 3.5 (SI Appendix, Fig. S4-10), respectively, compatible with the reported magnitude range (Mw 2.6 to 4.0) for LF and BB events (9, 48), indicating that mantle triplications are detectable by InSight.

Synthetic Waveform Modeling.

We used a waveform modeling approach based on a suite of mineral physics-derived radial models for Mars to fit the triplication data and determine the exact depth at which the 1,000 discontinuity occurs (SI Appendix, Fig. S2-2). These 1D models are derived from six different mantle compositions [DW85 (14), EH45 (17), KC08 (46), LF97 (15), TAY13 (19), and YM20 (20)] with mantle potential temperatures ($T_P$) ranging from 1,305 to 1,705 K. We selected a representative subset of 60 models out of a total of ∼2,000 models constructed to cover the entire temperature and composition model space. We used the spectral-element code AxiSEM (49) to compute high-frequency synthetic waveforms ($T_{\mathit{dom}}=2s$) that match the frequency content of the data (0.1 to 1.0 Hz). Since the focal mechanisms and depths of these events are unknown, we assumed a normal faulting source at 30 km depth in the waveform modeling, which is similar to seismic events located at Cerberus Fossae (50). The source depth assumption is supported by a recent body wave study (51) which suggests most high-quality LF and BB events are located at ∼30 km depth. Our synthetic tests performed for different focal mechanisms demonstrate that our waveform modeling results are not significantly affected by the assumption of a normal faulting source (SI Appendix, Fig. S4-2). Given the similar take-off angles at the source, triplicated waveforms do not have a high sensitivity to focal mechanisms, allowing us to model waveforms using a single source by adjusting polarities accordingly. We aligned the events with the synthetics that have similar $t_S-t_P$ values within the uncertainty range (SI Appendix, Table S1-1) using cross-correlations. Correlation coefficients (CCs) between data and synthetics were subsequently calculated to quantify misfits ($\mathit{misfit}=1-CC$).

Depth and Sharpness of the 1,000-km Discontinuity.

We explored a wide range of candidate mantle compositions and potential temperatures to identify models that matched the observed waveforms (SI Appendix, Fig. S1-15). We find two models that show a high correlation match for P (Fig. 2A, $\mathit{mean} CC=0.73$) and S triplications (Fig. 2B, $\mathit{mean} CC=0.77$), respectively. This high consistency between waveform matches supports the identification of these seismic phases as mantle triplications. The waveform modeling was carried out independently for P and S triplications, and each phase showed slightly different sensitivities to the 1,000 depth which were measured at the center depth of the phase transition. The two best-fitting models show consistent $T_P$ values (1,605 K) but different mantle compositions and 1,000 depths: KC08 and 1,013 km for P triplication (Fig. 2A) and YM20 and 1,027 km for S triplication (Fig. 2B). Triplicated waveforms are less sensitive to mantle composition; this topic is discussed in the next section. The ranges of 1,000 depths consistent with the observed P and S triplications are 994 $\pm$ 62 km (Fig. 3 A and C) and 1,008 $\pm$ 49 km (Fig. 3 B and D), respectively. These P- and S-derived 1,000 depths agree within two SDs; however, as both estimates rely on the same interior structure models, we must address covariance between the two estimates. We therefore computed the total misfit to search for best-fitting models in the overlapping depth range that can explain both triplications. We find a subset of 26 models that provide superior fits to the waveforms with misfits lower than 0.33 (Fig. 3E). By averaging these lower misfit models, the 1,000 depth range of the best fitting models is constrained to be 1,006 $\pm$ 40 km, and corresponds to $T_P=1605\pm 100$ K and temperatures between 1,670 and 1,892 K at the 1,000 depth. A recent study based on autocorrelation analysis of ambient noise suggested a mantle discontinuity in the 1,110 to 1,170 km depth range and interpreted this as the olivine–wadsleyite transition (52), which is considerably deeper than our estimate here. However, a detailed study of transient signals (glitches and spikes) arising from deformation internal to the InSight lander and SEIS instrument package showed that the seismic phase interpreted as a mantle discontinuity in that study is likely a signal processing artifact (53).

The 1,000 depth reported above is based upon the assumption of a 30-km-deep source. However, the source depths of these five events are poorly constrained, and there is a direct trade-off between the source depth and 1,000 depth. By assuming a linear areotherm, surface temperature (220 K), and heat flow (20 mW/m²), the maximum seismogenic depth for Mars is calculated to be ∼130 km based on a temperature threshold of 1,073 K (54). We would underestimate the 1,000 depth and $T_P$ by up to 33 km and 100 K, respectively, especially if the source depths of our events fall toward the deeper end of the 30 to 130 km depth estimate. Conversely, we would overestimate the 1,000 depth by up to 10 km if the seismic sources are all located near the surface (e.g., a meteoroid impact). Furthermore, if the events are located between 0 and 30 km depth, the depth phases (e.g., sS and pP) could therefore potentially interfere with the mantle triplications. However, we would expect this interference to have limited effects on the constraints of the 1,000 depth because the amplitudes of the depth phases are likely smaller due to crustal scattering, and their rectilinear particle motions could be obscured by mode conversions and scattering in the crust as well as noise contamination. Therefore, the interference of depth phases could be effectively suppressed by polarization filtering (SI Appendix, Fig. S4-16). Future studies that incorporate detailed inversions of the source depths and focal mechanisms of these five events or more future events could aid in improving the constraints on the 1,000 depth.

Mantle triplications are also sensitive to the sharpness of the 1,000 discontinuity which can be used to constrain the thickness of the postolivine transition. To investigate this, we performed a synthetic waveform sensitivity test by varying the discontinuity thickness between 0 and 200 km in the best-fitting model with lowest total misfit (Fig. 3B, where the 1,000 depth was fixed to 1,027 km). We find that the total misfit is minimized when the discontinuity thickness is between 20 and 100 km (SI Appendix, Fig. S4-5), corresponding to a pressure interval of 0.76 $\pm$ 0.51 GPa. This indicates that the postolivine transition occurs on Mars over a broader depth range compared to the sharp 410-km discontinuity observed on Earth (thickness <10 km or 0.4 GPa) (55). We also investigate the existence of opx to HP-cpx transition by removing the 800 discontinuity in the best-fitting model (Fig. 3B) and then repeating our waveform modeling. We find that the total misfit increased slightly from 0.25 to 0.26 after removing the relatively subtle discontinuity produced by the opx to HP-cpx transition (SI Appendix, Fig. S4-6), indicating the boundary plays a negligible role in the waveform fits. More seismic events from Mars will be required to reach any conclusion on the nature of opx to HP-cpx phase transition.

Polarization Filter.

Seismic events on Mars are often characterized by strong scattered waves as well as wind noise (4), thereby complicating the identifications of body waves. Here we used a polarization filtering technique to enhance rectilinearly polarized waves such as P and S waves, while suppressing other nonrectilinearly polarized waves including background noise or scattered waves. This technique is particularly useful to remove the noise with the same frequency content as the body waves. A polarization filter was previously applied to the detections of core-reflected phases in the Apollo seismic data (73) as well as surface-reflected waves (12, 74) and core phases (13) on Mars. We followed the method by Montalbetti and Kanasewich (43) to design a time-domain polarization filter for improving the SNR of teleseismic P and S waves. First, we computed the covariance matrix of the three-component data (R, T, Z) in a given time window (5 s) around t₀ and the eigenvectors and eigenvalues of this matrix. The rectilinearity at time t₀ is defined as the following equation:

where ${\lambda }_1$ and ${\lambda }_2$ are the largest and second largest eigenvalues, respectively, and n and J are the empirical exponents. The eigenvector of the principal axis is $\overrightarrow{E}$ = (e₁, e₂, e₃) with respect to the R, T, Z coordinate system. The direction function on each component at time t₀ is given by

where i = 1, 2, 3 (R, T, Z) and K is the empirical exponent. We chose these empirical values for the exponents: n = 0.5, J = 1, and K = 2 (43). The polarization filter for each component is defined as the product of rectilinearity and the corresponding direction function. The filtered three-component seismograms at any time t are given by

FDPA.

We applied FDPA on waveforms of our five LF/BB events used in the analysis. This method has been previously implemented to detect and identify the core-reflected ScS phases on Mars (13). To summarize, we computed the S transform of three component waveforms recorded by InSight SEIS Very Broadband (VBB) and obtained a 3 × 3 cross-spectral covariance matrix in 90% overlapping time windows whose duration varies inversely with frequency. The relative sizes of the eigenvalues of this covariance matrix are related to the degree of polarization of the particle motion, while the complex-valued components of the eigenvectors describe the particle motion ellipsoid in each time–frequency window (44). We focused on when the particle motion of our seismic signal was rectilinear so we could effectively describe the orientation of the particle motion with the inclination and azimuth. To pick the direct P and S arrivals in our seismic records, we used a pair of metrics that can emphasize horizontally (HRM) and vertically rectilinear motion (VRM) and searched for those arrivals whose average HRM or VRM across different frequencies were high (see more details in the supplement of ref. 13).

Comodulation Analysis.

In this section we adopt the comodulation analysis described in ref. 45 to evaluate the contributions of environmental signals into the seismic records and thus aid discrimination between weather-induced noise and seismic energy for each of the candidate mantle triplication events. Comodulation evaluates the correlation in signal power between the measured environmental variables (wind and pressure) and ground motion. It has been demonstrated particularly successful in identifying seismic energy that is in excess of the expected noise injected from the local weather (9, 45, 75).

The comodulation approach uses the method of moments to match the mean and variance of the seismic and environmental signals prior to the start time of the putative marsquake. We first estimate these first two moments for all signals for a total of 15 min prior to the marsquake event and remove the effect of any identified glitches in the record which contaminate the statistical moments. The method of moments is then used to match the mean and variance of the seismic acceleration and pressure rms envelopes to the measured wind speed signal. The seismic rms envelope is estimated as the total power of all three combined Z, N, and E components. Both seismic and pressure rms envelopes are extracted in discrete half-octave frequency bands with center frequencies from 1/26.5 up to 6.8 Hz, with the last frequency bin extending up to the anti-alias filter (8 Hz) of the continuous 20 samples per second VBB data.

In this approach, comodulation shows that the envelopes of the seismic power continuously track the pressure envelopes and wind speed, and therefore, any divergence from the expected noise match can then be interpreted as being from a seismic source. Each of the candidate events in the manuscript for mantle triplications located at epicentral distances between 60° and 85° is henceforth evaluated with respect to the comodulation analysis. For events that do not have weather data due to lander power constraints after sol 700, highly sensitive lander resonances excited by the wind and atmospheric phenomena are used as proxies for atmospheric noise injection (45, 76).

Data Processing Steps and Phase Identification.

Marsquakes are broadly classified as two families: 1) low-frequency (LF) and 2) high-frequency (HF) events based on their frequency content (9, 10). The dominant energy of LF events is below 1 Hz, while some events have more broadband energy extended above 1 Hz, which are classified as a subgroup: broadband (BB) events. The HF events have energy mostly above 1 Hz and can extend up to 6 to 12 Hz (40). MQS ranks these events from quality A to D (A being the highest quality) based on the clarity and polarization of the phases (9). In this study, we focused on the LF and BB events that travel deep into the mantle to search for triplications from the MTZ. We have identified five LF and BB events with $t_S-t_P$ values that cause them to fall in the target distance range for triplications: S1094b (SI Appendix, Fig. S1-1), S1102a (SI Appendix, Fig. S1-3), S0185a (SI Appendix, Fig. S1-5), S0234c (SI Appendix, Fig. S1-7), and S0345a (SI Appendix, Fig. S1-9).

We first removed the instrument response from the data to transfer seismograms to velocity, and the three SEIS VBB components were rotated from UVW configuration to the ZNE coordinate frame (8). Then, we band pass filtered the seismogram between 0.3 and 0.9 or 0.5 and 0.9 Hz to remove long-period noise and high-frequency resonances and eliminate the 1-Hz tick noise resulting from cross-talk between the temperature sensors and SEIS VBB seismometers (77). The choices of band pass filters are informed by the comodulation analysis and filter banks: we used the 0.5 to 0.9 Hz filter for the events (e.g., S1102a) that do not have clear excess seismic energy below 0.5 Hz; otherwise, we used the 0.3 to 0.9 Hz filter. We first applied polarization filtering in the ZNE coordinate to pick the P arrival of each event on the vertical component (BHZ) based upon the MQS reported start time of seismic energy (9). We then compared our picks to the FDPA results to examine the vertical polarizations for P waves (e.g., SI Appendix, Fig. S1-1A). We only picked the arrivals when these two methods show consistent peak energy, then assigned the uncertainties to include both peaks. To obtain event back-azimuth, we computed the average principal axis of P-wave particle motions, which is expressed as $\overrightarrow{P}=(P_E,P_N,P_Z)$, using a 5-s time window centered around the P wave. The back-azimuth of the event is then calculated using the following equation (e.g., ref. 78):

We used the P-wave–derived back-azimuth to rotate the horizontal components (BHN and BHE) to radial (BHR) and transverse (BHT) components. The event back-azimuths are summarized in SI Appendix, Table S1-1. However, due to the low amplitudes of P waves, the back-azimuths are often not well constrained, so we examined both the BHR and BHT components for our S-wave triplication study.

After rotation, we applied the polarization filter, now in the ZTR coordinate system, to pick the S arrivals with the guidance of MQS reported picks. Ideally, we would prefer to pick S arrival on the transverse component to avoid P to S converted phases. However, since the back-azimuths are poorly constrained, we picked S arrivals on both radial and transverse components and only kept the pick with the larger SNR. For example, the S wave of S1094b event is two times larger on the radial component than the transverse component (SI Appendix, Fig. S1-1B). Therefore, we picked this S arrival on the radial component. Similarly to P wave, we compared our S-arrival picks to the horizontal polarizations from FDPA to ensure these two methods show consistent peak energy within the uncertainty range (e.g., SI Appendix, Fig. S1-1B). All the P- and S-arrival picks were made at the maximum amplitudes of the polarization filtered envelopes. Finally, we examined our P- and S-arrival picks using the comodulation analysis: both picks must show excess energy above the environmental noise (e.g., SI Appendix, Fig S1-2); otherwise, they were discarded from our analysis. For example, we identified two other events (S0167a and S0167b) that are located in the distance range appropriate for mantle triplications, but their P waves do not show any excess energy (SI Appendix, Figs. S1-11 and S1-12). Therefore, we excluded these two events from our analysis due to the wind contamination. We repeated the same picking processes for the other four candidate events (SI Appendix, Figs. S1-3–S1-10). We calculated the SNR of P and S waves as the ratio between the maximum envelope and the average envelope amplitude in a noise time window (−150 to −50 s before P arrival). P waves typically have smaller SNRs between 2 and 7, whereas S waves show larger SNR between 3 and 51 (SI Appendix, Table S1-1).

SEIS is deployed on the surface of Mars where a harsh environment inevitably causes artifacts in the data. Although SEIS is protected under the wind and thermal shield, the daily internal temperature variations can still reach up to 15 K. The thermal variations lead to one type of instrument self-noise manifested as large amplitude deviations over a short time period, hereafter referred to as a “glitch”. A typical SEIS glitch appears as a large one-sided pulse accompanied by high-frequency spikes near the onset (79). We highlight some large glitches in SI Appendix, Figs. S1-1–S1-10, and they often appear as very broadband energy on the spectrogram. We avoided these known glitches identified by MQS during our phase picking by excluding any picks that fell within ±20s of a glitch signal (e.g., SI Appendix, Fig. S1-5). As informed by comodulation analysis, we also excluded the time windows that were contaminated by wind noise or pressure change (e.g., SI Appendix, Fig. S1-1).

Synthetic Waveform Modeling.

We derived a suite of seismic structure models based on six composition models [DW85 (14), EH45 (17), KC08 (46), LF97 (15), TAY13 (19), and YM20 (20)] with mantle potential temperature $T_P$ ranging from 1,300 to 1,700 K (see Structure Models of Mars for more details). We propagated synthetic seismic waves through these seismic structure models to simulate the waveforms of triplications for comparison with the data. We used the spectral element code AxiSEM (49) to generate synthetics with a dominant period of 2 s. The focal mechanisms and depths of these five triplication events are unknown, but the three Cerberus Fossae events (S0173a, S0235b, and S0183a) were constrained to have normal faulting sources between 33 and 40 km depth (50). Therefore, we assumed a similar tectonic setting for these five LF and BB events and selected a normal faulting source at 30 km for the waveform modeling. Note that a trade-off exists between the source depth and 1,000 depth (e.g., increasing the source depth by X km would cause the 1,000 depth to be ∼X/3 km deeper). We chose a moment tensor to allow for seismic energy in the receiver direction from all three types of waves: P wave, SH wave, and SV wave (SI Appendix, Fig. S1-13). The receivers were placed on the equator from 55° to 85° with a 0.1° spacing (5.9 km). We artificially increased the crustal attenuation ($Q_{\mu }=Q_{\kappa }=10$) to suppress the crustal conversions and reverberations after P and S arrivals, while the mantle attenuation was muted ($Q_{\mu }=Q_{\kappa }={10}^5$) to keep the high-frequency content of the triplicated waveforms. We suppressed the crustal phases as they likely vary between our observed events and are unconstrained in the models. In addition, we demonstrated that the polarization filtering can suppress the crustal phases in the data due to their low amplitudes, especially for events with lower SNR like S0234c and S0345a (SI Appendix, sections 4.5 and 4.6). We used these tests to ensure that the waveform modeling fits are not biased by crustal phases in both the data and the synthetics.

To show the characteristics of triplicated waveforms, we present six example models from our model space that span the range of possible 1,000 depths (938 to 1,066 km). We aligned the synthetics on the first P and S arrivals and ordered them by $t_S-t_P$ in the same way as we did with the observed seismic data (SI Appendix, Fig. S1-14). Triplicated waveforms display a K shape of travel time move-out after the alignment on P or S arrivals. The $t_S-t_P$ values at the center of the triplications generally increase with the 1,000 depth (or $T_P$). We observe crustal reverberations after the P wave on the vertical component (SI Appendix, Fig. S1-14A) and after the S wave on the radial component (SI Appendix, Fig. S1-14C), whereas the transverse component mostly contains triplications (SI Appendix, Fig. S1-14B). S to P conversions arrive ∼5 s before the first S arrival on the radial component (SI Appendix, Fig. S1-14C). These crustal phases also exhibit triplications but have much smaller amplitudes than the first P and S arrivals, thereby only contributing minimal misfits to the waveform modeling. Depth phases pP and sS are expected to arrive ∼12 and 22 s after the P and S waves, respectively, but they are eliminated in the synthetics as a result of the strong attenuation in the crustal model.

After computing the synthetics, we used them to interpret the triplicated waveforms. Since the epicentral distances of these events are unknown, we matched the synthetics and data based on $t_S-t_P$ values. For each event, we searched through the synthetics within the $\pm 5$ s uncertainty range of $t_S-t_P$ measurement (SI Appendix, Table S1-1) to find the best-fitting synthetic waveform. As the uncertainties of $t_S-t_P$ vary from the events, we also tested $\pm$7.5 and $\pm$10 s uncertainties for the waveform modeling. We obtained similar results despite using different uncertainties, thereby retaining the $\pm$5 s uncertainty range for the waveform modeling. We only investigated the vertical component for P triplication. Due to the less well-constrained back-azimuth, there is uncertainty in the horizontal rotation of the events. Therefore, we used both transverse and radial components for S triplication but only kept the one with a lower misfit. As the focal mechanisms of these events are unconstrained, we tested both positive and negative polarities for each synthetic waveform and again kept the one with a lower misfit. During the waveform modeling, we applied the following data processing steps for each synthetic waveform. 1) We normalized the amplitudes of the data and synthetics to unity. 2) We cut a time window around the first P and S arrivals: −5 to 8 s for P wave (SI Appendix, Fig. S1-15A) and −5 to 20 s for S wave (SI Appendix, Fig S1-15B). S0185a and S1094b events required extended time windows to account for possible later arrivals of S triplications. We used −5 to 30 s for the S waves of S0185a and S1094b events (SI Appendix, Fig. S1-15B). 3) We aligned the cut synthetics with the data using a cross-correlation method. 4) Finally, we computed the normalized CCs between the data $d(x,t)$ and synthetic waveforms $u(x,t)$ using the following equation (80):

We only kept the synthetics with the largest CC for each event in a given model. We then computed the mean value of CC in each model as follows:

where N represents the number of events ($N=5$). We define the misfit function for each model using the following equation:

We computed the weighted sum of the misfits of P and S triplications as the total misfits:

where ${\chi }_P$ and ${\chi }_S$ are the misfits of P and S triplications, respectively. We assigned a weighting factor of two to S triplications because the time window length for S triplication is about two times longer than P triplications, and the SNR of S waves are generally larger.

Here we show six example waveform modeling results from our model space (SI Appendix, Fig. S1-15). The synthetic waveforms are plotted on top of the data for direct comparisons. Misfits of S triplications (SI Appendix, Fig. S1-15A) show larger variations as a function of the 1,000 depth compared to P triplication (SI Appendix, Fig. S1-15B). Both misfits are minimized in the intermediate 1,000 depth range (1,000 to 1,030 km). In our model space, we found 26 models that can simultaneously fit both P and S triplications, including our best-fitting model ($T_P=1,605 K$, YM20).