Seismic Autocorrelation Analysis of Deep Mars

The InSight mission deployed one seismic station on Mars, providing a chance to apply single‐station‐based autocorrelation analysis to investigate Martian subsurface structures. However, recent analysis indicated the low‐frequency autocorrelation signals may originate from quasi‐periodic high‐amplitude instrumental “glitches” rather than the reflection response of deep Mars. In this study, we detected and removed these high‐amplitude glitches in raw seismic data and employed autocorrelation on the clean vertical component waveforms filtered between 0.05 and 0.1 Hz. We observed signals at the expected times for the olivine‐wadsleyite transition and core‐mantle boundary (CMB) as estimated by other methods. This result suggests that the low‐frequency autocorrelation signals are the reflection response from the olivine‐wadsleyite transition in the mantle and the Martian CMB region, rather than a noise phenomena. A grid search method to fit the observed PcP waveform was used to identify a layer intermediate in velocity between the Martian mantle and core at the Martian CMB.


10.1029/2023GL105046
2 of 9 recordings (e.g., ambient noise, coda waves) has been widely employed to extract body-wave reflection phases to reveal the subsurface structures of the Earth (Delph et al., 2019;Gorbatov et al., 2013;Kennett, 2015;Oren & Nowack, 2017;Phạm & Tkalčić, 2018, 2021;Qin et al., 2020;She et al., 2022;Tibuleac & von Seggern, 2012;Zhou & Zhang, 2021) and Moon (Nishitsuji et al., 2016(Nishitsuji et al., , 2020)).Recent autocorrelation analysis on Mars detected a strong P-wave reflection phase at 10.5-11.5 s two-way traveltime (Compaire et al., 2021;Deng & Levander, 2020;Schimmel et al., 2021).This 10.5-11.5 s signal can be interpreted as the reflection from either a mid-crust discontinuity or the Martian Moho at approximately 35 km (Deng & Levander, 2020).To examine deeper seismic structure of Mars, Deng and Levander (2020), hereafter referred to as DL2020, filtered the stacked autocorrelations in the low-frequency 0.05 and 0.1 Hz band and interpreted signals at ∼280 and ∼375 s as P-wave reflections from the Martian olivine-wadsleyite transition and CMB, respectively.Depth conversion using different reference velocity models determined the Martian core radius in the range of 1,790-1,870 km (Deng & Levander, 2020), consistent with the core radius estimation from other studies (Khan et al., 2018;Rivoldini et al., 2011;Yoder et al., 2003).However, Lognonné et al. (2020) identified high-amplitude seismic glitches within raw SEIS data which probably originate from stress relaxation within the InSight seismometer.Scholz et al. (2020) subsequently analyzed the glitch properties.These glitches may contaminate the final results of seismic analysis if not removed (Compaire et al., 2021).Barkaoui et al. (2021) analyzed the quasi-periodic recurrence time of these seismic glitches and suggest that it coincides with the arrival times of identified in low-frequency autocorrelations in DL2020.Kim et al. (2021) also did autocorrelation tests that the signals interpreted as olivine-wadsleyite transition and CMB in DL2020 may be contaminated by the quasi-periodic high-amplitude glitches.As more seismic data on Mars have been collected since DL2020, further analysis of deep Martian structure is warranted.
In this work, we designed a test to resolve whether the low-frequency autocorrelation signals in DL2020 are real seismic reflection signals or result from instrumental noise.The high-amplitude glitches in raw SEIS data were first detected to generate three data sets for further analysis: (a) the raw continuous data with glitches; (b) the data with only detected glitches; (c) the clean continuous data after the application of glitch removal algorithms, hereafter referred to as the "cleaned" or "deglitched" data set (Scholz et al., 2020).We autocorrelated the ambient noise recorded on the vertical-components of these three data sets, raw, glitch only, and cleaned, to demonstrate that the low-frequency autocorrelation signals reported by DL2020 are the reflection response from deep Mars rather than the recurrence time of quasi-periodic glitches.Using the velocity model of Stähler et al. (2021) we calculated synthetic PcP seismograms to model the complexity of the PcP phase observed in the cleaned autocorrelation.A layer of intermediate velocity between the Martian mantle and core at the Martian CMB was identified using a grid search to minimize the misfit between the synthetic and observed PcP phases.

Data and Methods
The SEIS data were collected with a 10 and 20 Hz sampling rate in two different periods.Prior to preprocessing, the open-source package SEISglitch (Scholz et al., 2020) was applied to extract and remove the high-amplitude glitches from the raw continuous SEIS waveforms to produce cleaned and glitch-only data.Figure 1a-1c provide a 1-day comparison of original and cleaned continuous U-, V-and W-component data.The glitches were mostly detected at night time during periods of relatively low ground temperature (Figure 1d) and windspeed (Figure 1e) because of the lower wind-and pressure-driven noise level during the nighttime (Scholz et al., 2020).Figures 1g  and 1h show the spectrograms of original and cleaned W-component waveforms in a 3-hr window (Figure 1f).The glitches are recorded as high-amplitude vertical peaks in spectrogram (Figure 1g).The glitch energy is effectively reduced by the glitch removal algorithm within SEISglitch package (Figure 1h).Figures S1 and S2 in Supporting Information S1 describe parameter selection test for glitch detection and removal, which suggest that parameter values used in SEISglitch have modest impact on the final autocorrelation results.The preprocessing procedures included the instrument response removal from the raw, clean and glitch-only continuous SEIS data with a broad, 0.01-3.5 Hz, bandpass filter.Data sampled at 20 Hz were decimated to 10 Hz.The continuous waveforms were cut into 4-hr-long segments and rotated from U-V-W into orthogonal Z-N-E channels (Compaire et al., 2021).We removed the mean and trend and applied a taper to each 4-hr window.
We used a workflow similar to Bensen et al. (2007) to calculate ACFs for the 3 data sets.Temporal balance and spectral whitening were employed to improve the interpretability of the stacked ambient noise ACFs (Bensen et al., 2007;Oren & Nowack, 2017).Each individual 4-hr trace was then filtered from 0.05 to 0.1 Hz.We computed and normalized the vertical component autocorrelograms by the amplitude at zero lag time for each 4-hr trace.We applied both linear (Linearly stacked (LS)) and phase-weighted (PWS) stacking (Schimmel & Paulssen, 1997) to produce stacked autocorrelograms.The power used for PWS, which defines the significance of coherency measure (Schimmel & Paulssen, 1997), was empirically selected as 2 (Korenaga, 2014;Niu & Chen, 2008;Wookey & Helffrich, 2008).

Autocorrelation Results
Figure 2 compares the LS and PWS stacked autocorrelation of raw, cleaned, glitch-only and 80% glitch-only data sets.The 80% glitch-only data set consists of ∼80% of largest glitches within glitch-only data set.We made this

10.1029/2023GL105046
4 of 9 comparison between glitch-only and 80% glitch to ascertain the importance of small glitches left by the glitch removal algorithm (Scholz et al., 2020).The stacked ACFs from these two data sets are almost identical for both LS (Figure 2a) and PWS (Figure 2b), which suggests the smaller amplitude glitches remaining in the clean data set will have modest impact on final ACFs.
Two signals at ∼285 and ∼385 s are clearly observed in the stacked ACF.The signal at ∼285 s in the clean ACF (red curves in Figures 2a and 2b) is consistent with the traveltime of the P-wave reflection from the olivine-wadsleyite transition as calculated for a number of Martian velocity models (Khan et al., 2023;Stähler et al., 2021).
We turn our attention to the analysis of the interpreted PcP phase at ∼385 s, leaving further discussion of the ∼285 s signal to the Supporting Information S1 (Text S4 and Figure S5).The arrival time of this signal is consistent with the prominent phase identified as the P-wave reflection from CMB reported in DL2020.The bootstrap calculation also indicates the observation of this PcP signal is robust (Figure S4 in Supporting Information S1).The ACFs of glitch-only and 80% glitch-only data do not show a prominent response at ∼385 s.We take this as indicating the reflection response at ∼385 s originated from Martian CMB rather than an unfortunate correlation of glitches suggested by Barkaoui et al. (2021) and Kim et al. (2021).The noise sources to generate this deep Mars reflection signal may come from atmospheric phenomena (Nishikawa et al., 2019) or coda waves of Marsquakes (Lognonné, Schimmel, et al., 2023), which are less attenuated during wave propagation as Mars is a dry, cold and tectonically inactive planet compared with Earth (Menina et al., 2021).In the stacked raw data ACF there is a signal at ∼385 s, but the reflection pulse is much more complicated than in the clean waveform because it is contaminated by the high-amplitude seismic glitches in the raw data.The PWS stacked autocorrelation of the deglitched ambient noise data (red curve in Figure 2b) will be used for the comparison with synthetic results since the observed reflection phases are clearer and simpler than the linear stack.

Comparison With Synthetic PcP Phases
We simulated the synthetic normal-incident PcP phases with Thomson-Haskell matrix method (Haskell, 1962;Thomson, 1950) for a set of velocity models (Figure 3a) from a probabilistic inversion of seismic traveltimes, Love number k 2 and moment of inertia (Stähler et al., 2021).The comparison between synthetic and real PcP signals (Figure 3b) shows that the velocity model 3 standard deviations lower than the mean model of Stähler et al. (2021) best matches the timing of the observed PcP phase.This observation is consistent with the Mars orbiting surface wave results retrieved from ambient noise autocorrelation (Deng & Levander, 2022).We will use Stähler et al. ( 2021)'s mean-3σ model (Figure 3a) as a reference in later analysis.

Grid Search for CMB Transition Zone Modeling
The PcP reflection waveform at ∼385 s is more complicated than that from the potential olivine-wadsleyite transition at ∼ 285 s (red curve in Figure 2b), therefore we model the CMB as a transition zone (TZ), with the velocity V TZ = αV mantle + (1 − α)V core , where 0 ≤ α ≤ 1.For α = 1, this is the mean-3σ velocity model shown in Figure 3a which best fits the observed PcP phase among the Stähler et al. (2021) suite of models (Figure 3b).Two parameters, the TZ thickness and α value, were determined by grid search to find the best combination that fits the observed PcP waveform shown in Figure 2b.The TZ thickness is set as 40, 60, 80, 100, and 120 km and the α value ranges from 0.05 to 1 with 0.05 grid size.For each TZ thickness and α value combination, we simulated the synthetic zero-offset PcP phase by Thomson-Haskell matrix (Haskell, 1962;Thomson, 1950) and calculated the correlation coefficients between the synthetic and observed PcP phases.Figures S6 to S10 in Supporting Information S1 showed the correlation coefficients and waveform comparisons with the observed PcP phase for 40, 60, 80, 100, and 120 km thick CMB transition zone models, respectively.Figure 4a illustrates the 40, 60, 80, 100, and 120 km thick CMB transition zone models that best match with the observed PcP phase.Figure 4b makes the comparisons between the observed PcP and the synthetics of the velocity models shown in Figure 4b.We can see the 60 km thick transition zone model (green model in Figure 4a) at the Martian CMB with V TZ = 8.05 km/s best fits the observed PcP phase (Figure 4b).Since a 60 km thick 1-layer transition zone produces a waveform correlation coefficient of 0.984 (Figure 4b), hence increasing the model complexity of CMB transition zone is unnecessary.

Possible Cause of CMB Transition Zone
A similar ultralow-velocity zone (ULVZ) near CMB has been observed on Earth across different regions, such as Southeast Asia (Yao & Wen, 2014), Hawaii (Cottaar & Romanowicz, 2012), Mexico (Havens & Revenaugh, 2001), Africa and Eastern Atlantic (Helmberger et al., 2000).ULVZs normally are several tens of kilometers thick and may yield P-wave and S-wave velocity reduction by up to ∼10% and ∼30%, respectively (McNamara, 2019).In general, two theories are widely acknowledged to interpret the existence of ULVZ at Earth's lower mantle: (a) partial melt occurs at the location of ULVZ other than the ambient mantle (Ohtani & Maeda, 2001;Wen & Helmberger, 1998;Williams & Garnero, 1996); (b) core rust is formed when iron interacts with water or hydroxide minerals at the lower mantle (Hu et al., 2016;Liu et al., 2017;Mao et al., 2017).The low-velocity layer near the Martian CMB observed in our study (Figure 4) may follow the same theories as the ULVZ observed on Earth.
Two recent studies have presented several different Martian velocity models with a basal melt layer consisting of a low velocity silicate melt, 100-175 km in thickness and 4.5-5.5 km/s in P-wave velocity (Khan et al., 2023;Samuel et al., 2023), where the transition zone thickness is thicker than our estimation of 60 km and transition zone P-wave velocity is lower than our estimation of 8.05 km/s.To better constrain the structure and understand the cause of the CMB transition zone on Mars, advanced studies on high-pressure experiment, mineralogical simulation, and thermodynamic modeling are required to cross-validate with the seismic observations.Kim et al. (2021) and Barkaoui et al. (2021) suggested that the autocorrelation signals identified as the reflection responses from olivine-wadsleyite transition and CMB (Deng & Levander, 2020) may originate from the quasi-periodic glitches rather than real seismic discontinuities.In this study, we applied the autocorrelation method on raw, clean, and glitch-only vertical-component data to investigate these claims.The prominent signal at ∼385 s was extracted from the autocorrelation of raw and clean waveforms but not from glitch-only data (Figure 2).This suggests the signal at ∼385 s was the reflection response from seismic boundaries within Mars rather than the recurrence time of high-amplitude glitches in raw SEIS data.We interpret the ∼385 s signal as the CMB.We interpret an earlier phase at ∼285 s (Figure 2) as the P-wave reflection at olivine-wadsleyite transition as in DL2020 (Text S4 in Supporting Information S1).The noise sources to produce deep Mars reflection signals may originate from atmospheric phenomena (Nishikawa et al., 2019) and/or Marsquake events (Lognonné, Schimmel, et al., 2023).The clarity of these signals likely results from the low seismic attenuation of the Martian interior (Menina et al., 2021).Timing of the PcP phase (∼385 s) is best matched using Stähler et al. (2021)'s mean-3σ model (Figure 3b), consistent with Mars orbiting surface wave R2 results (Deng & Levander, 2022).We then calculated synthetic PcP phases for a velocity transition zone at the CMB in which we varied velocity and layer thickness to model the observed complexity of the PcP waveform (Figure 4).The grid search shows 10.1029/2023GL105046 7 of 9 that a 60 km thick layer with V TZ = 8.05 km/s (green model in Figure 4a) produces synthetic seismograms that best match the observed PcP waveform (Figure 4b).Samuel et al. (2023) and Khan et al. (2023) derived a similar structure using InSight seismic data, which is interpreted as the molten silicate layer at the base of Martian mantle.The observed CMB transition zone on Mars (Figure 4a) can potentially enhance our understanding of geological history and dynamic evolution of Mars (Zhang et al., 2023).This study was funded by the Department of Earth, Environmental and Planetary Sciences at Rice University.We thank Doyeon Kim, Paul Davis, Mark Panning, John-Robert Scholz, Martin Schimmel and other members of the InSight mission for sharing their understandings about the autocorrelation results and the instrumental noise of SEIS data.We thank Simon Stähler for providing the velocity models in Stähler et al. (2021).We also thank Fenglin Niu at Rice University for discussions related to this work.Constructive comments from Editor Kevin Lewis, Reviewer Philippe Lognonné and an anonymous reviewer are appreciated.

Figure 1 .
Figure 1.(a) Original (blue) and Cleaned (orange) U-component waveforms on 3 January 2020.(b) Same as (a) but for V-component.(c) Same as (a) but for W-component.(d) Corresponding ground temperature data and (e) windspeed data on.(f) Original (blue) and Cleaned (orange) W-component waveforms in a 3-hr window (3 January 2020 03:00:00 UTC to 3 January 2020 06:00:00 UTC) shown as the red box in (c).(g) The spectrogram of original W-component waveform in (f).The black arrows identify glitch examples, which are characterized as broadband (i.e., vertical) features within the spectrogram.(h) The spectrogram of cleaned W-component waveform in (f), where the reduction of glitch energy is clear in comparison to the spectrogram in (g).

Figure 3 .
Figure 3. (a) Mean velocity model and 3 standard deviations above and below the mean model in Stähler et al. (2021).The solid lines are P-wave velocity profiles, and the dashed lines are S-wave velocity profiles.(b) Comparison between the synthetic PcP phases for the velocity models shown in (a) and the observed PcP phase from the autocorrelation of the clean ambient noise data (red curve in Figure 2b).The lag time represents the correlation time shift between the synthetic and observed PcP phases, where the positive values represent the observed PcP travels slower than synthetic PcP, and vice versa.

Figure 4 .
Figure 4. (a) The models for 40, 60, 80, 100 and 120 km thick transition zones that best match with the observed PcP phase (red curves in Figure 2b).(b) The comparison between the synthetic PcP phase of the velocity models shown in (a) and the observed PcP phase.The correlation coefficients between the synthetic and observed PcP are shown on the right side.The blue box marks the model with highest correlation coefficient among all models shown in (a) (The 60 km thick 1-layer core-mantle boundary transition zone model, with V TZ = 8.05 km/s).