Seismic Noise Autocorrelations on Mars

Mars is the first extraterrestrial planet with seismometers (Seismic Experiment for Interior Structure, SEIS) deployed directly on its surface in the framework of the Interior Exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) mission. The lack of strong Marsquakes, however, strengthens the need of seismic noise studies to additionally constrain the Martian structure. Seismic noise autocorrelations of single‐station recordings permit the determination of the zero‐offset reflection response underneath SEIS. We present a new autocorrelation study which employs state‐of‐the‐art approaches to determine a robust reflection response by avoiding bias from aseismic signals which are recorded together with seismic waves due to unfavorable deployment and environmental conditions. Data selection and segmentation is performed in a data‐adaptive manner which takes the data root‐mean‐square amplitude variability into account. We further use the amplitude‐unbiased phase cross‐correlation and work in the 1.2–8.9 Hz frequency band. The main target are crustal scale reflections, their robustness and convergence. The strongest signal appears at 10.6 s, and, if interpreted as a P‐wave reflection, would correspond to a discontinuity at about 21 km depth. This signal is a likely candidate for a reflection from the base of the Martian crust due to its strength, polarity, and stability. Additionally we identify, among the stable signals, a signal at about 6.15 s that can be interpreted as the P‐wave reflection from the mid‐crust at about 9.5 km depth.

SEIS has been deployed on the ground of the Homestead hollow, 1.81 m south of the closest lander foot . Homestead hollow is a degraded impact crater of about 27 m diameter in Elysium Planitia with a smooth sandy, granule-and pebble-rich surface with few rocks (Golombek et al., 2020). After leveling and functionality checks SEIS has been covered with a wind and thermal shield (WTS) to further protect the sensors from environmental factors such as laminar and turbulent winds which can reach about 20 m/s, low atmospheric temperature down to below −100°C, daily temperature variability in the order of 80°C and atmospheric pressure changes of a few Pa . The planetary boundary layer height is several kilometers in daytime conditions, in the nighttime its height collapses to a couple of 100 m, with turbulence switching from being thermally driven to being shear-driven under highly stable conditions owing to surface radiative cooling.
In spite of these harsh conditions, 174 mostly small marsquakes  have been identified until 30 September 2019 thanks to the sensitivity of SEIS and its careful deployment on the Martian surface. Nevertheless, due to the unfavorable deployment conditions different aseismic signals are being recorded together with events of the seismic origin Lognonne et al., 2020;Scholz et al., 2020;Stutzmann et al., 2021). That is, the seismic recordings contain a wealth of features caused by operational activities at the lander, noise induced by the lander due to atmospheric perturbations as laminar and turbulent winds, and artifacts caused by the response of the instruments to the variability of the severe climatic factors as pressure and temperature. Data processing requires therefore special care and adapted approaches to avoid any bias in the results and even misinterpretation of signals.
So far, no unambiguous surface waves have been detected in the SEIS data (Clinton et al., 2020;Giardini et al., 2020;Stutzmann et al., 2021). The very shallow structure down to 20 m depth at the InSight landing site has been probed using the travel time of multiple hammer strokes of the HP 3 (Heat Flow and Physical Properties Package) instrument and compliance observations (Kenda et al., 2020;Lognonne et al., 2020). The compliance approach employs seismic measurements of the surface deformation in response to pressure loading by dust-devils. First average crustal seismic attenuation and scattering results have been obtained from the body-wave coda shape of the best recorded marsquakes while first crustal layering has been inferred identifying body wave conversions in receiver functions . These results point to a crustal attenuation which is 3 times larger than on Moon and a highly variable upper crust layer of about 8-11 km thick with S-wave velocities of 1.7-2.1 km/s. The exact thickness of the crust is still the object of investigation, but receiver functions, event coda, and noise autocorrelation results suggest that the crust is either about 15-25 or 31-47 km thick (Knapmeyer-Endrun et al., 2021). Deng and Levander (2020) have also computed noise autocorrelations for Mars and observe persistent signals which they interpret as Moho, upper-mantle transition, and core-mantle boundary reflections. Their mantle signals have not yet been reproduced by other groups and can possibly be due to repeated glitches in the SEIS-data or interference of lander resonances (Kim et al., in preparation). High-frequency autocorrelations (5-7 Hz) have been presented by Suemoto et al. (2020) who showed robust signals which indicate the presence of two shallow reflectors in the first hundreds of meters. Noise and Marsquake autocorrelations have been also presented by Compaire et al. (2021), who find stable autocorrelations with signals related to the crustal structure only during low-noise periods around 2.4 Hz.
In this study, we use the VBB SEIS data to further analyze noise autocorrelations through independent approaches. Special care is given to avoid imprints of aseismic signals and artifacts on the autocorrelation results. The main objective is to extract a stable and clean reflection response from the noise recordings. We focus on crustal scale signals from vertical component recordings within the 1-9 Hz frequency band, and use differently selected subsidiary data sets. The data selection strategy is novel and we show that a stable noise response is built up within a few Martian days of the data and that it consists of several signals. Some of the signals are stable over a broad frequency band which makes them possible candidates for reflections. In any case, the autocorrelations are still difficult to interpret and possible ambiguities in the signal interpretations are discussed.

Determination of P-Wave Reflection Response at SEIS
On Earth, seismic ambient noise correlations have been successfully used to map crustal scale discontinuities (e.g., Becker & Knapmeyer-Endrun, 2018;Buffoni et al., 2019;Gorbatov et al., 2013;Kennett et al., 2015;Oren & Nowack, 2016;Romero & Schimmel, 2018;Ruigrok et al., 2011;Saygin et al., 2017;Taylor et al., 2016;Tibuleac & von Seggern, 2012). Strictly, the cross-correlation of a diffuse wavefield recorded at two sensors provides the Green's function (GF) at one of the sensors for a virtual source placed at the other sensor (Derode et al., 2003;Lobkis & Weaver, 2001;Snieder, 2004;Wapenaar, 2004, among others). When both recordings are the same then the cross-correlation becomes an autocorrelation and provides a zero-offset GF. The high-frequency zero-offset GF is the reflection response which mainly consists of body wave primary reflections and multiples from seismic discontinuities beneath the station. This principle has been presented already by Claerbout (1968) for one-dimensional (1-D) media as he demonstrated that the full reflection response can be obtained from the autocorrelation of a plane wavefield which is transmitted from below the reflecting structure. Many years later, his finding has been extended to two-dimensional (2-D) and three-dimensional (3-D) acoustic and elastic media (e.g., Wapenaar, 2004).
The numerical determination of the zero-offset reflection response is based on the computation of autocorrelations. For this purpose the continuous noise recordings are cut into short segments to compute the autocorrelations. Finally, the autocorrelations of all segments are being stacked to provide the reflection response. In practice, wave fields are generally not diffuse and the exact GF is unknown or only partly reconstructed. We therefore now refer to empirical GF (EGF) or reflection response rather than GF to emphasize its approximate character.
Note that outliers in the amplitudes of the data can affect correlation results and it is therefore common to preprocess the data before computing autocorrelations (e.g., Bensen et al., 2007;Schimmel et al., 2011). Different strategies exist to balance the amplitudes in the time and frequency domain. One-bit normalization and spectral whitening are often employed methods (Bensen et al., 2007). A disadvantage of these approaches is that they reduce the information content of the signal in a non-unique manner which can cause a less efficient EGF extraction as shown in Schimmel et al. (2018). Here, we avoid amplitude normalization strategies by replacing the conventional cross-correlation and linear stack with the amplitude unbiased phase cross-correlation (PCC, Schimmel, 1999) and time-frequency phase weighted stack (tf-PWS, Schimmel & Gallart, 2007;Schimmel et al., 2011).
In the following, we highlight the main methods used to determine the reflection response at the InSight landing site and describe our data selection and processing.

Correlations
The autocorrelation measures the self-similarity of a time series as a function of lag time. It is a special case of the cross-correlation where both time series are the same. There exist different strategies to measure the self-similarity of time series (e.g., Schimmel et al., 2018). Here, we employ the phase cross-correlation c PC-C (τ) (PCC, Schimmel, 1999) which is based on the phase coherence of instantaneous phases as determined from analytic signal theory. PCC is expressed as follows: SCHIMMEL ET AL.

10.1029/2021EA001755
Φ(t) and Ψ(t) are the instantaneous phases of the two input time series, τ is the lag time of the second time series with respect to the first one, T is the length of the time window, and ν a parameter which permits to tune the sensitivity of the PCC. e iΦ(t) is the envelope normalized analytic signal of the first time series. The analytic signal can be build with the Hilbert Transform and is a unique representation of a real-valued time series in the complex number space. PCC benefits from the fact that thanks to the analytic signal theory a real time series can be decomposed into an instantaneous phase and amplitude (envelope) function. PCC uses only the instantaneous phases and is therefore explicitly signal amplitude unbiased.
The phase autocorrelation is obtained using Φ(t) = Ψ(t). Further, we use ν = 2 since it permits to simplify the equations for fast computations (Ventosa et al., 2019). c PCC (τ) is, in analogy to a classical correlation, a real numbered functional with values ranging between −1 and 1. One of the main differences to the classical correlation is that waveform similarity is measured through the amount of phase-coherent samples rather than the sum of amplitude products.

Stacking
The obtained autocorrelations need to be averaged over larger time spans to achieve a stable noise response. Stacking data over larger time spans improve the azimuthal coverage of the noise wavefield and the cancellation of cross terms (e.g., Medeiros et al., 2015;Snieder, 2004). Both are necessary conditions in seismic interferometry. Here we employ linear stacks and tf-PWS (Schimmel & Gallart, 2007;Schimmel et al., 2011). The tf-PWS are based on the instantaneous phase coherence (in analogy to PCC) which is being measured in the time-frequency domain during stacking. The weights range between 0 and 1 to attenuate less coherent signals in the time-frequency domain. Finally, linear stacks are weighted by the time-frequency phase coherence to build the tf-PWS.
The time-frequency representation of the data is obtained with the S-Transform by Stockwell et al. (1996) which is based on Fourier theory and frequency-dependent Gaussian-shaped windows. The obtained time-frequency representation can be made an analytic signal as shown in Schimmel and Gallart (2007) to enable the use of an instantaneous phase coherence measure which is the backbone of tf-PWS. Alternatively, tf-PWS can be presented using the wavelet transform (Ventosa et al., 2017) to reduce redundancies and to increase computational efficiency. This modification is based on the fact that the S-Transform with Gaussian-shaped windows can be derived from the Morlet wavelet transform (Ventosa et al., 2008). Both strategies provide exactly the same result.

Data and Time Frame
The VBB SEIS data are recorded on three oblique axes. Data from the three axes are needed to build the vertical and horizontal components (Lognonne et al., 2019). For this study, we use vertical components, recorded from May 22, 2019 to June 30, 2020 at 20 samples-per-second and instrument corrected to ground velocity from 0.01 to 10 Hz. The correction is performed for the three oblique axes before rotation to the vertical and horizontal components. The late starting date has been chosen to avoid the initial deployment and testing phase. Key dates related to the seismic part of the mission are listed in Table 1

Data Segmentation and Selection
We filter the seismic recordings to a broad frequency band, 1.2-9.8 Hz, and determine the relative RMS variability for sliding data windows. A maximum RMS variability threshold is then defined to obtain a mask SCHIMMEL ET AL.

10.1029/2021EA001755
and to extract data segments with RMS variability below the threshold through a minimum time duration. We use this procedure to build subsidiary data sets which contain less data problems such as glitches and outlying amplitude events. The approach and the subsidiary data sets are described in the following.
The moving window RMS r is determined for window i and length N using where n stands for the time index and a n for the seismic record. i marks the center sample of the sliding window. In a next step the relative variance s 2 of the RMS amplitudes is computed during a second moving window analysis with window length M > N. The relative variance s 2 of window j is expressed as where R j is the non-zero mean RMS amplitude within the analysis window j. We define the relative variance s 2 as the variance divided by the square of the mean to obtain a dimensionless (scale invariant) measurement of variability. This definition permits to compare data segments with variable RMS mean and to adjust to different seismic background noise levels which may happen along the Sol.
In a next step, data segments are determined with  2 2 j s S , where S 2 is a chosen maximum RMS variability. A new subsidiary data base is then built accepting only data chunks with time duration larger than a minimum time duration, say t min .
Thus, the free parameters of this data selection and segmentation approach are the window length (N and M of Equations 2 and 3) and step intervals (i and j) for the two analysis windows, the maximum RMS-variability S 2 and the minimum data segment length t min .
On one hand, the analysis windows should be several times longer than the largest period as defined by the lowest corner frequency of the bandpass filter to warrant statistical significance. On the other hand, the windows should not be too long to permit tracking RMS-variability changes with time. Throughout this study, we used window length and step interval of 5 and 0.1 s for the sliding-window RMS measurements (Equation 2) and 20 and 1 s for the RMS-variability determination (Equation 3). Further, the shortest admitted data segment length (t min ) is 300 s and the chosen RMS-variability thresholds (S 2 ) are 0.1 and 0.2. Using these parameters two subsidiary data sets have been set up which consist of data volumes of about 3% (RMS-variability threshold 0.1) and 30% (RMS-variability threshold 0.2) of the total data base. In the following, we will refer to these data sets as the total or 100%, 30%, and 3% data sets.
Different other parameters have been tested. The thresholds on data segment lengths and RMS-variability have a direct influence on the amount of the selected data while varying the analysis window parameters SCHIMMEL ET AL.

Table 1 Summary of Key Dates and Works Related to SEIS Deployment and Operation on Mars
within reasonable bounds has a minor or insignificant impact on the subsidiary database. In any case, there exists a trade-off between the selected amount of the data and tolerance of RMS-variability changes. A fine tuning with the parameters may make sense for data sets where outlying noise is clearly separated from the rest of the data. In our case, aseismic noise exists at any amplitude level and therefore, it seems easiest to just analyze the results for more or less restrictive subsidiary data sets.
We expect that a small RMS variability is more compatible with a diffuse wavefield than a high RMS variability. A high RMS variability indicates non-stationarity due to transients (in time localized seismic or aseismic signals) rather than intensity fluctuations in a diffusive wave propagation regime due to scattering and mode conversions. Thus, the benefits of this data selection approach can be multi-fold as we discard data segments with non-stationary noise, ballistic signals, artifacts, and instrument disturbances which may bias the noise correlation function and as the selected data are likely closer to the theoretical requirement of a diffuse wavefield. Our data selection based on RMS variability may also be useful for analyzing the Earth seismic data of poor quality containing many glitches or other artificial transient signals.
Examples of the data selection and segmentation procedure are shown in Figure  2) to find the data segments with low RMS-variability. Segments with a duration larger than t min = 300 s are drawn in red and form part of the 30% data set.
It can be seen from Figure 1 that the selected red segments are data stretches with low-amplitude variability. That is, outlying amplitude signals have been avoided by the approach. Figure 1c also shows that segments with different amplitudes are selected inherent to the use of relative variance 2 j s . Further examples are shown in Figure S1. Figure 2a shows the evolution of the relative RMS variability as a function of Sol. Plotted are mean (red dots), one standard deviation uncertainty (gray bars), minimum, and maximum value (black triangles) per Martian-day. All these values stay stable until about Sol 450. From about Sol 450 a slight but systematic increase of the maximum and mean RMS variability is being observed. A more rapid increase of the minimum RMS variability is evident from about Sol 500 to Sol 530. At the end of the analyzed period all values are increased with respect to the beginning of the studied period. The RMS variability increase correlates with a systematic increase in bad weather attributed to the start of local winter with more storms (Spiga et al., 2020). For example, Figure 2b illustrates the measured mean wind speed (red dots), its one standard deviation uncertainty (gray bars), and minimum and maximum wind speed as function of sol. It can be seen that the minimum wind speed increases fast from about Sol 500 on, which is in concordance with the rise of the minimum RMS variability caused by the presence of winds at the lander site at all times. The highest wind speeds and standard deviation values decrease from about Sol 420 which reflects that the wind characteristics and distribution changed at the lander as also recorded by SEIS.
The just mentioned rise of RMS variability causes a decreased amount of the selected data toward the last analyzed Sols as also documented in Figure S2. This figure shows the LMST of the selected data segments as a function of Sol for the 30% and 3% data sets. It can be seen that most of the selected data segments are from the evening and morning hours. This is expected since lander and wind activity are low during the night. Still, glitches and other aseismic abrupt signals may happen during night which makes our data adaptive selection and segmentation approach different from any time dependent (e.g., day/night) selection. It is also seen that less data has been selected during the last 60 Sols owing to the increase of bad weather during this period.
The data segment length distribution for the 30% subsidiary data set is shown in Figure 3a. One third of the data segments has lengths smaller than 10 min. The number of segments decays quickly as a function of the segment length. Figure 3b illustrates a normalized RMS amplitude distribution. The vertical axis now shows the total time duration rather than the number of individual segments. The gray histogram is for the 100% data set while the red histogram corresponds to the 30% data set. The RMS normalization is the same for both data sets. The absolute RMS amplitude is not important here, but the comparison of both SCHIMMEL ET AL.  histograms shows that the 30% data set consists mostly of small-amplitude segments. The 100% data set RMS amplitudes have been obtained for 2-h windows where each window likely contains sections of high RMS amplitude which influence the RMS of that data window. Conversely, the 30% data set is a subsidiary data set with variable window length as data stretches of high amplitude variability (which means high RMS) are excluded due to the data selection. This exclusion translates into a lowering of the RMS for the 30% data set and explains the increased total time duration for small RMS.

Phase Autocorrelation Spectra for Selected Data Sets
The Fourier Transform of the conventional autocorrelation of a time series equals its energy spectral density (ESD) for a time series which have finite duration and are square integrable. In analogy to the ESD, we compute here the Fourier amplitude spectra of phase autocorrelations. The obtained spectra do generally not equal the ESD of the data as we use PCC to determine autocorrelations. Dissimilarities are often small and inherent to the different measure of signal coherence with each autocorrelation approach. However, outlying amplitude signals do not bias PCC which permits a fast convergence to the EGF and stable autocorrelation spectra as shown in Schimmel et al. (2018). Figure 4 shows the Fourier amplitude spectra of the linear stack of the phase autocorrelations using the 100% (black line), 30% (red line), and 3% (blue line) data sets. The data have been band-pass filtered between 0.8 and 9.5 Hz before computation of the autocorrelations. Different spectral lines and bumps are visible in these spectra. The spectral lines at integer frequencies are related to the "tick noise" which is an aseismic signal caused by electrical coupling during the SEIS temperature acquisition at 1 sps (sample per second). The spectral amplitude bumps are related to lander modes and/or resonance of the seismic subsurface structure (e.g., Ceylan et al., 2021;Dahmen et al., 2021). The main bumps have been labeled along the spectrum for the 100% data set.
Bump 1 and spectral lines at integer frequencies appear to be enhanced in the data sets with less RMS amplitude variability, that is, for the 30% and specially 3% data sets. This is partly due to the attenuation of other signals and noise with higher RMS amplitude variability. The origin of bump 1 is still under investigation.
It is centered at about 2.4 Hz and could be caused by the superposition of lander and structural resonances. The corresponding lander modes at bump 1 are expected to arise as a response of the solar array oscillations during laminar wind flow, that is, at wind speeds below about 3 m/s. This may explain why this bump 1 appears to be enhanced for the data sets with less RMS amplitude variability. Conversely, bumps 2-6 are systematically attenuated in the 30% and 3% data sets. These data sets consist mostly of seismic recordings SCHIMMEL ET AL.
10.1029/2021EA001755 8 of 22 during the local evening hours ( Figure S2) which is when there are less wind and less air pressure perturbations (e.g., Lognonne et al., 2020;Stutzmann et al., 2021). Bumps 2 to 6 correspond to resonances that are related to different lander modes which are mostly excited by turbulent winds above a threshold speed of about 3 m/s (Dahmen et al., 2021). In the following we will work in the 1.2-8.9 Hz frequency band to avoid resonances at 1 and 9 Hz. Figure 5a shows vertical-component phase autocorrelations stacks as a function of the autocorrelation lag time and recording date. The entire data set has been used without any selection criteria, cut into 2-h data segments and frequency band passed between 1.2 and 8.9 Hz before computation of the phase autocorrelations. The autocorrelations have been stacked linearly for recording dates within non-overlapping sliding 3-Sol data windows. In this figure, the positive amplitudes have been filled with red color and the positively correlated signals seen at every full second are tick noise (Figure 5a). We have chosen a window at large lag time to avoid the interference with reflections from shallow discontinuities to permit a visual judgment on the tick noise attenuation. Throughout this work, time is increasing downward as generally done in the field of seismic exploration since travel time is a proxy of depth.

Attenuation of Tick Noise, Lander, and Shallow Structure Resonances
The tick noise is caused by electro-magnetic coupling during the SEIS temperature acquisition at 1 sps with a recorded waveform which does not resemble an impulsive signal (Compaire et al., 2021). On vertical components and within the here considered frequency band, most of its energy is at 4, 7, and 2 Hz (Figure 4).  . Normalized Fourier amplitude spectra of autocorrelation stacks based on phase cross-correlation and linear stacking. Black, red, and blue lines show the spectra for the total data set, the 30% and 3% subsidiary data sets, respectively. Italic numbers 1-6 mark features discussed in the main text.
subsurface layer due to strong impedance contrasts at the limits of the layer. Only little energy leaks to deeper layers which therefore can inhibit the imaging of deeper structures with frequencies in the structural resonance band.
Here, we attenuate the tick noise together with other resonances using band-rejection filters. First, we try two band-rejection filters (3.9-4.4 Hz and 6.8-7.2 Hz) to remove the two strongest tick noise peaks at 4 and 7 Hz and neighboring lander modes (bumps 4 and 6). Figure 5b shows in full analogy to Figure 5a the autocorrelations after the application of the two just mentioned band-rejection filters. Bump 3 has not been included into the first rejection band since this signal is already attenuated for the 30% data set (Figure 4, middle panel).
It can be seen from Figure 5b that there exist still some lower-frequent tick noise which justifies employing a third band-rejection filter from 1.9 to 2.5 Hz. The resulting autocorrelation section Figure 5c is finally not dominated by tick noise anymore. Besides, the data has been cleaned from other resonances (bumps 1, 4, and 6).
An alternative approach is to measure the tick waveform through stacking of 1 s data segments for subsequent subtraction from the data (Compaire et al., 2021). Here we chose a different and independent approach to further show stability in our results.
SCHIMMEL ET AL.
10.1029/2021EA001755 10 of 22  Figure 6 shows autocorrelations computed using PCC and stacked with the tf-PWS approach. The employed seismic recordings belong to the 30% data set. We bandpass filtered these data from 1.2 to 8.9 Hz and applied the three band-rejection filters explained in the previous section. The autocorrelations have been stacked within non-overlapping 3-Sol windows and negative amplitudes are plotted in blue. Note that P-wave reflections at seismic discontinuities with an impedance increase are expected to have negative amplitudes. The autocorrelation section of Figure 6 has been split into three lag-time windows to apply an independent amplitude normalization in each window for visual purposes. The earliest lag time starts at 5 s as the first seconds are dominated by zero-lag autocorrelation sidelobes (e.g., Romero & Schimmel, 2018;Ruigrok & Wapenaar, 2012). The amount of stacked data is shown to the top of each autocorrelation. The maximum duration of the data used per autocorrelation is about 1 Sol rather than 3 Sol which is the length of the data window. This is because of the applied data selection which reduced the data used here to about 30%. The corresponding linearly stacked autocorrelation section is displayed in Figure S4 of the Supporting Information. As also shown further below, the tf-PWS approach provides cleaner results as incoherent signals and noise have been attenuated during the stacking.

Reflection Response and Interpretations
The autocorrelation section of Figure 6 contains signals repeated on most autocorrelations at specific lag times. These signals are stable as they show up for independent autocorrelations along the entire analyzed time interval. One of the most outstanding signals is seen at about 10.6 s. We measured its lag time on the individual autocorrelations ( Figure 6) and determined a mean lag time and standard deviation of 10.612 and 0.054 s, respectively. Similar values are found for the 100% (10.627 s, 0.061 s) and 3% data sets (10.63 s, 0.04 s) using the processing of Figure 6. Standard deviations are in the order of the 0.05 s sample interval.
If we interpret this arrival as a P-wave reflection and assume an average crustal P-wave velocity of 4 km/s then the reflector is expected to be at about 21 km depth. Note that the assumed average P-wave velocity is slower than for Earth as indicated by receiver function inversions (Knapmeyer-Endrun et al., 2021; Lognonne et al., 2020).
As the raw data consists also of different other aseismic signals, an autocorrelation signal is not a synonym of a seismic reflection. Some of the signals are expected to be part of the seismic reflection response while other signals can be due to any coherent and repeated feature in the data or due to beating of resonating frequency components from the lander. Resonances with slightly different frequencies f 1 and f 2 > f 1 cause an interference pattern (wellknown in music) with beat frequency f b = f 2 − f 1 (e.g., Kinsler et al., 1999). That is, the addition of two cosine waves of different frequency is equivalent to a cosine wave with carrier frequency 0.5(f 1 + f 2 ) whose amplitude varies at the frequency 0.5(f 2 − f 1 ). The envelope modulation, in the sense of intensity of the wave, oscillates at f b = f 2 − f 1 and is in music the frequency at which beats are heard. The beats of the interfering signals manifest in time-domain autocorrelations as distinct and repeating signals of maximum amplitude at lag times t n = n/f b with n ∈ N which bring the beats to overlap. Similarly, glitches or other signals which are repeated systematically with the same time interval are expected to show up in autocorrelations at lag time which corresponds to the time separation of the repeated features.
Indeed, the interpretation of autocorrelations without any further a priori information is difficult. However, receiver function studies of marsquakes (Knapmeyer-Endrun et al., 2021;Lognonne et al., 2020) provide the first a priori information of the underlying structure and expected crustal-scale P-wave reflections on Mars. In this respect, the strong arrival at 10.6 s lag time interpreted as a P-wave reflection seems to be consistent with P-to-S wave conversions observed with receiver functions (Knapmeyer-Endrun et al., 2021). Indeed, the 10.6 s signal is quite robust and a likely candidate for a P-wave reflection from the crust-mantle boundary. Note that Compaire et al. (2021) also report a signal at 10.6 s using a different approach.
The 10.6 s signal is discernible in linear and tf-PWSs of the three data sets as shown in Figure 7. This figure contains the linear stacks in black and the tf-PWSs in red for all autocorrelations from the 100% (top traces), 30% (middle traces), and 3% (bottom traces) data sets. The traces have been normalized by their RMS amplitudes measured between 8 and 20 s lag time. Further, amplitudes are clipped during the first seconds to favor the visibility of later arriving signals. It can be seen that the 10.6 s signal is visible with amplitude of about twice the amplitude of neighboring signals and noise. Its shape resembles more a zero-phase wavelet Figure 6. Vertical-component autocorrelation stacks for sliding 3-Sol data windows. The frequency band is 1.2-8.9 Hz and data windows do not overlap. Shown are time-frequency phase weighted stacks of phase autocorrelations computed for the 30% data set. Blue marks negative amplitudes. The three lag-time windows have been used to improve the visibility through independent amplitude normalization. Top panel shows the total duration of the selected data used to compute autocorrelations within each of the 3-Sol data windows. Linear stacks for the same data are shown in Figure S4. than a wave train. The direct comparison of the stacks also shows that the tf-PWS approach provides cleaner autocorrelations through the attenuation of incoherent signals and noise.

Convergence
We now look at the convergence of autocorrelation stacks toward a stable seismic noise response. In seismic interferometry, the minimum amount of data needed is controlled by the cancellation of noise-cross terms (e.g., Snieder, 2004) and the abundance, distribution, quality, and duration of signals to build up a stable EGF.
Under the assumption of stationarity and minimum frequency of 1 Hz, cross-terms cancel out within the confidence level ϵ ≈ 0.01 after about 3 h (Equation 9 in Medeiros et al., 2015). The SEIS data are highly variable and the assumption of stationarity is not valid. Nevertheless, the estimated 3 h can be taken as a guide value. This value is shorter than the used data length of each autocorrelation shown in the figures for the 30% and 100% data sets.
The total data length needed to extract a robust reflection response can be estimated by a convergence analysis. In this approach the evolution of autocorrelation stacks of random subsidiary data sets is being inspected as function of amount of data used for the autocorrelations in each of the stacks. For this purpose the individual autocorrelation stacks are being compared with a reference waveform obtained by stacking all available data. The waveform similarity of the two stacks can be determined with the zero-lag cross-correlation of the two time-series. Here, we employ PCC and the geometrically normalized cross-correlation (CCGN, e.g., Schimmel, 1999). The reference traces are those shown in Figure 7.
There exist different data sampling strategies to draw the random autocorrelation stacks. We tested two strategies, random sampling without replacement and Bootstrapping which uses random sampling with replacement (Efron & Tibshirani, 1986). Without replacement means that no autocorrelogram is used more than one time in a draw. We repeated the random drawing 40 times to enable the determination of mean and standard deviation for the similarity to the reference waveform. Figure 8 and Figure S5 demonstrate our results using the Bootstrapping resampling, PCC, and CCGN, respectively, to measure waveform similarity in moving windows of 0.5 s length. The 30% data set has been used for Figures 8a and 8b. Figure 8a shows the waveform similarity while Figure 8b contains the standard deviation of the similarity. The similarity (red line) and corresponding standard deviation (gray error bars) at 10.6 and 19.5 s are shown in Figures 8c and 8d. Nineteen and five tenths seconds has been chosen arbitrarily for comparison. A plateau in the similarity curve means that stacking more autocorrelations is not changing significantly the waveform similarity. The figures illustrate that the 10.6 s signal converges faster than a signal at 19.5 s lag time. Using CCGN ( Figure S5) an apparent overall faster waveform similarity is measured than using PCC. This is expected as PCC is the more waveform sensitive measure due to the employed instantaneous phase coherence.
From Figure 8a and Figure S5a we also observe that signals at shorter lag time tend to converge faster than signals at larger lag time. For a P-wave reflection response a larger lag time means longer wave paths owing to a deeper discontinuity or multiple reverberation. Such signals may need more data to stabilize in the EGF since the corresponding waves are more susceptible to the different types of attenuation.
A comparison of the similarities for the 100%, 30%, and 3% data sets is shown in Figures 8e and 8f and Figures S5e and S5f. It can be seen that slowest signal convergence is achieved using the 100% data set (black curves). This means that a high amplitude variability is not contributing significantly to the extraction of the EGF. Note that this result does not change when using only the data until Sol 410. That is, the increasing RMS variability after Sol 450 (Figure 2) is not causing the slower convergence.
The 3% and 30% data sets seem to provide more similar results. The 3% and 30% data sets have a total length of 151.8 and 1557.1 h. The 3% data set results are therefore obtained with a higher replacement rate. SCHIMMEL ET AL.  The linear stacks are plotted in black and the corresponding time-frequency phase weighted stack in red. The amplitudes of each stack have been normalized for visual purposes. The 10.6 s signal (marked by an arrow) appears with in impulsive shape for the three data sets and stands out from its immediate surroundings. Figure 8. (a) Autocorrelation waveform convergence as function of lag time and total length of data using the 30% data set. Colors are used to contour the mean similarity of randomly stacked autocorrelations with the reference trace. The reference trace is the linear stack of all autocorrelations from the 30% data set (black trace, middle panel of Figure 7). The total data length is 2350.2 h. (b) Standard deviation of the similarity shown in (a). (c) Similarity (red) and standard deviation (gray) of the 10.6 s signal for the 30% data set. (d) Same as (c) but for a signal at 19.5 s lag time. (e) Same as (c) but showing also the results for the 3% (blue) and 100% data sets (black). The reference traces are based on a total data length of 201.2 and 8355.6 h for the 3% and 100% data sets, respectively. (f) Same as (e) but for the 19.5 s signal.
Drawing the random data sets without replacement (not shown here) illustrates that the similarity curves for the 3% and 30% data are still more similar than for the 100% data set.

Stability of Reflection Response Over Time
The autocorrelations for the three data sets are shown in Figure 9. The only differences to Figure 6 are the respective data sets (from left to right: 100%, 30%, and 3%) and the use of non-overlapping 30-Sol windows rather than 3-Sol windows. The top panels show the amount of the data used for each of the stacks. The vertical axes need to be multiplied by 0.3 and 0.03 to obtain the duration in Sol for the 30% and 3% data sets. 30-Sol windows have been employed to reduce the number of autocorrelations to permit a comparison of SCHIMMEL ET AL.  Data set: 100% x 0.3 30% x 0.03 3% the three autocorrelation sections at one glance. A welcomed side effects is that these autocorrelations are more robust resulting from the increased amount of the data used in each of the stacks. The corresponding linear stacks are demonstrated in Figure S6 for completeness.
The most striking observation is that the 100% autocorrelations show systematically changing noise responses, that is, the autocorrelations from about Sol 450 toward the end of the section look different with respect to the autocorrelations obtained for the first Sols. This variability is likely due to changing climatic conditions with stronger storms and corresponding responses of the lander. In fact, the statistical properties of the seismic recordings start to change simultaneously at about Sol 450 as testified with the observed daily RMS variability shown in Figure 2. In contrast to the 100% data, the 30% autocorrelations (middle panel of Figure 9) show stable signals for all Sols. This is because the RMS variability threshold retained most of the seismic recordings with imprints of climate and lander response variability. This result reinforces that a thorough data selection is essential to warrant a minimum influence of aseismic signals which may bias any autocorrelation result. In this respect, a data adaptive approach is better than simply selecting the data through LMST slots.
Amplitude spectra of the just discussed autocorrelations are shown in Figure 10. The amplitudes have been normalized at 6 Hz for visual purposes. Spectra for linearly stacked autocorrelations are shown in Figure S7 for comparison. It can be seen from both figures that the 100% data set suffers at about 3.8 Hz and Sols larger than 500 from energy which appears at the corner of one of the applied band-rejection filters. These signals are more pronounced in the spectra of the linear stacks and are likely lander resonances in response to the bad weather. Amplitude spectra for linearly stacked autocorrelations computed for data without application of band-rejection filters ( Figure S8) reveal that the strongest amplitude changes indeed occur around 4 Hz.
SCHIMMEL ET AL.

10.1029/2021EA001755
16 of 22 Figure 10. Amplitude spectra of time-frequency phase weighted stack autocorrelations shown in Figure 9 as function of Sol for the 100% (left panel), 30% (middle panel), and 3% (right panel) data sets. Each amplitude spectrum is placed at its window center time and has been normalized at 6 Hz. Spectra of linear stacks are shown in Figure S7.
Broader band-rejection filters could attenuate these signals, but unnecessarily reduce the band width of the autocorrelations for the data recorded before Sol 450.
As expected from the time-domain autocorrelations of Figure 9 the corresponding spectra ( Figure 10) for the 30% and 3% data sets show less variability than for the 100% data set. The strongest resonances have still been attenuated by the data selection approach and three band-rejection filters.

Stability of Reflection Signals Over a Broad Frequency Band
Finally, seismic signals related to structural discontinuities are expected to be detected for a rather broad than narrow frequency band, that is, possible frequency-dependent reflectivity generally varies smoothly.
Besides, at small lag time, sidelobes of the zero-lag autocorrelation peak may affect early structural arrivals. The autocorrelation sidelobes are caused by the autocorrelation of the noise source time function (e.g., Ruigrok & Wapenaar, 2012) and are expected to change faster with frequency than a reflection pulse (see Figure 6 in Romero & Schimmel, 2018).
We therefore inspect the signal stability across different frequency bands. The first three panels in Figure 11 show autocorrelations for the 30% data set and different one-octave frequency bands. From top to bottom these are 1.5-3.0 Hz, 2.4-4.8 Hz, and 3.6-7.2 Hz. The effective bandwidth, nevertheless, is smaller than an octave owing to the employed band-rejection filters. The tf-PWS of all autocorrelations is illustrated in red while the black lines show tf-PWSs of 10% randomly drawn autocorrelations. The fourth panel shows the red color stacks of the top three panels plotted on top of each other. The lowermost panel contains the tf-PWS (red line) and linear stack (blue line) for the wide 1.2-8.9 Hz band.
It can be seen from this figure that some signals are persistent or in a phase along different frequency bands. Labels a-f mark some of these signals at lag times of 6. 15, 6.85, 8.66, 10.6, 16.85, 23.8, and 24.49 s. The dashed vertical lines have been added to aid the visual inspection. Signals a and c are not seen for the lowest frequency band and have therefore been marked at the base of the figure. None of the signals is being observed for only one frequency band. For example, signal b is in phase for the three frequency bands while signal d is in phase for the two lower bands. The latter signal, however, appears only in the middle frequency band as a clear signal. It is seen from the top two panels that the amplitude maxima of signal d are not coherent while its minimum amplitude occurs at the same lag time to build up the observed negative polarity when filtering over a broader frequency band. Features a and c are in phase for the two higher frequency SCHIMMEL ET AL. bands. Signals e and f are also in phase and show also up as slightly distinctive features in the 1.2-8.9 Hz correlations. e and f are less impulsive as d, that is, in the 1.2-8.9 Hz correlations they seem to be part of a larger waveform. F has been marked twice, at the first minimum amplitude and 0.6 s later at the lowest minimum. We use only one letter, as it is not clear if this are several signals which arrive at a similar time.
In any case, a visible interference of zero-lag sidelobes or their misinterpretation is not expected owing to the stationarity of these signals over a broader frequency band (Romero & Schimmel, 2018). That is, zero-lag sidelobes are expected to be shifted in time for different frequency bands while reflections from sharp discontinuities appear at the same time lag.

Discussion and Conclusion
The 10.6 s signal, if interpreted as P-wave reflection, would point to a discontinuity at about 21 km depth considering expected Mars's crustal average seismic velocities of 4 km/s. We observe the 10.6 s signal robustly on the 100%, 30%, and 3% data sets with amplitudes which can reach twice the amplitude of other nearby signals and noise. The shape of the waveform resembles more a zero-phase wavelet with negative polarity than a wave train. The negative polarity is consistent with a reflection from a discontinuity with a seismic impedance increase while the wavelet shape can be explained by a sharp discontinuity and lack of other interfering signals as reflections from strong nearby discontinuities.
Indeed, the 10.6 s signal is a likely candidate for a reflection from the base of the crust owing to its two-way travel time, negative polarity, strengths, and robustness. Energy at similar lag time has also been reported by Deng and Levander (2020) and Compaire et al. (2021) and interpreted as P-wave Moho reflection. Compaire et al. (2021) also observe energy on autocorrelations for the North components in the frequency band around 2.4 Hz which could be S-wave reflections from the same discontinuity if one assumes a P-to-S wave velocity ratio of 2.1.
The 10.6 s signal seems to be consistent with models obtained for receiver functions from the marsquake data (Knapmeyer-Endrun et al., 2021;Lognonne et al., 2020). Knapmeyer-Endrun et al. (2021) present two model families resulting from receiver function inversions where the 10.6 s signal can be explained either as a P-wave Moho or internal crustal discontinuity reflection. More events from different distances or other geophysical constraints are needed to further limit the model space and to aid the identification of autocorrelation signals. We prefer the first model class which can explain the 10.6 s signal as a Moho reflection since we observe no other, at least equally strong signal at a later time to explain a deeper Moho. This is based on the assumption that the crust-mantle boundary is stronger than a mid-crust discontinuity.
The seismic data acquisition conditions are harsh on Mars and there exist a wealth of different aseismic signals Lognonne et al., 2020;Scholz et al., 2020;Stutzmann et al., 2021) which also may show up on autocorrelations (Kim et al., in preparation). For instance, regularly repeated glitches (long duration pulses typically with frequencies below 1 Hz) and donks (short duration pulses typically at frequencies above 8 Hz, Ceylan et al. 2021) with similar waveforms can manifest, depending on their abundance and coherence, at their repeat time. To attenuate the effect of glitches and donks, Deng and Levander (2020) determine the absolute mean amplitude in moving data windows of 100 s length which then are used to weight the seismic data (e.g., Bensen et al., 2007). Finally, they observe at 11.5 s a signal which appears robustly on their Martian daytime and night-time data and which they interpret as Moho P-wave reflection.
Conversely, Compaire et al. (2021) removed the tick noise and attenuated glitches using an algorithm by Scholz et al. (2020). High-frequency components of glitches remain in the data and the algorithm does not remove donks. Overall, they use data recorded during the evening hours for which autocorrelations have the highest signal-to-noise ratio (SNR) and find a signal at 10.6 and at 21 s ,which they interpret as reflections from crustal layers. They also show that glitches cannot interfere with these reflections as they are separated by more than 30 s. The repeat time of donks, however, can be much smaller with maximum at about 10 s around 17-18 h LMST (Compaire et al., 2021). The authors argue that the observed stability of the 10.6 s signal with time is not correlated with the overall distribution of donks, which finally makes a donk origin of these reflections unlikely.
Here, we use a different data processing strategy. For instance, we do not de-glitch the records as high frequency glitch signatures seem to remain in the cleaned data. Further, we refrain from selecting data by choosing fixed time slots and design a data-adaptive approach to segment and pick seismic records based on their RMS amplitude variability. Our 30% and 3% data sets have been built selecting stretches of data with low RMS variability which also led to an attenuation of most lander resonances as seen from the spectra in Figure 4. We further attenuated strong resonances using three band-rejection filters (1.9-2.5, 3.9-4.4, and 6.8-7.2 Hz). These band-rejection filters were also chosen to remove jointly the strongest components of the tick noise.
Our data selection procedure is expected to avoid large amplitude donks and glitches and should be less affected by aseismic features, which are related to bad weather lander and SEIS responses. However, small amplitude donks which do not affect the ,RMS variability can still be present in the selected data. Concerning the 10.6 s signal, a donk origin would require a sharp repeat time distribution of opposite polarity donks to explain the negative correlation peak and zero-phase shaped wavelet. A smooth distribution of donks is expected either to cancel out or to produce wave train shaped signals at lag times which correspond to their repeat times. In any case, at present we are not aware of hidden donks, but can also not rule out their existence.
Finally, large amplitude donks and glitches do not seem to have a dominant imprint on the autocorrelations as can be seen by comparing the autocorrelation sections for the 100% and 30% data sets ( Figure S9). This is because randomly distributed or incoherent glitches and donks, as well as other outlying signals and noise do not affect the amplitude unbiased PCC. Therefore, PCC is an adequate approach for the data from unfavorable deployment and detection conditions as shown in Schimmel et al. (2018).
Another source of aseismic signals can be beats produced through interfering resonances (Kinsler et al., 1999). Beats can manifest in autocorrelations at lag time, which equals the inverse of the difference in frequency of the two resonances and multiples. That is, a 10.6 s autocorrelation signal could be due to a 0.094 Hz beat. The 10.6 s signal does not show a seasonal variability and would therefore require stable resonances throughout the analyzed data. We believe that the autocorrelations are not dominated by strong beats as the main lander resonances have been attenuated through rejection filters and data selection on frequency band pass filtered data. Still, the presence of weak amplitude beats in the data cannot be ruled out.
Further, the fast convergence of the 10.6 s signal which can be extracted with about 2 Sol of the data is remarkable. There are no oceans on Mars, and all seismic energy is supposed to be released by weak magnitude events (the largest observed Marsquake has a magnitude of 4; Clinton et al., 2020) and through atmospheric phenomena. In spite of the absence of strong sources, the low level of natural Martian noise  may favor the detectability of low-amplitude signals and therefore aid the convergence of the reflection response.
We focused the discussion on the 10.6 s signal, as it appears to be a likely candidate for a reflection from the base of the crust. Nevertheless, the obtained reflection response contains also other signals which likely are reflections and multiples from other discontinuities. Some of them have been marked in Figure 11 (labels a to f) owing to their appearance over a broader frequency band. The lag times are 6.15, 6.85, 8.66, 10.6, 16.85, 23.8, and 24.49 s. Specifically for the earlier arrivals the frequency stability is an important criterion to avoid misinterpretations or interferences with zero-lag autocorrelation sidelobes as shown in Romero and Schimmel (2018).
Concerning inter-crustal discontinuities, receiver function studies (Knapmeyer-Endrun et al., 2021;Lognonne et al., 2020) provide stable estimates on a discontinuity at about 8-11 km depths with S-wave velocities Vs of about 1.7-2.1 km/s. If we assume a mean discontinuity depth of 9.5 km and that our signal a at 6.15 s lag time ( Figure 11) is a reflection from this discontinuity then the expected P-wave velocity Vp and Vp/Vs ratio for the 9.5 km layer are 3.1 km/s and 1.48-1.82 considering the Vs range. Under same assumptions, the close by signal b at 6.85 s lag time, would lead to Vp = 2.8 km/s and Vp/Vs ratio of about 1.33-1.65. We believe that signal a provides a more realistic Vp/Vs ratio than signal b and could therefore likely be the reflection from the base of the upper crust layer. The lower Vp/Vs bound of 1.48 is very low with respect to expected values. For comparison, the Vp/Vs ratio for the very shallow regolith layer has been estimated to 1.67 . Nevertheless assuming the upper bound of 11 km for the intercrustal discontinuity (Knapmeyer-Endrun et al., 2021;Lognonne et al., 2020) and our 6.15 s lag time leads to a more realistic Vp/Vs range of 1.7-2.1 based on the Vs values estimated from receiver functions.
Finally, the estimated lag time dependent convergence of the autocorrelations shows that more data are needed to build the reflection response at later lag times. This makes sense, as waves propagating over a longer time interval are more prone to scattering and attenuation which needs to be compensated by stacking more autocorrelations. In any case, we see no reason for aseismic signals to show a systematic lag time dependent convergence.