Autocorrelation of the ground vibrations recorded by the SEIS-InSight seismometer on Mars

32 Since early February 2019, the SEIS seismometer deployed at the surface of Mars

manuscript submitted to JGR: Planets Abstract Since early February 2019, the SEIS seismometer deployed at the surface of Mars in the framework of the InSight mission has been continuously recording the ground motion at Elysium Planitia.In this work, we take advantage of this exceptional dataset to put constraints on the crustal properties of Mars using seismic interferometry (SI).To carry out this task, we first examine the continuous records from the very broadband seismometer (SEIS-VBB).Several deterministic sources of environmental noise are identified and specific pre-processing strategies are presented to mitigate their influence.Applying the principles of SI to the single-station configuration of InSight, we compute, for each Sol and each hour of the martian day, the diagonal elements of the time-domain correlation tensor of random ambient vibrations recorded by SEIS.A similar computation is performed on the diffuse waveforms generated by more than a hundred Marsquakes.
A careful signal-to-noise ratio analysis and an inter-comparison between the two datasets suggest that the results from SI are most reliable in a narrow frequency band around 2.4Hz, where an amplification of both ambient vibrations and seismic events is observed.The average autocorrelation functions (ACFs) contain well identifiable seismic arrivals, that are very consistent between the two datasets.Interpreting the vertical and horizontal ACFs as, respectively, the P-and S-seismic reflectivity below InSight, we propose a simple stratified velocity model of the crust, which is mostly compatible with previous results from Receiver Function analysis.Our results are discussed and compared to recent works from the literature.

Plain Language Summary
The correlation of seismic records is the basis of seismic interferometry methods.
These methods use seismic waves, either from ambient vibrations of the planet or from quakes, that are scattered in the medium in order to recover information about the structure between two seismic sensors.The method is implemented to compute the auto-correlation functions of the three components of the ground motion recorded by the SEIS seismometer.The comparison of the results obtained from earthquake data to the ones obtained from ambient vibrations demonstrates that the ambient seismic vibration is clearly above the self-noise of SEIS during early night hours around a specific frequency (2.4Hz).The seismic vibrations appear to be amplified at this frequency by an unknown mechanism.Some seismic energy arrivals appear consistently in the auto-correlation functions, at specific propagation times, independent of the data sets and processing parameters tested.
These arrivals are interpreted as vertically propagating seismic waves which are reflected on top of crustal layers.Their propagation times can be used to constrain a model of Mars crustal structure.

INTRODUCTION
NASA's InSight (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport) mission landed on November 26, 2018 near the Martian equator in Elysium Planitia (Banerdt et al., 2020).The seismological records provided by the main instrument -SEIS (Seismic Experiment for Interior Structure)-constitute a dataset of unprecedented quality for planetary seismology.The seismometer sensitivity and the extremely low amplitude of the ambient ground vibrations make the martian seismic signals very different from those collected on Earth (Lognonné et al., 2019(Lognonné et al., , 2020)).Not all the features of the ground velocity records can be interpreted as seismic signals, and all the details of the InSight mission system have to be taken into account to correctly interpret the data.Among these, we mention the lander-related noise, the electrical noise, the atmospheric noise and all the mechanical resonances that are temperature dependent (Murdoch, Mimoun, et al., 2017;Murdoch, Kenda, et al., 2017;Lognonné et al., 2020;Garcia et al., 2020).Henceforth, we refer to these as "environmental noise".
Planetary bodies like the Moon or Mars have a much lower seismicity-rate compared to Earth due to the absence of plate tectonics.As a consequence, the scarcity and low amplitude of seismic sources can be an issue for seismological applications.Fortunately, the last twenty years have seen the rapid development of new processing methods which allow seismologists to extract meaningful seismological signals from passive records of random vibrations generated by natural processes such as winds, earthquakes, etc. (Nakata et al., 2019).These methods are often referred to collectively as "seismic interferometry".They are based on the close relation between the mean correlation function C AB of diffuse seismic wavefields recorded at any two stations A and B and the Green's function G AB between A and B (Lobkis & Weaver, 2001;Campillo & Paul, 2003).We recall that the Green's function is the seismic response recorded at A due to an impulsive source operating at B. Because seismic wavefields are rarely perfectly diffuse, the signal reconstructed by cross-correlating random seismic signals is referred to as "Empirical" Green's function as it will in general deviate from the exact Green's function of the medium (Weaver & Lobkis, 2005).We will also adopt this terminology.
Seismic interferometry has been applied with success to Lunar seismic data (Larose et al., 2005).These authors demonstrated the possibility to virtually reconstruct propagating Rayleigh waves by cross-correlating ambient vibrations recorded on an array of geophones deployed in the framework of the Apollo Lunar Seismic Profiling Experiment.
Using the dispersion properties of surface waves, they were able to put constraints on the shear wave speed profile in the Lunar regolith, thereby gleaning new insights from an already prolific dataset (Garcia et al., 2019;Nunn et al., 2020).The experimental conditions of InSight are more specific in the sense that only a single (6-axis) sensor is at our disposal.In this situation, the most direct observables are Auto-Correlation Functions (ACF).This is by no means a severe limitation.Indeed, it is known since the works of Claerbout (1968), that the reflection response of a stack of layers located beneath a single seismometer can be retrieved from its transmission response.In other words, if the medium is illuminated from below, computing the ACF of the transmitted wavefield is essentially equivalent to virtually activating a source at the location of the seismometer and recording the seismic response.The results of Claerbout (1968), valid for a 1Dhorizontally layered medium, have been extended to 3D-inhomogeneous media by Wapenaar (2003), based again on a perfect illumination of the structure.Following the pioneering works of Claerbout (1968) andWapenaar (2003), a number of studies reported the successful passive reconstruction of the crustal reflection response from the AC of field data acquired on Earth (Tibuleac & von Seggern, 2012;Ito & Shiomi, 2012;Gorbatov et al., 2013;Kennett et al., 2015;Saygin et al., 2017;Pha . m & Tkalčić, 2017;Oren & Nowack, 2017;Romero & Schimmel, 2018) and on the Moon (Nishitsuji et al., 2016).Suemoto et al. (2020) were the first to extract propagating signals from AC of ambient vibrations recorded by SEIS.They considered high frequencies and short propa-gation times.They observed stable phases on the ACF of the Vertical, East and North components of the seimic wavefield.From this observation, they deduced the presence of a very shallow interface underneath InSight.More recent works by Deng and Levander (2020) focused on the deep seismic structure of Mars.Their results will be discussed and compared to ours in section 4. In the present work, we apply the seismic interferometry method to two types of diffuse wavefields -ambient vibrations and the diffuse part of seismic events-and investigate the possibility to retrieve the deep vertical reflection response of the Martian crust beneath InSight.In section 2, we describe the particular features of the continuous records of SEIS-VBB and the various types of seismic events.
In section 3, we present a pre-processing strategy to mitigate environmental noise and compute empirical Green's function (EGF) by correlating records of both ambient vibrations and diffuse seismic events.We show excellent agreement between the two datasets in a specific frequency band.Based on these results, some simple layering models below InSight are proposed in section 4, and compared to previous results from the literature.
Section 5 summarizes our main findings and proposes directions for future works.

DATA
2.1 Overview of the SEIS instrument SEIS is the main instruments of the InSight mission (Banerdt et al., 2020).It is a six axes seismometer composed of three Very Broad Band (VBB) sensors sensitive to frequencies between 0.01Hz and 10Hz and three Short Period (SP) sensors sensitive to frequencies between 0.1Hz and 50Hz (Lognonné et al., 2019(Lognonné et al., , 2020)).The three VBB sensors are installed in a vacuum sphere as a first thermal protection.The sphere and the three SP sensors rest on a three-legged levelling system (LVL) that couples with the Martian ground.The whole constitutes the sensor assembly (SA).The SA is connected to the lander via the tether.A relaxation loop (LSA) has been installed at the junction between the SA and the tether.Thus, when the tether contracts or expands under the effect of temperature the movement is not transmitted to the SA.SEIS was deployed on the Martian surface with the robotic arm of the InSight lander on December 19, 2018.
It is covered by the Wind and Thermal Shield (WTS) since February 2, 2019 (Lognonné et al., 2020).SEIS is continuously recording ground motion at the InSight's landing site.
We summarize the acronyms relevant to this study in Table 1.The orientations of the six sensors of SEIS in the Martian geographical system are extracted from SEIS dataless information and provided in Table 2.More information on the location and orientation of the various instruments is provided by Golombek et al. (2020).In this study we focus our analysis on the frequency band below 10Hz and chose to use VBB data at 20 sample-per-second (sps) (InSight Mars SEIS Data Service, 2019).The continuous records of SEIS present several features that have been described in Lognonné et al. (2020) and Ceylan et al. (2021).We briefly summarize the main characteristics here.We show in Figure 1 from the fact that the lander is assembled from different mechanical parts.The various components resonate in different frequency bands, that depend in particular on their elastic properties (Murdoch, Mimoun, et al., 2017).Hence, the periodic modulation of the resonance spectrum is a direct consequence of the thermo-elastic response of the lander induced by daily temperature variations.Finally, we note that the vast majority of lander modes are found in the 1-50Hz frequency range and that most of them produce strongly polarized signals that are detected by SEIS-VBB.As a consequence, the lander modes rarely contaminate simultaneously all components of the signal.
The spectrogram of Figure 1.A also shows a sharp peak of energy at exactly 1Hz frequency.This feature is more clearly visible on the mean power spectrum shown in Fig- ure 1.B.In addition to the main peak at 1Hz, we observe various harmonics at 2Hz, 3Hz, etc.These harmonics are not all visible on the spectrogram due to the low image resolution.The signal at the origin of the peaks at 1Hz, 2Hz, ... seen on the spectrogram is called "tick noise" and corresponds to a periodic cross-talk induced by the temperature acquisition at 1sps.Because the temperature sensors and the seismometer axes share the same clock, the cross-talk signal (tick noise) is perfectly synchronized with the acquisition of the SEIS records.
Figure 1 (black rectangles) also reveals a permanent excitation of a continuum of frequencies between 2Hz and 3Hz, that are not modulated in time in sharp contrast with the Lander modes.This feature is called the "2.4Hz resonance" and is interpreted by Giardini There also exist two families of transient perturbations, presented in Figure 2, that deserve attention.One of these transients dominates the low frequency band and is referred to as "glitch" (panels 2.A and 2.B).The other dominates the high frequency band and is referred to as "donk" (panels 2.C and 2.D).The most likely origin of these signals is the activation of presumably pre-existing cracks in the various mechanical parts of the InSight station including the lander, the Sensor Assembly (SA) of SEIS and the tether between SEIS and the lander (Ceylan et al., 2021).Indeed, the temperature variations between day and night at the InSight landing site can reach 100 K.All the mechanical parts are thus subject to high thermal stresses and the elastic energy accumulated is partly released in the form of low-and high-frequency acoustic emission.For a detailed analysis of glitch signals we refer the reader to Scholz et al. (2020).Because glitches, tick noise and donks can adversely affect the results of seismic interferometry, section 3.1 describes several data procedures to mitigate their effect.

Seismic Events
Since the beginning of the operation phase SEIS has recorded more than a hundred seismic events (Giardini et al., 2020).The event nomenclature used in this study refers to the catalogue provided by the MarsQuake Service (MQS) (InSight Marsquake Service, 2020) described in Giardini et al. (2020).The events are classified into five types and have an assigned quality score from A (best) to D (worst).The description of the five event types is presented in Table 3. Up-to-date information on SEIS events and the Mars quake catalog are provided by Clinton et al. (2021).
In this study, we focus our analysis on the High Frequency events (HF), Very High Frequency events (VF) and 2.4Hz events.As described by Lognonné et al. (2020) and Giardini et al. (2020), the high-frequency seismic events (HF, VF and 2.4Hz events) waveforms have a diffusive character.The abundance and long duration of high-frequency seismic events offer a unique opportunity to apply seismic interferometry to the retrieval of the local seismic response below InSight.All in all, we have used forty-one HF, fourteen VF and sixty-nine 2.4Hz events.Only seismic events with quality better than or equal to C were selected.The complete list of events used is given in Appendix A. For all these events we were able to exploit the entire waveform as previous polarization and multiple-scattering analyses strongly suggest that the diffusive regime sets in almost immediately after signal onset (Lognonné et al., 2020).
Table 3. Characteristics of the different types of events as defined by the MQS (InSight Marsquake Service, 2020)

Event type Description
Low Frequency (LF) Energy exclusively below the 2.4Hz resonance Broadband (BB) Excite the 2.4Hz resonance but with the major part of their energy at lower frequencies 2.4Hz event Excite only the 2.4Hz resonance High Frequency (HF) Excite the 2.4Hz resonance but also the higher frequencies Very High Frequency (VF) Same as HF events but with larger energy on the horizontal components than on the vertical component

PROCESSING AND RESULTS
3.1 Pre-Processing

Tick Noise Removal
The tick noise is an electrical disturbance (cross-talk) resulting from the acquisition of the temperature inside SEIS.The frequency content of this deterministic noise depends mostly on the sampling rate of the temperature sensor.It is important to note that the tick noise waveform differs on the three components U, V and W of the VBB.
As the acquisitions of the temperature and seismic channels are synchronized, the tick noise repeats periodically every N samples, where N is the sampling rate of the SEIS channel.This implies in particular that it is not sensitive to the temporal drift of the SEIS clock.As observed in Figure 1, the tick noise peaks in the frequency domain at 1Hz but also exhibits a non-negligible amplitude at each harmonics (2Hz, 3Hz, ...).In order to remove this noise from the raw data at 20sps, we stack non-overlapping contiguous 20-samples windows.Under the assumption that the background noise is random, the stack should converge towards a good estimate of the tick noise waveform.To obtain the estimates presented in Fig. 3.A, we stacked two months 20sps VBB records acquired between 18:00 and 22:00 Local Mean Solar Time (LMST).We only employ evening data because they sample the quietest period of the day on Mars.Finally, to remove the tick noise from the data we first determine its lag-time with respect to the raw time-series using a simple cross-correlation algorithm.We then appropriately shift the tick noise to align it with the raw time-series and remove the former from the latter in contiguous windows of 20 samples by subtraction.Figure 3.B shows the spectrograms before and after the tick noise removal around a frequency of 1Hz for each component.Figure 3.C and 3.D show the time-series and the spectral contents of VBB-V during the evening of Sol 183 before and after the tick noise removal.We see that our procedure successfully removes the tick noise at the target frequency, while preserving the rest of the spectral content.

Glitch suppression
As described in Lognonné et al. (2020), the term "glitch" refers to a particular type of signal in the seismic channels of SEIS whose waveform is similar to the response of the instrument to a step in acceleration.Glitches are broadband signals but most of their energy is localized in the low frequency domain (< 1Hz).Rarely, glitches can be preceded by a high-frequency precursor.Glitch amplitude extends over six orders of magnitude.They may happen anytime during a Sol but those with highest amplitudes appear to occur when the temperature exceeds certain specific values (Scholz et al., 2020).
Polarity analysis suggests that glitches have preferential polarization in the directions of the feet of the WTS, the feet of SEIS Leveling System (LVL) and the LSA/tether (Lognonné et al., 2020;Scholz et al., 2020).The fact that a large number of glitches are thermally activated can be an issue for passive seismic applications.Notwithstanding the fact that some glitches have very high amplitude, if their temporal distribution shows some regularity and reproduces at fixed temperature conditions during each Sol, these transients may eventually deteriorate the ACFs and mask interesting signals.To mitigate the risk of contamination by glitches, we apply a glitch-correction algorithm (Scholz et al., 2020) to the raw data, that detects and removes the low-frequency waveform of the most energetic glitches.Due to the variability in their waveforms, some glitches may leave a small imprint in our data, even after correction.For this reason, their impact on our results will be further discussed in section 3.4.2.
In the next sections, we present two methods to retrieve empirically the seismic response below InSight from the pre-processed continuous data.The first method is based on the well established identity between the temporal correlation function of a diffuse wavefield and the Green's function of an elastic medium (Lobkis & Weaver, 2001;Campillo & Paul, 2003).The second method proceeds in the frequency domain and exploits the Wiener-Khintchine relation between the power spectral density of random signals and their autocorrelation function (Yaglom, 2004).rious phases related to environmental noise, whose excitation depends strongly on the daily variation of wind and temperature (Scholz et al., 2020).

Computation of Autocorrelation
A band-pass filter is applied to each 1-hour trace.We subsequently subdivide each trace into windows of 60 seconds duration with 70% overlap.To mitigate the impact of energetic transients, a 1-bit normalization is applied to enforce stationarity of the signal amplitude, thereby improving the signal-to-noise ratio (SNR) (Bensen et al., 2007;Ito & Shiomi, 2012).Following De Plaen et al. ( 2016), we do not apply any spectral whitening.Depending on the bandwidth and the component involved, we also apply several notch filters to remove the lander modes resonances.The low-frequency band 0.4-1 Hz does not appear to contain lander modes.The 1-2Hz and 1-3Hz bands contain a lander mode around 1.6Hz but this one is mainly polarized on the horizontal components.As a consequence, a notch filter centered at 1.6 Hz, which corresponds to the frequency of the lander mode averaged over one Sol, is applied to the East and North components.The 3-6Hz band contains two lander modes at 3.3Hz and 4Hz visible on the three components.
Two notch filters are consequently applied to the Z, N and E traces during the pre-processing.
We used a second-order infinite impulse response notch filter with a quality factor Q = 30.The width of the rejected band at -3dB is approximately one thirtieth of the center frequency.Finally, we compute the full normalized autocorrelations (ZZ, NN and EE) for each of these 60 second-long traces, and we stack them to obtain the ACF for the given LMST and the given Sol.For the computation of ACFs from seismic events waveforms, the processing is almost identical.Different seismic events are simply considered as different Sols and the first subdivision into LMST is not applied.Note nevertheless that a large number of seismic events are recorded during the evening.
In order to check the stability of the various phases observed in the ACFs, we use the definition of the Signal-to-Noise Ratio (SNR) given by Clarke et al. (2011) based on a method of Larose et al. (2007).This SNR is a function of N , the number of realizations, and t, the correlation time-lag.In our case N represent either the number of Sols or the number of seismic events.The SNR is given by: With s(N, t) = ACF (t) +iH( ACF (t) ) and σ(N, t) = Here with i the index labelling the Sol (or the event).
We denote s(N, t) the envelope of the stacked autocorrelation function ACF (t) , H the Hilbert transform and σ(N, t) the amplitude of the residual fluctuations normalized by the number of realizations minus one to avoid biasing.We smooth the time-dependent SNR using a moving time window with the following frequency-dependent duration: 2.5s in the 0.4-1Hz band, 0.5s in the 1-2Hz and 1-3Hz bands and 0.25s in the 3-6Hz and 4.5-7Hz bands.
We show in Figure 4.A the SNR estimated from 149 ambient vibration ACFs derived from data recorded from Sol 222 to Sol 399.We remark that the phases that are visible on individual correlograms as well as on the stacked waveform between 5s and 8s, 10s and 13s and around 21s also correspond to peaks of the SNR.In Figure 4.B, we show the SNR for each hour of the day expressed in LMST.We see that the most energetic arrivals (5-8s, 10-13s and 21s) are visible and stable during nighttime (i.e. from 17:00 LMST to 06:00 LMST) and are particularly clear during the evening (from 17:00 LMST to 23:00 LMST).This period of the day corresponds to the lowest atmospheric activity on Mars.By contrast, when atmospheric noise is particularly strong during daytime (from 06:00 LMST to 17:00 LMST), we cannot discern any clear seismic phase in the ACF.This observation suggests that the signals recorded by SEIS during daytime are not seismic waves but rather ground deformation induced by atmospheric forcing.

Reflectivity via Power Spectral Density (PSD) estimate
By the Wiener-Khintchine theorem, it is known that the PSD of a stationary random signal contains the same information as its ACF, the two being related by a Fourier transform.This property suggests an alternative method to process our data set.First, we compute the average power spectral density of each 1-hour long data segment using the Welch method (Welch, 1967) with sub-windows of 60 seconds duration and 70% overlap.Note that no 1-bit normalization is applied to the original dataset, in sharp contrast with the method described in the previous section.The application of Welch's approach to non-normalized data has been recommended in the framework of ambient seismic vibration processing by Seats et al. (2012).All the PSD estimated between 18:00 and 23:00 LMST on the Z-component are shown in the 0.5-4Hz frequency band in Figure 5.A.At the top of Figure 5.A we show the PSD averaged over all Sols.In Figure 5.A, we observe the broad 2.4Hz resonance first reported by Giardini et al. (2020).Superposed on the broad resonance, we observe characteristic oscillations of the PSD on smaller frequency scales.In order to improve the SNR ratio around the 2.4Hz resonance, we have normalized the PSDs by their maximum in the 1-3Hz frequency band prior to the stack.This procedure downweights data from Sols with anomalously high energies (as around Sol 320 in Figure 5.A) and therefore improves the ensemble averaging.Furthermore, although it may not be apparent on the individual PSDs, there may remain a small imprint of the tick-noise, even after application of the denoising algorithm.To mitigate this issue, we simply replace the value of the PSD at 1Hz (and higher harmonics) by the mean value of the PSD at the two nearest samples.The averaged, normalized PSD is shown in In the following we refer to this method as the "Welch Method".
The deconvolution stage is often used in seismic interferometry to mitigate the effect of the source time function of ambient vibrations, which creates strong side lobes near 0 lag-time (Ruigrok et al., 2010;Casas et al., 2020).In this operation, the band-width of the smoothing window used to estimate the PSD of the source function is a key parameter.In particular, we remark that both the amplitude and the arrival time of the wave packets in the first 4s depend on the smoothing parameter.As a consequence, we will refrain from interpreting any arrival reconstructed by the Welch method at lag time smaller than 4s.

ACFs and SNR analysis
The autocorrelation functions (ACFs) of the ambient vibrations have been computed in four different frequency bands: 0.4-1Hz, 1-2Hz, 1-3Hz and 3-6Hz.The SNR and power spectra of the resulting ACFs are shown in Figure 6  We recall that the evening is known to be the period of the Sol where the atmospheric disturbances are the lowest.Hence, it is only during this time of the Martian day that seismic ambient vibrations may be expected to predominate.The second peculiarity of the 1-3Hz frequency band is the presence of the 2.4Hz resonance.We see on the power spectra of Figure 6.F that this resonance dominates the frequency content of the ACF.
Our interpretation is that the 2.4Hz resonance amplifies the ambient vibration at the In-Sight's landing site, thereby allowing a better reconstruction of the zero-offset reflection response.In the following, we therefore focus our analysis on the results obtained in the 1-3Hz frequency band.
Our main results are shown in Figure 7 where we represent all the ACF waveforms computed from ambient vibrations from Sol 222 to Sol 399 for the Z, E and N components (upper panels), as well as the Z, E and N ACFs waveforms obtained from 124 seismic events (HF, VF and 2.4Hz events, lower panels).In addition, Figures 8, 9 When the envelope of the wave packet presents a clear maximum, common to the three methods, we mark its arrival time with a vertical dark grey line.
The excellent agreement between the three correlograms is worth noting.The close similarity between the waveforms derived from ACFs of seismic events and ambient vibrations is particularly compelling.Furthermore, the agreement between these two ACFs and the results from the Welch method, which proceeds in the frequency domain, con-

Potential contamination by Glitches and Donks
We have seen in section 3.1.2that only the low frequency part of a given glitch is removed by the detection-correction algorithm.This implies that the temporal distribution of glitches can still have an impact on the correlation analysis.In the same way, the distribution of donks is susceptible to contaminate the ACFs.These possible issues are now critically examined.In Figure 11 Having detected clear arrivals in the ACFs, that are furthermore consistent between two different datasets, we now explore the possibility that these phases correspond to the reflection of body waves on deep interfaces under InSight.Other interpretations based on scattered surface waves are of course possible but are left for future works.

Possible reflectivity structure below InSight
We begin our analysis with the vertical ACFs which show high S/N ratio (> 4) up to lapse time of 30s.Following the basic principles of seismic interferometry, we interpret the identified phases in terms of P-wave reflectivity.The most energetic arrival (besides the side lobe at time t = 0) seen on both ambient vibration and event ACFs arrives at a two-way traveltime of 5.6s.This arrival is most directly interpreted as a reflection from an upper-crustal interface that was also detected by Receiver Function analysis.As reported by Lognonné et al. (2020), the plausible depth of this interface is 9.6km (±1.8km), corresponding to an arrival time in the range 3.5s-8.2s.A series of later arrivals is detected at lapse-time 10.6s, 12.6s, 21.s, whose amplitudes decay gradually.As an illustration, the amplitude ratio between the 5.6s and 21s phases is typically in the range 1.5-2.This makes it rather unlikely that any of the late arrivals corresponds to a multiple reflection.Assuming the same range of P-wavespeeds in the lower crust as Lognonné et al. ( 2020) (reported in Table 4), our observations suggest the presence of several deep interfaces at approximate depths of 21.6km ±4.8km, 26.4km ±6km and 46.5±11.1km,corresponding respectively to arrivals at 10.6s, 12.6s and 21.0s.The uncertainty on the depth of the interfaces is directly inherited from our assumption on the plausible range of P-wave velocities in the crust.This analysis of the ZZ ACFs is summarized in Table 4.
We now examine the horizontal ACFs and interpret them in terms of shear wave reflectivity.In doing so, we try to identify arrivals that possibly correspond to reflections from the interfaces deduced from the vertical ACF interpretation.This task is complicated by the fact that the arrival times of phases seen on the EE and NN ACFs do not necessarily match.Comparing the Figures 9 and 10 we observe that the NN component presents more similarity up to 40s lag-time between the three methods than the EE component.Since the agreement between ambient vibration and event ACFs is by far better on the NN component than on the EE component, we focus our interpretation on the former.We observe a broad high-amplitude arrival between 4s and 8s lapse-time.This arrival corresponds to lapse-time that are too short to be unambiguously interpreted on the ZZ ACFs.The dispersed nature of the signal may be the signature of scattering from the surface regolith which hampers the identification of an interface.A series of later and lower amplitude arrivals can be discerned simultaneously on the event and ambient vibration NN ACFs at lapse-times of 11.9s, 14.4s, 16.5s and 22.4s.The simplest interpretation of these arrivals, compatible with the P-wave reflectivity is the following one.The 11.9s phase could correspond to the reflection of S waves from the upper crustal interface detected on the ZZ ACF.This allows us to determine the P-to-S wavespeed ratio 1 which in turn puts the first interface at a depth of 11.6km ±0.9km.The reduction in uncertainty with respect to the ZZ analysis comes from a far better constraint of the S-wave velocities by Lognonné et al. (2020) (reported in Table 5) than the P-wave velocities.The next phase that we identify arrives at 22.4s and corresponds to the 10.6s reflection seen on the ZZ ACF.Taking into account the range of shear wave velocities deduced from RF by Lognonné et al. (2020) and shown in Table 5, this implies that the second interface could be located at 25.2km ±2.5km depth.We summarize this analysis in Table 5.This interpretation is not entirely satisfactory, as we observe two arrivals at respectively 14.4s and 16.5s that are difficult to reconcile with the P-reflectivity profile.The two uninterpreted arrivals on the NN ACF correspond to phases arriving between 5.6s and 10.5s lapse-time for P-waves.In this time window, the ZZ ACF shows a broad wave packet where individual arrivals are difficult to discern.Interestingly, the 14.4s arrival is visible on both the EE and NN ACF, which indeed suggests the presence of an interface between 11.6km and 25.2km depth.

Discussion
We now critically examine our interpretation as well as those proposed in the recent literature.The work by Deng and Levander (2020) is particularly relevant to ours.
These authors used vertical ACFs at long and short period to propose a stratified view of Mars interior.We will focus the discussion on the short period band that is common to both studies.Deng and Levander (2020) reported the observations of two arrivals at 11.5s and 21.s.The most notable difference with our findings is that instead of a single arrival at 11.5s, we observe two arrivals at 10.6s and 12.6s.We therefore see the interest in considering both ambient vibrations and coda waves in the processing.
The plausibility that the reconstructed phases correspond to deep reflections depends crucially on the crustal attenuation below InSight.As reported in Lognonné et al. (2020), the scattering attenuation as quantified by the seismic diffusivity D = 90km 2 /s is rather moderate.This diffusivity corresponds to a scattering attenuation length of about 90kms (Lognonné et al., 2020) so that ballistic waves propagating two-way through the crust would see their amplitude reduced by a factor of 3 due to scattering attenuation.
The absorption quality factor Q has also been estimated to be in the range 1250-1325 for S waves at 2.5 Hz and is probably much higher for P waves.Thus absorption likely affects negligibly the ballistic wave amplitudes.All these estimates are of course subject to large uncertainties but they do not rule out our interpretation of ACFs as body waves reflections.
Although scattering is moderate in the Martian crust, we have observational evidence that the wavefield of high-frequency Martian events is diffuse (Lognonné et al., 2020).Since the minimal hypocentral distance of the events used in our study is larger than 500km (Giardini et al., 2020), the diffuse character is acquired by long range propagation through an inhomogeneous crust rather than by strong local scattering (Lognonné et al., 2020).Coda records are more likely to possess the diffusive character required for the application of seismic interferometry.Therefore, by contrast with Deng and Levander (2020) we would rather interpret the phase detected at 21s as a PmP rather than as an SmS.We also note that the horizontal ACFs, which are supposed to be more faithful to the S reflectivity profile, lack a clear reflection associated to the 21s phase seen on the ZZ ACF.Therefore, we believe that the depth of the Moho under InSight is still to be confirmed by additional studies.To carry out this task, a joint interpretation of the vertical and horizontal ACFs would be necessary.This will require more advanced modeling approaches such as full-waveform inversion.Even more importantly, additional constraints could come from the analysis of converted phases as in RF analysis, which in addition provide constraints on the absolute velocities in the crust.

CONCLUSION
The ground velocity records of the SEIS instrument have been analyzed with seismic interferometry methods.The stability analysis of the autocorrelation functions of SEIS components demonstrates that the ambient seismic vibration is most reliably observed in a specific frequency band (2.4Hz resonance) and only when the environmental noise is lowest (17:00 to 23:00 LMST from Sol 222 to Sol 399).Based on the Power Spectral Density (PSD) of ambient vibrations, we show that the 2.4Hz resonance has an oscillating structure which is perfectly stable with time and which is directly related to the Empirical Green's functions reconstructed in the 1-3Hz band by autocorrelation analysis.The good agreement between the autocorrelation functions computed on ambient seismic vibration in the 1-3Hz range and the autocorrelation functions computed on the Marsquake waveforms is consistent with the interpretation of the 2.4Hz resonance as a local ground amplification due to the shallow structure beneath the InSight's landing site.
The autocorrelation functions present seismic energy arrivals in the 4 to 30 seconds lagtime range that are validated by Signal to Noise Ratio (SNR) analysis and inter-comparisons between results obtained from ambient vibration and seismic event records.We report the possible detection of vertically propagating P-waves reflected on internal discontinuities with two-way travel times of 5.6s, 10.6s, 12.6s and 21s.A clear identification of the corresponding S-waves reflections is more speculative due to the poor agreement between the EE and NN components.Nevertheless, based on the more reliable component (NN) we suggest that the 11.9s and 22.4s arrivals could be the S-waves reflections corresponding to the P-waves reflections at 5.6s and 10.6s.Two internal structure models deduced from these travel times are presented, but they must be further constrained by other seismic analysis such as receiver functions, in order to obtain a reliable internal structure model.

HF events
Name Quality    -21-manuscript submitted to JGR: Planets .A a spectrogram of raw records of the VBB-V axis at 20 samples-per-second (sps) between Sol 183 and Sol 190.The time windows with highest energy correspond to sunlight periods.The clear difference between daytime-and nighttimerecords is due to atmospheric processes.As described in Lognonné et al. (2020), atmospheric noise entails elastic deformations induced by pressure effects on the ground, tilt of the lander under wind, and lander vibrations under wind.The changes in the speed and turbulent flow of the wind are the main drivers of the SEIS background signal.Next, we remark that certain narrow frequency bands are noticeably more energetic than others.Careful examination reveals the excitation of resonances whose central frequency exhibits a daily modulation.Some of these resonances are clustered in the frequency domain between 3.2Hz and 4.5Hz (as delimited by the red rectangle in Fig. 1.A).This particular frequency band contains a series of modes of resonance of the InSight lander which are continuously excited by the wind.The complexity of the modal distribution stems manuscript submitted to JGR: Planets et al. (2020) as a local ground resonance.Nevertheless, the exact way this resonance is generated and excited remains unknown.
Functions (ACFs) and Signal-to-Noise Ratio (SNR) analysis To compute the velocity ACF of Martian ambient vibrations, we follow the workflow described by Bensen et al. (2007) with minor modifications.After tick noise and glitch removal, the instrument response is removed and the traces are rotated onto the local geographical coordinate system Z, N, E, with Z the upward vertical, N the horizontal North component and E the horizontal East component.The signal is subsequently cut into segments of one Sol duration which are processed independently.Due to the highly non-stationary character of both the ambient Martian vibrations and the perturbations caused by the environment, we further process independently each hour (in Local Mean Solar Time) of the continuous record.This hourly processing allows us to check the stability of the ACFs and to remove segments that are potentially contaminated by spu- Figure 5.B.Using both numerical simulations and real data examples,Oren and Nowack (2017) have suggested that the rapid oscillations of the PSD may contain information on the reflection of propagating waves from deep interfaces.Therefore, to extract the reflectivity structure under InSight, it is natural to deconvolve the PSD from the smooth 2.4Hz resonance.In the example shown in Figure5.C, this operation is performed by dividing the original PSD by its smoothed version using a 0.5Hz window.In order to remove contamination from unwanted spectral oscillations caused by Lander modes, we also normalize the PSD to one outside the frequency band of interest (Figure5.D).The final step of the process consists in taking the inverse Fourier transform of the deconvolved PSD to reveal the arrival times of the possible seismic phases at the origin of the spectral oscillations.
for the ZZ component.The high values of SNR (in red) near 0 lag-time are caused by the noise source time function since no source deconvolution has been applied during the processing.We see in Figure 6 that outside of the time window near 0 lag-time, the SNR is close to zero everywhere except in the 1-3Hz bandwidth between 17:00 LMST and 23:00 LMST (Figure 6.E).
and 10show the average waveforms for the three approaches outlined above: seismic events, Welch method, and ambient vibration for the ZZ, EE and NN components respectively.The ACFs referring to "Seismic Events" have been obtained by stacking linearly the individual ACFs of all HF events, VF events and 2.4Hz events (124 diffuse seismic events in total) shown in Figure7.D.E.F.The ACFs referring to "ambient vibrations" have been obtained by stacking linearly the hourly ACFs during the evening period (17:00 to 23:00 LMST) of each Sol and by subsequently averaging the results from Sol 222 to Sol 399 shown in Figure7.A.B.C.Note that because source deconvolution has been applied in the Welch method, the waveform amplitudes shown in panel B are not directly comparable with the one displayed in panels A and C. In Figure8, 9 and 10 we highlight in gray the wave packets that are simultaneously visible on the three different types of ACFs.
firms that the arrival time of the different wavepackets is indeed encoded in the periodical oscillations of the PSD around the 2.4Hz resonance.Again, this should come as no surprise by the Wiener-Khintchine theorem.Examining Figures 7.A.B.C, we also notice that the individual ambient vibration ACFs computed over a single Sol already contain all the arrivals that show up in the complete stack.By contrast, the individual ACFs computed from seismic events waveforms (Figure 7.D.E.F) exhibit much more variability.Quite remarkably, the average ACF nevertheless converges toward the same waveform as ob-tained from the ambient seismic vibration.Interestingly, Hillers and Campillo (2016) observed the same variability of cross-correlation waveforms derived from earthquake coda waves on Earth.This comforts us in our interpretation of the arrivals visible in the ACFs as true propagating seismic signals.
.A, we show the statistics of the lapse-time between pairs of glitches detected in the 20 sps VBB data.A clear pattern can be distinguished on the plots, where each glitch is represented as a point in the (LMST, Sol number) plane.In particular, a high density of glitches can be observed along certain curves suggesting that glitches are driven by environmental forcing.Indeed,Scholz et al. (2020) have shown that the curves correspond to particular values of the temperature, which varies seasonally.The regularity and repeatability observed in the glitch distribution leads us to take into consideration the time delay between two consecutive glitches.This parameter is key to unravel a potential contamination of ACFs.In Figure11.A, we show the histogram of the time delays between two consecutive glitches on the V component at different LMST.We focus our analysis on the evening part of each Sol because it is the period during which we obtain the best results.We see on the histogram that the smallest delay is approximately 30 seconds and that the mean delay between two consecutive glitches is around 200 seconds.These values are too large to explain the arrivals visible on the ACFs.Moreover we see that the distribution of the delays depends on the LMST.This observation is not compatible with the stability of the arrivals over the nighttime windows.We show in Figure11.D the temporal distribution of the donks detected between 17:00 and 23:00 LMST from Sol 180 to Sol 261.This period is delimited by the red rectangle in Figure11.B.As the 100 sps SP data are not always available, we perform the detection of donks on a composite SP channel called ESTASP (Energy Short Term Average -SP) which is available in the continuous data stream.The output of this channel is defined as the root mean square of the raw vertical SP components (SP1), filtered in the 12-14Hz bandwidth and averaged over 1 second(Lognonné et al., 2019).The peaks visible on the ESTASP channel time series shown in Figure11.E are markers of donks.Application of a simple STA/LTA (1s/25s) criterion thus permits the automatic detection of the vast majority of donks from ESTASP records.In Figure11.D we may observe three time windows with a particularly high density of donks at the beginning of the evening.Each time window opens at a particular Sol and extends gradually Sol after Sol.We hypothesize that the appearance of a high number of donks is related to the seasonal reactivation of lander cracks.We see on the histogram of Figure11.C that the typical delay between two donks is much smaller than the delay between glitches.A broad maximum in the distribution of delays is visible around 10 seconds between 17:00 and 18:00 LMST.Such a peak could cast doubts on the physical origin of the arrival seen around 10.5s in the ZZ ACFs.Nevertheless, there are a number of arguments that favor the interpretation of this arrival in terms of wave propagation.(1) The phases seen in the ACF's are stable during the whole evening, while donks activity strongly decreases after 18:00 LMST.(2) Donks leak very little energy into the 2-3Hz bandwidth.(3) All the arrivals seen on ambient vibration ACFs are also present in the events ACFS.The latter are tens to thousands of times more energetic than the ambient vibration, making donk contamination rather unlikely.

Figure 1 .
Figure 1.A) Spectrogram of a 6 Sols-long record of raw 20sps VBB data (V component in count).The spectrogram is computed using a moving time window of 15 minutes without overlap.The vertical bands of high amplitude correspond to daytime windows, i.e. between 5:00 and 16:00 LMST (Local Mean Solar Time).B) Mean power spectrum of the same data (in dB).The red rectangles delimit specific frequency ranges where lander modes are clustered.The black rectangles indicate the broad resonance around 2.4Hz.

Figure 2 .
Figure 2. A) Spectrogram of the raw VBB-V component at 20sps at the beginning of the evening of Sol 348.The spectrogram is computed using a moving time window of 100 seconds with 20% overlap.Each high energy peak in the low frequency band (< 0.1Hz) is the signature of a glitch.Glitches are also clearly visible in the time-series (red line) superposed on the spectrogram.The black rectangle delimits the time window of the glitch presented in B).In inset, the whole VBB-V record is shown for Sol 348.The red rectangle delimits the time window used for the spectrogram computation.B) Typical raw glitch waveform.C) Spectrogram of the raw SP2 (horizontal) component at 100sps during the evening of Sol 348.The spectrogram is computed using a moving time window of 5 seconds with 20% overlap.Each high-energy peak in the high frequency band is the signature of a donk.Donks are also clearly visible in the time-series (red line) superposed onto the spectrogram.The black rectangle delimits the time window of the donk presented in sub-Figure D).The inset in sub-Figure C shows the complete SP2 record for Sol 348.The red rectangle delimits the time window used for the spectrogram computation.D) Typical raw donk waveform.

Figure 3 .
Figure 3. A) Estimated waveforms of the tick noise on the U, V and W channels of SEIS-VBB at 20sps.B) Spectrograms of the U, V and W components of SEIS-VBB before (raw) and after (corrected) tick-noise removal.The spectrograms have a temporal resolution of 1 hour and are represented in the 0.999Hz-1.001Hzfrequency band.C) Raw VBB-V record between 19:00 and 21:00 LMST (black line) on Sol 183.Zooms into the time-window delimited by the red rectangle are provided in insets to highlight the effect of the correction.D) Amplitude Spectral Density (ASD) of the Raw VBB-V record shown in C) before and after tick-noise removal.

Figure 4 .
Figure 4. A) From top to bottom : ZZ-Correlograms computed from Sol 222 to Sol 399 between 19:00 and 20:00 LMST in the 1-3Hz frequency band; Resulting stacked ACF; SNR as a function of the correlation lag-time, following the method of Clarke et al. (2011).B) SNR (see colorbar) as a function of LMST (Local Mean Solar Time) and correlation lag-time for the ZZ component filtered between 1Hz and 3Hz.

Figure 5 .Figure 6 .
Figure 5. A) Spectrogram showing the power spectral density (PSD) in the 0.5-4Hz frequency band during the evening period (18:00 -23:00 LMST) from Sol 222 to Sol 399.The spectrogram is computed using a moving time window of 1 minute with 70% overlap.The averaged PSD is shown at the top of the spectrograms.B) Averaged and normalized PSD estimated between 18:00 and 23:00 LMST and its smoothed version.C) Deconvolved PSD.D) Deconvolved PSD after frequency band selection.E) Inverse Fourier Transform of the PSD shown in D).

Figure 7 .Figure 8 .
Figure 7. A) ZZ, B) EE and C) NN autocorrelation functions of ambient vibrations computed between 17:00 and 23:00 LMST from Sol 222 to Sol 399 in the 1-3Hz frequency band.D) ZZ, E) EE and F) NN autocorrelation functions of HF, VF and 2.4Hz events in the 1-3Hz frequency band, ordered by event number.

Figure 9 .Figure 10 .
Figure9.Same as Figure8for EE ACFs in the 1-3Hz frequency band.Note that in the Welch method, an additional source deconvolution is applied.

Table 1 .
Summary of acronyms used in the text.

Table 2 .
Azimuth (Az)and Dip (angular deviation from horizontal) defining the spatial orientation of the 6 sensors composing the SEIS instrument.

Table 4 .
Extension of the crustal model by Lognonné et al. (2020) derived from the main arrivals observed in the ZZ ACF.

Table 5 .
Extension of the crustal model by Lognonné et al. (2020) based on a combined interpretation of ZZ and NN ACFs.

Table A1 .
List of seismic events used in this study referring to the catalog V3 by the MQS(InSight Marsquake Service, 2020) S0194c B provement of seismic records, and characterization of long period atmospheric waves from ground displacements.Journal of Geophysical Research: Planets,