Detection, Analysis, and Removal of Glitches From InSight's Seismic Data From Mars

The instrument package SEIS (Seismic Experiment for Internal Structure) with the three very broadband and three short‐period seismic sensors is installed on the surface on Mars as part of NASA's InSight Discovery mission. When compared to terrestrial installations, SEIS is deployed in a very harsh wind and temperature environment that leads to inevitable degradation of the quality of the recorded data. One ubiquitous artifact in the raw data is an abundance of transient one‐sided pulses often accompanied by high‐frequency spikes. These pulses, which we term “glitches”, can be modeled as the response of the instrument to a step in acceleration, while the spikes can be modeled as the response to a simultaneous step in displacement. We attribute the glitches primarily to SEIS‐internal stress relaxations caused by the large temperature variations to which the instrument is exposed during a Martian day. Only a small fraction of glitches correspond to a motion of the SEIS package as a whole caused by minuscule tilts of either the instrument or the ground. In this study, we focus on the analysis of the glitch+spike phenomenon and present how these signals can be automatically detected and removed from SEIS's raw data. As glitches affect many standard seismological analysis methods such as receiver functions, spectral decomposition and source inversions, we anticipate that studies of the Martian seismicity as well as studies of Mars' internal structure should benefit from deglitched seismic data.


Introduction
InSight (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport) landed successfully on Mars on 26 November 2018 (Sol 0, a sol is a Martian day with around 24 hr 40 min). Since 9 February 2019 (Sol 73), InSight's main scientific instrument SEIS (Seismic Experiment for Internal Structure) is recording seismic data in its operational configuration . The SEIS package (Lognonné et al., 2019), which is located at 135.6234°E and 4.5023°N in Elysium Planitia  and whose network and station code for the scientific data is XB.ELYSE, consists of 2 three-component seismometers; one being very broadband (VBB) with a corner period of 16 s, and one being short-period (SP) with a corner period of 35 s. The noise floor of the two instruments is equivalent only above 5 Hz while it is about ∼30 dB lower for the VBB at frequencies below 0.1 Hz. It is this frequency dependence of the seismometers' self-noise that determines their names as established for the InSight project (Lognonné et al., 2019), even though the naming convention does not follow terrestrial standards (e.g., Ahern et al., 2012). Due to their different noise floors, the VBB is the main instrument to detect distant marsquakes, while the SP is used to cover the frequency range of ∼5-50 Hz for more detailed analysis of regional events and lander-induced signals. Both seismometers have nonorthogonal sensor orientations (Figures 1a and 1c). To date, all six seismic components as well as the acquisition system have functioned nominally, exceeded mission requirements, and delivered unprecedented seismic data from the surface of Mars (InSight Mars SEIS Data Service, 2019). In addition to seismic signals of natural and artificial origins, that is, marsquakes Giardini et al., 2020;Clinton et al., 2020; InSight Marsquake Service, 2020 for event catalog) and records from the HP 3 instrument hammering sessions (Spohn et al., 2018), respectively, these data show a variety of nonseismic signals . Among the most prominent and abundant types of these nonseismic signals are what we termed a "glitch". Glitches influence many of the standard seismological methods such as receiver functions, polarization analyses, source inversions, and spectral decomposition. On Earth, data influenced by such disturbances are often discarded especially when coinciding with earthquake phase arrivals (e.g., Zahradnik & Plesinger, 2005). This obviously represents no valid option for the seismic data returned from Mars, hence the correct treatment of the glitches is of upmost importance for the scientific analyses. Within the InSight team, four groups therefore independently investigated the glitch phenomenon to ensure redundancy of the effort, reproducibility of the outcomes, and convergence toward a complete understanding. To some extent these efforts also reflect the scientific objectives across the different groups as well as institutional responsibilities within the InSight project. The present study is the current culmination of this teamwork and focuses on the detection, analysis, and removal of glitches. It extends Supplement V of Lognonné et al. (2020).

Glitches
In the literature (e.g., Iwan et al., 1985;Zahradnik & Plesinger, 2005;Vacka et al., 2015) the phenomenon we are investigating here is sometimes referred to as "long-period disturbances", "acceleration offsets", or even "mice", all generally describing the same type of data disturbance. Throughout the present publication, however, we choose to apply the term "glitch" to these disturbances as it has been established as such since their first observations in InSight's seismic data and hence been communicated so to a wider audience on various occasions. While we are aware that the word glitch is typically associated to more general data artifacts and alike, we indeed use it here to refer to specific, clearly defined disturbances in the data. A glitch (Figures 1b  and 1d), thus, is a particular type of transient instrumental self-noise in the raw time series data that appears as a high-amplitude, one-sided pulse with a duration controlled by the seismometer's transfer function. For the VBB sensors, which have 76% of critical damping, glitches have a fast rise time followed by an exponential decay with a ∼9% overshoot before almost returning to the baseline after ∼25 s. For the SP sensors, that are overdamped with 110% of critical damping, glitches have a similar rise time followed by a decay before almost returning to the baseline after ∼50 s. Glitches may also occur before a previous glitch has sufficiently decayed. The highest order of such "polyglitches" we observe to date is four. Glitches (and polyglitches) can occur on all three VBB and all three SP sensors simultaneously, but for VBB about one-third occurs on only one component while for SP about two-thirds occur on only one component. Glitches occur at all times of the sol but are observed more frequently during the quiet parts in the early evening and night (Figure 2a). This is due to the decreased seismic noise level driven by diurnal wind and pressure variations. The largest glitches reach amplitudes of 10 −7 m/s and more. We observe a few of these per sol, while for amplitudes of 10 −8 m/s we can observe already hundreds per sol. Especially in the early evening, when the wind and pressure variations have calmed down, we observe a period with many consecutive glitches mostly of lower amplitude (Figure 2b). Certain types of glitches can furthermore repeat over many consecutive sols at the same local time, thus indicating a driving process behind their generation. In the frequency domain, glitches range from lowest frequencies up to almost 1 Hz, thus influencing analyses of seismic records especially for longer periods.
Premission tests of the SEIS instrument showed no indication of glitches in the data records. For example, in the heat chamber, where temperature cycles as expected on the surface of Mars were generated, the seismic noise was too elevated to identify glitches as such, if they occurred. At the Black Forest Observatory in Germany, among the quietest seismic sites for long-period noise in the world, the thermal environment was very stable and the tests did neither suggest the presence of glitches. Glitches and especially their abundance were thus not expected prior to the InSight mission.  Earth and Space Science SCHOLZ ET AL.

Glitch Spikes
Many glitches, furthermore, show a high-frequency signal at their glitch beginning that lasts around 40 samples regardless of the data sampling frequency. We refer to these initial oscillations as "glitch spikes". These spikes occur simultaneously with the glitch onset for both VBB and SP (Figures 1b and 1d). In the frequency domain, glitch spikes can range from ∼1 Hz up to the data's Nyquist frequency. They do not represent artifacts caused by the on-board analog or digital electronics but instead-as we demonstrate-are closely linked to glitches and their causes.

Glitch Detection
To automatically detect glitches on SEIS's VBB and SP raw data, several groups (MPS, ISAE, UCLA, and IPGP) independently developed algorithms in the Python and MATLAB programming languages. The group acronyms stand for the affiliation of the group's leading analyst, that is, Max-Planck-Institute for Solar System Research (MPS), Institut Supérieur de l'Aéronautique et de l'Espace SUPAERO (ISAE), Department of Earth, Planetary, and Space Sciences, University of California Los Angeles (UCLA), and Université de Paris, Institut de physique du Globe de Paris (IPGP). We describe each approach in the following. The common detection idea and working hypothesis of this study is that glitches in the raw data represent steps in acceleration convolved with the seismometer's instrument response, while spikes represent steps in displacement convolved with the seismometer's instrument response. The lists of detected glitches in 2019 can be found in the Supporting Information S1.

Glitch Detection by Instrument Response Deconvolution (MPS)
This detection algorithm, implemented in Python (Rossum, 1995) and ObsPy (Beyreuther et al., 2010;Krischer et al., 2015), performs the following processing steps on a given period of three-component seismic data (components U, V, W): (i) decimate the data to two samples per second (sps), allowing all data per seismometer to be run with the same parameters and enabling faster computations; (ii) deconvolve the Color-coded symbols correspond to glitches for the different groups that are not common to all. Those common to subgroups are plotted on top of each other and so the last plotted is shown. (b) Expanded section showing that as the threshold for declaring a glitch is lowered, either in terms of signal-to-noise or correlation with the template, the results differ and a few, small candidate glitches have been missed.

10.1029/2020EA001317
Earth and Space Science SCHOLZ ET AL.
instrument response on each component and convert to acceleration; (ii) band-pass filter the acceleration data (e.g., 10-1,000 s), so the steps in acceleration emerge more clearly; (iv) calculate the time derivative of the filtered acceleration data so the acceleration steps become impulse-like signals; and (v) on this time-derivative, trigger glitches based on a constant threshold. To avoid triggering on subsequent samples also exceeding the threshold but belonging to the same glitch, we introduce a window length in which no further glitch can be triggered. This parameter can be thought of as glitch minimum length. We note that this parameter is smaller than the typical glitch length for VBB and SP, allowing our detection algorithm to detect polyglitches.
A glitch simultaneously occurring on multiple components is detected on each affected component but the respective start times may slightly differ. However, after modeling of the full glitch waveform (section 4) we can retrospectively establish that such glitches occur at the same time to within milliseconds. This holds true for all multicomponent glitches observed to date on either VBB or SP, also for data with the highest available sampling frequency of 100 Hz. Therefore, we declare as glitch start time the earliest time detected across the UVW components. The list of unified glitch starts contains still many false-positive triggers caused by nonglitches with a steep enough acceleration change to be triggered. This is because we choose to apply a constant threshold to the time derivative of the filtered acceleration, rather than a threshold based on the current seismic noise level that undergoes strong diurnal changes (amplitudes varying by a factor of 100 and more) dominated by meteorological influences (e.g., Banfield et al., 2020;Lognonné et al., 2020). To circumvent, we rotate the gain-corrected UVW raw data of the glitch windows into the geographical reference frame (ZNE components) and perform a 3-D principle component analysis (e.g., Scholz et al., 2017). Theoretically, a glitch is linearly polarized as the associated vector of acceleration change is not varying, however slightly altered only by seismic noise. Indeed, most glitches exhibit a high linear polarization >0.9 which we use to discriminate against other triggered signals. The polarization analysis further allows to obtain the apparent glitch azimuth and incidence angles which we use to associate glitches with particular glitch sources (section 3). Visual inspection reveals the detected glitch onsets are usually accurate to within ±1 s with respect to the theoretical glitch onsets (Figures 1b and 1d).

Glitch Detection by Cross Correlation With Impulse Response Function (ISAE)
The principle of this MATLAB-implemented detection algorithm is cross correlation. It performs the following processing steps on a given period of three-component raw seismic data (components U, V, W): (i) a synthetic glitch is constructed by convolving the poles and zeros of the transfer function of the VBB and SP sensors with a step in acceleration. To increase the temporal resolution to subsample range, we synthesize several glitches each with a different subsample time shift; (ii) while the frequencies above 2 Hz are filtered, the long period variations of the data are extracted using a low-pass filter with 10 −3 normalized cutoff frequency for VBB and 0.25 × 10 −4 normalized cutoff frequency for SP. These are then subtracted from the signal (and added back at the end), before (iii) the synthetic glitch is cross-correlated with the data. A glitch detection is triggered for the maxima of the cross-correlation function that exceed a threshold a on a given component.
Another step is added to prevent nondetection of glitches or false-positives, depending on the correlation threshold. For that, two thresholds are chosen: threshold a and threshold b, with a ≥ b. The first step presented above is done for each component, with threshold a. Then, for each component, a second cross correlation with threshold b is implemented. For the times of every maximum of cross correlation exceeding threshold b, we come back to the glitches detected on the other components during the first step. If a glitch had indeed been detected at that specific time on another component, a new glitch is declared on the component under study. We can therefore detect small glitches with low signal-to-noise ratio when a strong glitch is detected at the same time on some other component. In addition, in order to be able to detect polyglitches, a second iteration of the detection algorithm is performed after the glitches from the first iteration have been removed from the data.

Hierarchical Glitch Detection (UCLA)
This MATLAB-based method took into account that glitch amplitudes follow a power law distribution with many more very small glitches than larger ones (see Figure SI2-3 in Text S2 in the supporting information). Therefore, the strategy was to remove the largest glitches first and repeat the process on the smaller ones.

10.1029/2020EA001317
Earth and Space Science SCHOLZ ET AL.
In this method, the raw UVW VEL channel data are inspected for glitches and their spikes. The instrument response to a step in acceleration was termed Green's function. The 20 sps data were decimated to 2 sps and each channel was tested for correlation with the response function as follows.
An inverse filter was designed that turned glitches into narrow Gaussians with rise times equal to the glitch so that each glitch represented one peak without the overshoot. This enables detection of multiple close-spaced glitches. An STA/LTA (short time average/long time average) ratio was found using convolution of the data with two boxcar functions separated by more than a glitch window. The absolute value of band-passed data was tested for peaks above the STA/LTA threshold. For the first iteration the STA/LTA was set large to remove the largest glitches. The Green's function was correlated with the data spanning a peak and if the correlation coefficient was above 0.90 the detection was registered. If multiple peaks occurred close together, multiple Green's functions were fit to the data using nonlinear least squares. The data were then cleaned by removing the glitches. The process was then repeated lowering the STA/LTA threshold ¼ 7, and the new glitches removed from the data. For the last iteration the STA/LTA threshold was set to 3, that is, lowered again and the correlation threshold was also lowered to 0.8. This removed many of the small glitches. Our glitch detection is applicable to SEIS's VBB and SP sensors in both low and high gain modes.

Triple-Source-Based Glitch Detection (IPGP)
Implemented in MATLAB, this glitch detection method processes mostly 2 sps continuous data and is therefore focused on long period continuous signals. It first removes the aseismic signals of each raw axis by subtracting the trend and the first 12 integer multiples of one cycle per sol (first 12 "sol-harmonics", that is, 12 frequencies between ∼0.01 and 0.1 mHz). Then the three axes are equalized in digital units by convolving the V and W channels by the convolution ratio of the U/V and U/W transfer functions, in order to correct for the gain and transfer function differences between U, V, and W. Note that this process also transforms an impulse response in time on V and W into an impulse response with the U transfer function. As the inversion is a linear one, the glitch search and deglitching can be done either on the UVW or on the ZNE rotated channels, with practically no differences for the inverted glitches.
The glitch detection is done first by identifying all extrema in the signal and then least square testing them for the occurrence of a glitch using a modeled glitch. To model a glitch, we convolve a step in acceleration not only for one sample (as all other methods) but for three consecutive data samples, that is, we model three glitches within 1 s of data. As we have equalized all components beforehand, we only use the poles and zeros of the U component for this step. Continuity of the signal is forced at the beginning and at the end of the glitch window by Lagrangian multipliers. The signal is then considered a glitch when the variance residual after glitch removal is less than 1-2% of the original data squared energy over a running window of 50 s, starting 5 s before the glitch maximum. To remove the glitch spikes after the glitch removal, a delta impulse is then searched around the glitch onset time and removed if associated with a 50% variance reduction of the signal in a window of width ±3 s. Glitches and spikes amplitudes are inverted on the three axes. We use these amplitudes to calculate dip, azimuth, and amplitudes of the spikes that we use to potentially locate glitch sources (section 6.1). An average of about 170 glitches per sol is found for 1% of variance residual and about 100 glitches per sol for 0.5% of variance residual. For the former case, about 40% are detected on the three components while the other are on single VBB components. As this approach is detecting the glitches through the success of the functions' fit with the data, glitch removal is a subproduct of this method.

Performance of Glitch Detection Algorithms
A 24-hr comparison of our glitch detection algorithms is illustrated in Figure 2. The detection threshold for some methods was set low in order to examine differences in the detections close to the ambient seismic noise levels. For example, ISAE and UCLA used a correlation coefficient threshold of 0.8 which opens the possibility that some of the detections may be noise. Approximately 250 detections were made by UCLA and IPGP, and 140 by MPS and ISAE; however, the latter two detected fewer glitches during the daytime noise. Figure 2a shows the 73 glitches that were common to all four groups, which correspond to those with the largest amplitude. Table 1 shows the number of detected glitches common to pairs of groups. The noncommon glitches are plotted color-coded according to each group. An expanded section (Figure 2b) reveals that the various criteria detect mutually exclusive glitches as the noise level is approached. We note that the Marsquake Service (MQS, Clinton et al., 2018Clinton et al., , 2020 continuously monitors InSight's seismic data to detect and catalog seismic events (InSight Marsquake Service, 2020). As part of their routine they manually seek and annotate glitches with principal focus on time windows of seismic events. Our detection methods generally compare well with these manual annotations both in amount and onsets of glitches, especially for larger ones. For smaller glitches, that is, less than 10 −8 m/s in amplitude, we find that each detection method, if the parameters are chosen sensitively enough, delivers satisfying results with the amount of false detections only slightly increased. However, not each annotated glitch is detected as the noise level is approached and the signal-to-noise ratio hence decreases. Nevertheless, our comparisons show that our algorithms for glitch detection are reliable in most circumstances.

Glitch Analysis
Our working hypothesis is that glitches in SEIS's time series data represent sudden steps in the sensed acceleration convolved with the instrument response of the respective seismometer, either VBB or SP. We can use that assumption to constrain the physical mechanism that led to the glitch. When interpreted as an inertial acceleration of the seismometer frame, a step in acceleration translates to a unlimited, linear change of velocity. This of course becomes quickly nonphysical and can be ruled out because it implied that SEIS would have left its landing location by now. On the other hand, accelerometers like the VBB or SP are also sensitive to changes in gravity. One way this can occur is by tilting the instrument, thus changing the projection of the local gravity vector onto the directions of the sensitive sensor axes. For small tilt angles α, this translates into a first-order effect for the horizontal components (∼sin(α) ≈ α) but only into a second-order effect for the vertical component ( ∼ ½1 − cosðαÞ ≈ α 2 =2). The vector sum of acceleration changes in U, V, and W due to a tilt of the SEIS sensor assembly (including the leveling system) will therefore point in the horizontal direction. This is true for both SP and VBB. Any other direction cannot be explained by a rigid motion of SEIS and must be due to instrumental artifacts.
It is useful to recall the sign convention for accelerometers: A positive output signal corresponds to a positive acceleration of the frame in the sensitivity direction, not the direction in which the proof mass moved. Therefore, if one analyzes the apparent glitch azimuth and incidence angles under consideration of the actual sensor orientations as well as the behavior of these angles over time, one can draw conclusions on possible glitch origins. The analysis of apparent glitch polarizations is therefore our method of choice.
The determination of the apparent glitch azimuth and incidence angles is implemented in our glitch detection algorithm (section 2.1) and based on a 3-D principle component analysis. To resolve the 180°ambiguity of the azimuths inherent to that method, we used the fact that glitches have a clear one-sided pulse ( Figures  1b and 1d); a glitch of positive polarity on the N component is associated with a step in acceleration acting in this direction, its respective azimuth is therefore ≈0°(assuming there is no glitch on the E component). The same consideration holds true for a glitch showing on the (reconstructed) vertical component. In Text S2, we have detailed our theoretical considerations of apparent glitch polarizations especially with respect to the nonorthogonal sensor orientations of both the VBB and SP seismometers. There we demonstrate that our polarization analysis is correct and that some resulting angles may not be intuitive for special cases for VBB and SP.
Figures 3-5 demonstrate the polarization analysis of the VBB and SP glitches for 2019. The plots incorporate the VBB data at 20 and 10 sps (channels 02.BH? and 03.BH?, respectively, "?" stands for the components U, V, and W), and the SP data at 20 and 10 sps (channels 67.SH? and 68.SH?, respectively). These are the data rates that, depending on the actual satellite down-link capacities, are continuously returned to Earth. Besides some minor data gaps in the continuous operation, there is a large period with no data return between Sols 267-288 (28 August 2019 to 19 September 2019). This is due to the solar conjunction period where Earth-Mars communications were obscured by the Sun as consequence of their relative orbital positions. With respect to the Local Mean Solar Time (LMST, local InSight time, e.g., Allison & McEwen, 2000), the polarization patterns prevail over many sols and we discuss some of them in the following to understand the glitch behavior in more detail. First, we consider glitches occurring on only one VBB or SP component before building our arguments for multicomponent glitches. We conclude the section by looking at glitches that occurred simultaneously on VBB and SP. Note that all details concerning the SEIS sensor assembly and available SEIS channels can be found in Lognonné et al. (2019).

Glitches on Only One Seismometer Component
For  Note the rate change of nighttime glitches (red colors) after heater activation (Sol 168), (d) SP glitch azimuths. Pattern 5, that also occurs on VBB, is marked. The blue dots mostly refer to false glitch detections caused by HP 3 -hammering sessions and InSight's robotic arm movements, (e) SP glitch incidence angles, almost exclusively pointing in the horizontal plane. This demonstrates that multicomponent SP glitches occur only among the horizontal SP V and SP W components. Color code is same as in Figure 3.

10.1029/2020EA001317
Earth and Space Science SCHOLZ ET AL.
can only be related to instrumental artifacts such as (but not limited to) thermally driven stress relaxations in the suspension spring or pivot, displacement of one of the fixed plates of the displacement transducer, voltage offsets in the individual feedback electronics, or tilting of the individual sensor within the SEIS frame. Figures 3a and 3b show the VBB one-component glitches. For most identifiable patterns we find their behavior clearly changed either when the SEIS heaters were turned on (these are mounted on the leveling ring, see Lognonné et al., 2019) on Sol 168 (19 May 2019), or after the solar conjunction period in which the heaters were off and the SEIS instrument cooled down. This plus the fact these glitch patterns emerge due to their recurrence with respect to the local time, that is, repetitively at the same time of the sol, leads us to conclude that they are indeed thermally driven. What we suspect is that the enormous 10.1029/2020EA001317

Earth and Space Science
Martian surface temperature changes, that can reach up 100 K each sol, introduce stresses into the material-possibly within the Evacuated Container. Even though the temperatures inside SEIS do not vary as much as outside, the stresses grow and are released once a critical temperature is reached, thereby producing a glitch. When the heaters are on, the SEIS's thermal regime exhibits essentially higher temperatures and, in second-order, lower diurnal amplitudes and thermal spatial gradients. This contributes to minimize thermal stresses in this complex assembly, thus diminishing or at least altering glitch production. We demonstrate heater-related glitch behavior in more detail in the next section 3.2 for multicomponent glitches. We have no good explanation why we observe so many more glitches on VBB W compared to the other two VBB components, especially after the conjunction period during which the SEIS heaters were off. Only after ≈100 sols after the conjunction the number of one-component glitches (mostly constituted by glitches on VBB W) returns to the preconjunction level ( Figure 3b).
For SP, a glitch occurring on only one single component could potentially be interpreted as the SEIS instrument tilting if the glitch shows one of the two horizontal components, SP V (2) or SP W (3). The tilt direction must furthermore be orthogonal to the other horizontal component so the glitch could only be seen on one component. More plausible than being caused by SEIS tilt is that these glitches are also thermally driven. Figure 3c demonstrates that the horizontal one-component SP glitches change their behavior/occurrence with heater activation. For SP U, oriented almost vertically, a one-component glitch cannot be explained by instrument tilt because it does not point in the horizontal plane. These glitches therefore must relate to effects on the sensor level. Interestingly, Figure 3d demonstrates that SP U glitches that occur during the morning hours, that is, when the environment becomes warmer, point upward while during the evening/night hours, that is, during the cooling cycle, the glitches point downward. We interpret this behavior as further evidence for the thermally driven nature of one-component glitches.

Glitches on Multiple Seismometer Components
The multicomponent glitches for VBB and SP are illustrated in Figure 4. Especially for VBB, for which we generally detect more glitches, clear patterns emerge over the period of 2019. We discuss five of these patterns in the following.
We observe a glitch pattern with associated acceleration change pointing toward North (blue dots, Pattern 1). These three-component glitches are often accompanied by glitch spikes and occur around 1800 LMST and thus when the local temperatures start dropping. The incidence angles are ≈90°(in the horizontal plane) and hence may represent the SEIS instrument tilting. For this glitch pattern, however, we observe an additional 4.2 Hz ringing in some cases for the duration of the glitch, something not expected for an unhindered SEIS tilt. This occasional ringing could be related to other short duration data artifacts ("donks", Ceylan et al., 2020) we observe mostly in data with higher sampling frequencies (>20 sps). Due to the apparent temperature dependence of this glitch pattern we currently favor the possibility that they are produced by the temperature decrease resulting in slight contractions of the tether and/or Load Shunt Assembly (LSA)located both at azimuths ≈15°and connecting SEIS with the InSight lander. This argument is supported by the fact that the heater activation on 19 May 2019 (Sol 168) seemed to have no significant effect on these glitches (Figure 4c), bearing in mind that the heaters are located within SEIS and the LSA/tether is not. Furthermore, the largest of these VBB glitches (amplitudes larger than 10 −7 m/s) are also observed on SP with agreeing glitch azimuths and incidence angles ( Figure 5) and the same 4.2 Hz ringing. It therefore could be concluded that this glitch pattern is indeed due the SEIS instrument tilting, caused by cooling effects of the tether and/or LSA that also cause the 4.2 Hz ringing. On the other hand, the glitch azimuths of Pattern 1 average to ≈0°and not ≈15°where the LSA/tether are located. Also, the acceleration changes associated with these glitches point northward and hence suggest SEIS tilting southward, something difficult to reconcile with, for example, the contracting tether "pulling" SEIS toward the northern direction. One may therefore suspect not the tether itself as possible glitch cause but instead its connection with SEIS. Interestingly, there is another glitch pattern (green dots, Pattern 2) with similar features: azimuths pointing consistently south (instead of north), incidence angles of ≈90°, often preceding glitch spikes, occurrence at about 0800 LMST (instead of 1800), occasional 4.2 Hz ringing during the glitch, no significant effect of heater activation on glitch amount, and the largest among them also visible on SP with coinciding azimuths and incidence angles ( Figure 5). This pattern could represent the counterpart to Pattern 1; in the warming cycle of the sol the glitch cause reverses.

Earth and Space Science
The glitches with azimuths ≈240°occurring around 2100 LMST (pink dots, Pattern 3) show clear indications of being thermally driven. These three-component glitches with accompanying glitch spikes, that are not seen on SP, appear just after SEIS heater activation while before they were absent. Their consistent incidence angles of ≈100°prohibit their interpretation of SEIS tilting but instead point toward a thermal effect acting on all VBB sensors. After the conjunction period, during which the heaters were off, they do not immediately reappear with the heater reactivation but only ≈30 sols later together with azimuths being more variable. Such conjunction-delayed behavior (before the preconjunction state is reached again) is also readily visible for other multicomponent patterns during the night time (red and pinks dots at azimuths of ≈40°). For these reasons, such glitch patterns are likely to represent SEIS-internal, thermal effects. This is further supported by the glitch histogram in Figure 4e that clearly shows reduced glitches for the nighttime just after heater activation (fewer red dots). We note that there is a similar pattern on SP at azimuths of ≈350°(red dots) that occurs at the same times as the corresponding VBB one.
Another prominent VBB multicomponent glitch pattern occurs in the early sol hours with azimuths mostly due East (yellow-orange dots, Pattern 4). These three-component glitches with accompanying glitch spikes, that are not seen on SP, happen during the diurnal cooling cycle. Although there seems to be no obvious influence by the heater activation (or reactivation after conjunction), with increasing sols they occur at earlier hours. This plus the fact that their incidence angles INC ≠ 90°exclude a rigid tilt of the SEIS instrument lets us conclude that for this pattern, too, thermal effects are the primary glitch cause.
There is another thermally-driven glitch pattern that appears on both VBB and SP in the early morning (yellow-orange-red dots, Pattern 5), which again leads to glitches on the vertical VBB component (INC ≠ 90°). It is discussed in detail in the next section 3.3.
Patterns 3-5 are therefore all associated with nonhorizontal incidence angles suggesting that the three VBB sensors are not detecting an overall instrument tilt. Instead, each of the three VBBs detects a different tilt that consequently leads to the nonzero glitch on the vertical axis. The VBB sensors are mounted on a titanium plate inside the Evacuated Container through three mounting bolts oriented at azimuths of 105°( IF1), 225°(IF2), and 345°(IF3). So the first one is pointing roughly due east, while the two other ones point due west and are symmetrically to one another with respect to the west. This configuration produces colder temperatures on the east side during the night than on the west side (and the opposite during the day), with larger gradients between IF1-IF2 or IF1-IF3 than between IF2-IF3. This is likely the primary source of these thermal glitch patterns. We note that the temperatures between the inside and outside of the Evacuated Container are out of phase with the outside being ahead by about 7-9 hr (Pou et al., 2019). Figure 5 shows all glitches that occurred within ±2 s on both VBB and SP. From these 638 glitches, 118 glitches reveal the same azimuths to within ±10°. Most of the glitches on VBB and SP that match in azimuth were discussed already in the previous section 3.2 (green and blue dots, parts of Patterns 1 and 2). As we pointed out, these glitches show incidence angles of ≈90°for both VBB and SP and therefore could signify the whole SEIS instrument tilting.

Glitches on Both VBB and SP
The most prominent glitch pattern in Figure 5 is the one at azimuths of ≈145°for VBB and ≈110°for SP (yellow-orange-red dots, Pattern 5). From the beginning of SEIS's operational mode, these relatively strong glitches occurred once every morning with persistent glitch azimuths throughout 2019. Between sols 80 and 167, so before SEIS's heater activation, their onset times shift each sol by on average 4 Martian minutes (around 2% longer than SI minutes). This can be interpreted as the glitches occurring at a critical temperature during the diurnal cooling cycle that is reached earlier every sol as the northern hemisphere (where InSight is) is entering the colder season. When the heaters were turned on, leading to SEIS being in a thermally mitigated state, the glitches continued drifting toward earlier times but now with an average rate of less than 2 min per sol. After the conjunction period, during which the heaters were turned off, we observe the same as for many other glitch patterns; a more diffuse signature of the glitch azimuths and incidence angles that seem to return to preconjunction states only around 100 sols later. Also, the onsets time now drift toward later times (red to yellow) each sol which interestingly coincides with the fact that the Martian solstice occurred just after the conjunction on Sol 308. For this pattern as a whole, we were able to clearly identify −54°C as the critical temperature around which the glitches occur (Figures 5d and 5f). In addition for 10.1029/2020EA001317

Earth and Space Science
VBB, there are more patterns with similar behavior for which we could find the critical temperatures; these correspond to Pattern 3 (red and pink dots, section 3.2). All this evidence once more supports the fact that most glitches are thermally caused. Note that the temperature sensor we used here is scientific temperature sensor A (SCIT A, channel 03.VKI), located at the northern, inner side of leveling support structure. The temperatures measured at this sensor can also occur elsewhere in the SEIS assembly at the same time.

Glitch Removal
Once a glitch has been detected (section 2), the raw waveforms are modeled as a linear combination of the glitch-the response of the seismometer to a step in acceleration-and the glitch spikethe response of the seismometer to a step in displacement. The two responses can be modeled from the poles and zeros of the transfer function of either the VBB or SP seismometer. Only the amplitudes and the precise timing of the source (which might be between two recorded samples) are to be inverted with such model. Due to the time-limited extent of glitches and spikes as opposed to permanent (everlasting) steps in acceleration and displacement, respectively, all methods prefer to correct the raw data rather than the data after conversion to physical units.
The MPS group models a glitch waveform for each detected glitch using three parameters: an amplitude scaling factor, an offset, and a linear trend parameter. To find the best fit, the model is iterated over the data for each (sub)sample and the best fit for the three parameters is determined using nonlinear least squares Figure 6. Automated glitch removal for VBB at 20 sps: (a) We fitted the glitches in the data (blue lines) with the nominal VBB responses to a step in acceleration (red lines). The deglitched data (black lines) were obtained by subtracting only the scaled version of the synthetic glitches from the original data, that is, without offset and linear trend parameters. (b) High-frequency spikes (red lines) were modeled with the nominal VBB responses to a step in displacement and fitted to the deglitched data of (a) (blue lines). Our glitch model allows to fit both the glitch and the glitch spikes very well, even if small mismatches remain. Dashed vertical lines: theoretical onsets identical for glitch and spike; a: calculated amplitudes of acceleration and displacement steps corresponding to the glitches and spikes, respectively; t: time difference between calculated glitch and spike onsets smaller than sampling period (subsample fitting); VR: achieved total variance reduction.
(NLSQ, via the Trust Region Reflective algorithm). The deglitched data then is obtained by subtracting the fitted glitch without the offset and linear trend from the original data. To avoid introducing tiny DC offsets in the data caused by the fact that glitches are not yet fully returned to their baseline after, for example, 30 s, the fitted glitch is not only removed for the fit windows but for time windows corresponding to 10,000 subsequent samples (independent of the data sampling period). The same procedure is done for glitch spikes once a glitch has been removed; however, the subsample search grid is finer than for glitches because it has greater impact on the goodness of fits. To prevent our method from removing data where the glitch fit is not good enough, that is, the model is fitted to data that are in fact no glitches or fitted to glitches that cannot fully be represented by our model of a step in acceleration, we correct glitches only for which we can achieve a variance reduction of, for example, >80% with respect to the glitch fit window. We find this threshold to generally permit the removal of all large glitches, leading to a major reduction of all acceleration steps associated to glitches ( Figures SI2-4a and SI2-4b). Small glitches are also removed if their waveforms represent that of the underlying model well. For cases where such glitch fits do not work well, we repeat the approach but allow for a finite rise time of the acceleration change (as opposed to a zero rise time acceleration step). This does not change the resulting waveform of the glitch model too much while improving the data fits in some cases. We note that this limited ramp, that is, usually less than 5 s in length, is linear, a Gaussian-like ramp does not improve the fits. For spikes and their corresponding steps in displacement such finite rise modeling should not be done as it changes the resulting spike waveforms drastically. The MPS method is implemented for all sampling frequencies. An example of its glitch removal is shown in Figure 6.
The UCLA group carries out glitch and spike removal on 10/20 sps data. Some glitches show symmetric or asymmetric broadening relative to the glitch template, suggesting the source function is more complicated than a Heaviside step in acceleration. As a first approximation, convolution with a unit Gaussian or All four methods remove the glitch sufficiently however not fully equally. We note that this glitch is a prime example of glitches that do not perfectly fit our step model of acceleration but show a slightly broader response that calls for adaptions in our removal algorithms (section 4). exponential decay, which adds an extra parameter, significantly improves the fit, but runs the risk of overfitting data. To minimize this effect, the approach is only applied to data that show >0.9 correlation coefficient with the glitch corresponding to our acceleration step model. Glitch (sometimes broadened) and spike templates were fit to the glitches and spikes, respectively, using NLSQ. Because of the delta-like The data are shown before and after temperature decorrelation (TD), the latter of which is needed to hunt for Phobos' tidal signal in the SEIS data (Pou et al., 2019;Van Hoolst et al., 2003). The deglitched data (DG, ISAE method) after TD show reduced spectral peaks caused by glitches. This is true for both time spans shown, indicating our deglitching is stable over different periods and improves the data quality. (b) Comparison of raw data and deglitched data (UCLA method, see also Figure 7) and their Ps receiver functions for marsquake S0173a. (top panels) Waveform data (vertical, radial, transverse components) around P wave onset of S0173a, band-pass filtered between 0.1 and 0.8 Hz. The prominent glitch around 20 s is still dominating the horizontal components. (gray boxes) Time window used for the deconvolution in Ps receiver functions in lower panels: the long-period contamination by the glitch becomes apparent after 8 s on the horizontal components, masks any later arrivals, and casts doubts on the reliability of earlier phases. For example, an additional arrival near 7 s is clearly visible on the radial component, a phase that is also observed in receiver functions for other marsquakes not contaminated by glitches (Lognonné et al., 2020, Supplement IV).
shape of the spike over one or two sample intervals, the starting model must find the location to within a fraction of a sample interval (e.g., 0.05 s). Glitches are easier to fit than spikes, being low frequency, and requirements on the starting model are less stringent. Spikes are much smaller in 2 sps data relative to glitch sizes. Thus 2 sps data were used to generate a glitch catalog (section 2). The starting parameters from the 2 sps fits were then used to fit glitches in the 20 sps data and residuals were calculated. The residuals were examined for the presence of a spike in the data before the glitch peak, by requiring its amplitude to be greater than 5 standard deviations of the residuals after the peak. If true, an iterative forward model was run by shifting the phase of the spike template about the corresponding peak in the residuals (in steps of sample interval/10), and finding the amplitude and phase of maximum cross correlation. The NLSQ was run again with both spike and glitch templates, and the result checked whether cross-correlation of data and model are above a threshold, and if so, the results are stored. At this stage, for polyglitches (one on top of another) we search for spikes throughout the sequence. Even though a number of spikes have been removed, there are residuals and transients that remain. Polyglitches can have several internal spikes, and extreme glitch overlap, making automatic procedures difficult, requiring manual fitting.
The removal algorithm of the ISAE group is basically described in section 2.2 (glitch detection). Once a glitch has been detected using cross correlations between the model and data, the model without linear trend and offset is subtracted from the data. This method is implemented for all sampling frequencies available. Spike removal and deviations from the simplified acceleration step model are not implemented.
The IPGP group inverts three consecutive acceleration step sources for the glitch which allows not only to invert for multicomponent glitches occurring within these three samples but also to invert for the phase delay through finite-difference approximation of the first and second time derivative. This linear approach allows the inversion to provide identical results in the U, V, W coordinates or in the Z, N, E coordinates, as the rotation between the two coordinates systems is a linear relation. Conversely, the three other methods, through their nonlinear part of the inversion or through the cross-correlation phase fitting, have built-in small reasons to provide different solutions depending on the coordinate systems.
In the end, all the proposed deglitching methods are nevertheless based on the same idea of assuming a step in acceleration and displacement to model a glitch and spike, respectively, by using the instrument impulse response of either the VBB or SP seismometer. Removal differences across the methods are mostly due to thresholds below which a glitch is removed or not, and by how these methods attempt to fit glitches that do not fully correspond to our acceleration step model. No general rule on the thresholds can be provided as they depend on the data processing target. As an example, all methods provide similar deglitching for the large glitches occurring during the cooling periods and during the night. More freedom is available for fitting longer source duration glitches during the day although some may represent the real response of SEIS to small pressure drops (section 6.3) which can generate nanotilts of the SEIS instrument. At the same time, while many glitch spikes are fitted by the templates, there are a significant number that have quite different morphology, longer ringing, or longer-period transient behavior. Caution must therefore be exerted when attempting to remove these as it may unintentionally lead to removal of small parts of higher-frequency content. To circumvent such effect, spike fitting should only be attempted within a few samples left and right near the theoretical glitch onset. Subsequent removal should only occur if the fits are good enough. Due to the spikes' delta-like overall shape, we argue that this procedure diminishes any unwanted removal of signal. It is also possible to filter the spikes out rather than remove them from the raw data; however, small artifacts depending on the exact case may remain. All these arguments combined is the reason we do not provide glitch and/or spike corrected data for all available periods but instead make our codes available, enabling own comparisons and removal choices to those interested. An example of glitch removal showing all four methods is demonstrated in Figure 7 for two glitches occurring during marsquake S0173a. In Figure 8 we demonstrate the positive effects deglitching can have for the case of VBB long-period spectra to hunt for Phobos' tidal signals and for the case of receiver functions of marsquake S0173a.
We lastly point out that we have discontinued our deglitching efforts using the stationary wavelet transform as described in the Supplement V of Lognonné et al. (2020). While this approach provided promising and correct results for a fair amount of cases (as far as one can tell), there is no underlying, physical model involved and the implicit data "correction" therefore seemed too arbitrary. For many cases this approach further introduced DC offsets in the deglitched data whose amplitudes and lengths depended on the length of data read (and therefore maximum decomposition level); an artifact that we could never manage to fully avoid.

Glitch Model
Throughout this paper we have assumed that glitches can be understood as steps in acceleration and glitch spikes as steps in displacement. This model allowed us to successfully detect, analyze, and remove one-component and multicomponent glitches for both VBB and SP. In the following we detail the theoretical considerations behind this simple model.
Let us assume glitches are caused by a small instantaneous tilt. By instantaneous we mean that the time history of the tilting is so short that it cannot be resolved with any given sampling frequency available to us (maximum 100 sps). We are thus allowed to idealize any step in time by a Heaviside function. Physically, such short instantaneous events can, for example, be the result of stick-slip events.
The small tilt is assumed to be the result of a rotation around a horizontal axis, a ! . Recall that the VBB is a pendulum seismometer where the (inverted) pendulum is constrained to rotate around a horizontal axis, b ! .
The sensitivity direction, s ! , of the pendulum is perpendicular to the b ! axis and is inclined relative to the horizontal plane by a dip angle of δ ¼ −29:3°. Let us also assume for simplicity that all the mass of the pendulum is concentrated in its center of gravity (CoG)-which would be the case for a mathematical pendulum. Now we can distinguish five cases which differ by the location of the accelerometer relative to the tilt axis, a ! : (1) the two axes a ! and b ! are parallel and a ! passes through CoG: In this case the accelerometer gets only reoriented relative to the gravity vector but the CoG stays in place.
(2) the two axes are parallel and a ! does not pass through CoG but is at the same height as the CoG: In this case the accelerometer gets displaced vertically and reoriented relative to the gravity vector. However this reorientation is negligible because it is only a second-order effect.
(3) the two axes are parallel and a ! does not pass through CoG. Furthermore, a line parallel to s ! passing through CoG intersects with a ! . In this case the accelerometer gets displaced vertically and reoriented.
However, the displacement is in the direction perpendicular to the axis and hence is not seen by the accelerometer. Only the reorientation is sensed. (4) For all other locations of the rotation axis a ! for which a ! and b ! are parallel, the accelerometer will see both a displacement and a reorientation relative to the gravity vector. (5) For the general case where a ! and b ! are not parallel, the same arguments can be made but the effect sensed for a given tilt angle will always be reduced relative to the case with parallel axes a ! and b ! since the tilting is reduced.
As soon as the accelerometer gets reoriented relative to the gravity vector we expect to see the response due to a step in acceleration, because the projection of the gravity vector into the sensitive direction is changed. In those cases where the accelerometer gets displaced we expect to see the response due to a step in displacement. The five cases then only differ in the relative size of the displacement and tilting.
What do these signals look like? In Figure 6 we have plotted the response of the VBB sensors to a step in acceleration and the response to a step in displacement, both including the effects of the limited pass-band and downsampling. To model the instrument responses to these steps, we take the full seed response and evaluate it at the frequencies corresponding to those of the Fourier transform of the input steps using evalresp-a piece of software provided by the Data Management Center of the Incorporated Research Institutions for Seismology (DMC/IRIS). Figure 6 also demonstrates how we can use the modeled glitch and spike to remove them from the data.
Can these signals explain the data? As Figure 6 demonstrates, the modeled responses have been shifted in time and scaled to match the data. The fit is excellent both for the low-frequency glitch and the high-frequency spike. We take this as confirmation that our simple model is capable of explaining the glitch waveform with four parameters: start time and amplitude of the step in acceleration plus the start time and amplitude of the step in displacement. In fact we could show that the start times of the acceleration and displacement steps coincide to the millisecond-which is what our model predicts. Thus, we only need three parameters: the start time and the amplitudes in displacement and acceleration. Determining the start time requires an excellent calibration of the high-frequency part of the sensors transfer functions, as well as high sampling rate. While deglitching on the 20 sps data is therefore much more precise and has been done for two of the described methods (MPS, ISAE), the deglitching on lower rate data, for example, 10 SPS (UCLA) or even 2 sps (IPGP) can be achieved, including for the spike amplitude, however, with the signal-to-noise ratio reduced by the frequency ratio of the bandwidth. Fitting the spike plus glitch with these three parameters implies determining the start time to subsample resolution. We provide a more mathematical description of our model for the glitch plus spikes phenomenon in Text S2.

Other Observations
In the following we briefly discuss other aspects of glitches and spikes that we encountered during our investigations. This section shall therefore complement our understanding of glitches and detail some more implications.

Possibly Locating SEIS-Internal Tilts
Our glitch model presented in section 5 is valid for rotations of the sensor assembly as a whole (e.g., caused by a change at one foot of the sensor assembly), for just the VBB sensors (e.g., caused by stick-slip events originating at the interface between the Evacuated Container and the leveling support structure), but also for an individual sensor (e.g., caused by stick-slip events originating at the sensor-support interface or at the fixed side of the pivot or spring). Each of these cases implies a different value of r: the distance between VBB U to the sensor assembly feet at 16 or 21 cm (Fayon et al., 2018), or the distance from the sensor's center of gravity to its pivot with 2.6 cm (Lognonné et al., 2019).
We illustrate this geometry with the glitch example of Figure 6 recall the glitch and spike characteristics in Table 2. This glitch has a vertical component and can therefore not represent the SEIS instrument tilting as a whole. The azimuth of the glitch opposite (opposite of acceleration) and of the spike (displacement) are 219°a nd 228°, respectively. These values average 223.5°, which is quite close to one of the plate's mounting bolts IF2, located at 225°. The opposite signs of the glitch amplitudes of VBB V and VBB W suggest a deformation relatively symmetrical with respect to the IF2 azimuth, while the low-amplitude glitch on VBB U suggests the latter to be much reduced between the two other IFs. This glitch is therefore compatible with a radial deformation of the mounting bolts IF2. Further analysis on the impact of the thermo-elastic stresses in the VBB sphere and the resultant glitch generation will however be demonstrated in a future publication.
During the night, a wide range rotation radii are found, likely resulting from internal deformation of the Evacuated Container triggered by thermal effects, as discussed previously. During the day, however, the rotation radii of the glitches are more stable and in the range 10-30 cm, suggesting an external source and therefore rigid tilt of SEIS, likely generated by the atmospheric activity.

Loading With the Arm
The InSight mission includes the Heatflow and Physical Properties Probe (HP 3 , Spohn et al., 2018) that carries a probe ("mole") intended to hammer itself 3-5 m into the Martian regolith. The mole has had difficulty getting started, and so the lander's Instrument Deployment Arm (IDA) has been pressed into service to help. On several occasions, the IDA has pushed down on either the regolith or the mole itself. When the IDA pushes down, it induces an elastic response in the regolith, deforming the surface into a funnel shape, inducing a tilt at the seismometer about 1.2 m away. This tilt of about 70 nrad is clearly observable on both the SP and VBB sensors in Figure 9 as steps in the horizontal accelerations.
In this example, at the start of the command sequence the IDA was pushing down lightly on the mole and was given four commands: (1) move up to get off the mole, (2) move radially outward slightly, (3) move down to just above the mole, and (4) move down again to reload the mole with a downward force. We see in the Earth and Space Science seismometer data the first move up and the resulting tilt up to the NE. The arm resonates after it loses contact with the mole, and we see that as the 4.2 Hz ringing in the seismometer data. The seismometer does not have a significant response to the radial outward move. Then on the third move, it appears that the IDA actually touched the mole while stopping and then rebounded and resonated while hovering in midair just above the mole. Finally, the IDA moves down to load the mole and we see a tilt down to the NE at the seismometer.
We also observe several glitches (accleration steps) that happen at the same time as the IDA motions (red circles). We interestingly observe that the BHE component shows steps of the same sign for both the arm loading and unloading. Two of the glitches further appear to involve the whole sensor assembly as they are seen on both the VBB and SP. Other glitches seem to be limited to one or more components of the VBB. This all points toward that these glitches are internally caused and only triggered by the IDA movement. Attempting to remove these IDA-induced glitches show convincing fits with our acceleration step model for the BHV and BHW components; however, for the BHU component the removal is more difficult also because of the additional 4.2 Hz ringing (Figure 9b, top panels). Nevertheless, IDA movements are limited and therefore this type of glitch does not represent a major contamination of InSight's seismic data.

10.1029/2020EA001317
Earth and Space Science pressure effects are generating signals on the SEIS components mostly from 0.5 mHz up to about 2 Hz, among which convective vortices are generating the largest physical signals observed by SEIS. Their dominant period, as seen both by atmospheric pressure sensor and SEIS, can be close to the one of the glitches depending on their size, distance to SEIS, and wind speed.
At frequencies lower than 0.1 Hz, the compliance response of the ground is dominated by tilt effects which are strongly impacting SEIS's horizontal components . These ground responses are usually more complicated than our simple acceleration step model (Murdoch et al., 2017). We instead often observe that the dust devils' pressure signal convolved with the instrument response of SEIS can match well in shape with the integrated raw waveforms of the observed SEIS glitches. Such ground responses are the reason SEIS signals induced by convective vortices may, wrongly, be detected as glitches. On top of these complexities, the ground deformations induced by convective vortices are sometimes generating real glitches (SEIS's raw data matching perfectly with our acceleration step model) that can even show on the vertical components. Discriminating between these various SEIS signals is therefore a challenge for all glitch detection methods.

Glitch Causes
As we established in section 3, there are likely two main categories of glitches: (1) instrument-related glitches caused by thermal variability across the different parts of the SEIS sensor assembly, and (2) glitches caused by true, rigid tilt of the SEIS instrument as a whole. In the following we discuss both types in more detail.
The majority of glitches is related to internal instrument effects. We could conclude so as these glitches have incidences angles pointing typically away from the horizontal plane which excludes their interpretation as true instrument tilt. Our analysis demonstrates that these glitches are closely linked to temperature or even specific temperatures. This comes as no surprise, on Mars the SEIS sensor assembly is installed in a harsh environment with the daily atmospheric temperature varying by 80-100 K. Although these temperature variations are attenuated by the Wind-Thermal Shield, Remote Warm Enclosure Box and Evacuated Container, and further mitigated by heaters and the Thermal-Compensation Device (the latter as part of the VBB sensors), daily temperature variations of 15 K are still reached inside the Evacuated Container hosting the three VBB sensors. These temperature fluctuations inevitably lead to thermal strains and thermally induced stresses at the contacts between materials with different thermal expansion coefficients. These stresses can be relaxed by a variety of irreversible mechanisms such as creep, diffusion of lattice dislocations, or stick-slip along mechanical contacts. While we do not know which actual mechanism or which combination of mechanisms is at play, we attribute thermally related glitches to intermittent stress relaxation events, for example, stick-slip events. Furthermore, as shown in Figure  10a, we observe nonhorizontally polarized glitches to occur more during the diurnal cooling cycle while in the warming cycle the amount of glitches reduces more and more until the temperatures start rising again. This correlation is closest with the temperatures of the VBB sensors (red curve).  Figure 11c for sketch), and the three VBB temperature sensors inside the Evacuated Container (purple, brown, pink). The background colors indicate the diurnal cooling cycle (blue) and warming cycle (orange) of the VBB sensors. (a) Only glitches with an INC ≠ 90°are shown, that is, not polarized in the horizontal plane and therefore only explainable by instrumental causes, most likely by thermal effects. During the VBB warming cycle the amount of glitches reduces while in the cooling cycle an increase is observed. (b) Only glitches with an INC ≈ 90°are shown, that is, polarized in the horizontal plane and therefore the only glitches that could by explained by rigid motion of the whole SEIS instrument. Closely linked to the atmospheric temperature, two peaks are visible at 0800 and 1800 local InSight hour, corresponding to Patterns 2 and 1 in Figure 4, respectively. We suspect these glitches to be related to the LSA/tether system pushing and pulling on SEIS, driven by diurnal changes of atmospheric wind, pressure, and/or temperature.
As explanation we suggest that in many places within SEIS layered materials are in frictional contact. When cooling, the outer layers cool first and hence contract against the inner layers. This enhances normal forces and makes for both larger and more frequent stick-slip events. Conversely, in the warming cycle, the outer layers heat up first and therefore expand away from the inner layers, reducing or possibly eliminating the normal forces implying smaller and possibly less frequent stick-slip events.
The question of whether external events can trigger glitches arises when we inspect marsquake S0173a (Figure 7), the VBB response to certain pressure drops, or the VBB response to ground loading experiments with the scoop of the instrument deployment arm (IDA, see Figure 9). In all these cases the seismic waveforms are contaminated by a glitch. We argue that external events alone do not cause thermally-related glitches. Instead, as the SEIS sensor assembly goes through the daily temperature cycle, internal stresses build up until a threshold is reached, and a stick-slip or another stress relaxation event occurs. In other words, an infinitesimally small additional stress may suffice to trigger a glitch if it occurs at the right time, that is, a time when thermal stresses have almost reached the critical threshold and a relaxation event is about to happen. Any additional external acceleration, be it a marsquake, the passage of a pressure drop, an IDA arm movement, or a soil loading experiment with the IDA scoop will make the glitch occur earlier than it would have without the external event. So in this view external events alone do not cause glitches, they merely advance their time of occurrence. To look at this closer, we analyzed the delays between arrivals of seismic events and glitches detected shortly after them. Since the broadband and low-frequency marsquakes are suspected to be due to a stationary Poisson process (personal communication Martin Knapmeyer) while glitches are distributed unevenly over the sol, one may suspect that triggered glitches occur already within a few seconds after an arrival if following our model. We found no obvious relation ( Figure SI2-2). While the number of events with clear P and S arrivals is small, and a more thorough reanalysis with a larger data set may be worthwhile, all our analyses combined still suggest that the timing of glitches generally has a strong stochastic component next to a deterministic component. This is further supported by the frequency-amplitude distributions of glitches that seemingly follow a Gutenberg-Richter relation ( Figure SI2-3), and the presence of the diurnal harmonic and all its integer multiples in a time series composed of modeled glitches ( Figure SI2-4c).
On the other hand, around 30% of all glitches exhibit quasi-horizontal polarization and thus could represent the whole SEIS instrument tilting. Some of these cases may indeed be rooted in the ground tilting and thus be real seismic signal, a scenario demonstrated by Zahradnik and Plesinger (2005). They found glitches (they use the term "long-period disturbances") during earthquakes phase arrivals solely to occur on the horizontal components, something we also observe for marsquakes but only for a few cases. They preferably interpreted such glitches as ground tilt, possibly caused by small-scale material instabilities beneath the station triggered by the incoming waves or thermally or chemically induced microcracks that would not require any incoming wave energy. These interpretations of tilt causes, however, are not unique and our investigations did not allow us to narrow down their causes as the InSight setup puts too many variables in question. For example, next to true ground tilt it is further conceivable that horizontally polarized glitches are caused by the SEIS instrument tilting either due to imperfect anchoring of its feet to the ground or by the Load Shunt Assembly (LSA)/tether pushing and pulling on SEIS as reaction to atmospheric changes in temperature, pressure, and wind. We have no clear observation that azimuths of such glitches cluster toward the feet of SEIS's leveling system (Figures 3-5); however, we cannot finally conclude that the anchoring may not cause such glitches at all. Instead, we find most of the instrument tilt indicating glitches to point either north or south, that is, either close to the surface-exposed LSA-tether system or diametrically opposed (Figure 4, Patterns 1 and 2). Indeed, as shown in Figure 10b, these glitches correlate more closely with the atmospheric temperature, peaking at 0800 (Pattern 2) and 1800 (Pattern 1) local InSight hour, that is, right when the atmospheric temperatures rise and fall, respectively. While the picture for these two glitch patterns is not fully conclusive with respect to their polarization signs (section 3.2), there remains the suspicion that the LSA-tether system or even the lander exert influence on SEIS and promote glitch production via mechanisms for which we have no unique interpretations.
The two common mechanisms suspected responsible for the generation of glitches are depicted in Figures 11a and 11b. In Figure 11c we schematically illustrate the SEIS instrument and its major components.

Glitch Mitigation
Given the abundance of glitches and their influence on the data analysis, the question arises how glitches could be mitigated for future installations. For thermally related glitches, the most obvious action is to decrease the thermal amplitudes the seismic sensors are exposed to. While for the SEIS instrument great care was taken to achieve just that (Wind-Thermal Shield, Remote Warm Enclosure Box, Evacuated Container, heaters, Thermal-Compensation Device; see Lognonné et al., 2019), the daily temperature cycle still exceeds those of fine terrestrial stations by 4 orders of magnitude, that is, around 15 K compared to a few mK. Given the harsh environments typically found on extraterrestrial, planetary bodies, it may not be easy to achieve . We suspect such thermally-driven effects to be the primary reason for multicomponent glitches such as shown in Patterns 3-5 ( Figure 4). (b) SEIS tilt α, corresponding to a true, rigid motion of the whole instrument thus producing glitches with INC ≈ 90°. Patterns 1 and 2 ( Figure 4) are likely caused by this scenario. Note that in both cases the VBB sensors may experience a tilt and a displacement (sections 5 and 6.1), and that the tilt α is depicted as 10°for both cases but is in reality in the order of nanoradians. Similar considerations apply for the SP sensors. Green lines: moving pendulum parts; P: proof mass; δ: VBB sensor dip ≈ −30°. (c) Simplified top-down view of the SEIS instrument; LSA: Load Shunt Assembly decoupling the tether from SEIS; gray ring: support structure carrying the the SP sensors, heaters, and the SCIT-A (scientific temperature A) sensor; IF1-IF3: mounting bolts connecting the Evacuated Container of the VBB sensors with the ring structure; numbers 1-3: SEIS's three protrudable feet integrated in the support structure (= leveling system, Fayon et al., 2018); TCDs: Thermal-Compensation Devices (one for each VBB seismic sensor). Not depicted is the Wind-Thermal Shield put over the Remote Warm Enclosure Box. For details on the SEIS instrument, see Lognonné et al. (2019). SI1: Lists of glitches detected by our different methods. SI2: Theoretical considerations for apparent glitch polarizations, mathematical description of glitch plus spike origins, and additional figures.
higher thermal stability; however, it should be considered by engineers. We can only speculate as to the exact sources of glitch production within the instrument. While we have good candidates such as places close to the mounting bolts of the Evacuated Container, the fact remains that InSight's seismometers, especially the VBB, are complex devices consisting of many materials, joints, and connections. One way to approach thermal glitch reduction may therefore be to use fewer materials and thus minimize potential thermal conductivity gradients, stresses, and expansions. A last, ultimate step to achieve thermal stability would be to completely bury the instrument and possibly even the tether, but this may not be feasible for many reasons. Future engineering efforts may also benefit from testing similar instruments in cold environments with preferably low seismic noise conditions, perhaps utilizing infrared cameras to investigate causes of thermal glitches in greater detail. No such experiments were performed with the flight or spare model of SEIS. For glitches indicating instrument tilt, one way to mitigate could be to improve on the feet anchoring (Fayon et al., 2018;Lognonné et al., 2019) by usage of even more specialized feet shapes possibly digging into the ground, bury the tether, and/or by deploying the instrument on hard rock as opposed to InSight which rests on regolith (e.g., .

Summary
On the seismic records returned from Mars many data disturbances are observed. The most prominent kind, which we termed "glitches" (section 1, Figure 1), significantly complicate the analyses of the seismic data. We have developed a physical model for the generation of glitches and their associated high-frequency spikes that occur simultaneously with the glitch onsets. In this model, glitches represent steps in the acceleration sensed by the individual sensors convolved with the instrument responses while glitch spikes represent steps in the displacement sensed by the individual sensors convolved with the instrument responses. We used this model to develop different algorithms for the glitch detection that are all able to identify most of the high amplitude glitches for both the VBB and SP seismometers (section 2, Figure 2). Based on the model we further demonstrated that many glitches are thermally driven (section 3, Figures 3-5, and section 7, Figure 10), some of which we showed to occur at specific temperatures, for example, at −54°C inside the SEIS instrument. Thermal glitches can also be triggered by external events such as convective vortices or movements of Insight's robotic arm (section 6, Figure 9) and likely represent SEIS-internal tilts that differ among the individual sensors and hence produce glitches on the vertical components, an observation that cannot be reconciled with the SEIS instrument physically tilting. Only a portion of all observed glitches can be explained by a tilt of the whole SEIS package, either related to true ground tilt, imperfect feet anchoring, or the Load Shunt Assembly/tether pushing and pulling on the SEIS instrument as response to atmospheric changes in wind, temperature, and pressure during the Martian day. We depicted the two most common cases of glitch production and a sketch of SEIS and its major components in Figure 11.
While terrestrial data influenced by glitches and glitch spikes may simply be discarded due to their difficult handling, this represents no valid option for the seismic data returned from Mars. Especially glitches influence virtually all single station methods available to seismologists and therefore can bias many analyses. We therefore devoted much of our efforts to develop code for the glitch and spike removal (section 4, Figures 6-8). Our algorithms have proven successful in many cases for both seismometers VBB and SP. Of course, there remain glitches and spikes especially of smaller amplitudes that we cannot sufficiently well fit and therefore confidently remove. To correct for such glitches, nevertheless, we have slightly deviated from our step model in acceleration, that is, we introduced fits for nonzero rise times (MPS), for a combination of multiple source functions (UCLA), and for three consecutive acceleration steps of varying amplitudes (IPGP). The resulting glitch waveforms of these adaptations are still close to those corresponding to a zero-rise time acceleration step, allowing however to fit for glitches varying from our simple step model. It is these adaptations where the removal codes diverge from each other, but the basic concept and general success is equal (if the parameters are chosen accordingly). Analyzing differences in the removals can be worthwhile for glitches deviating from our step-like acceleration model in more important time periods such as during marsquakes. Overall, the removal of glitches following the approaches presented in this study allows to improve on the quality of seismic data and hence helps to accomplish InSight's high scientific goals.
As no glitch removal algorithm can warrant a perfect clean-up of all glitches and their spikes, we prefer to not provide a deglitched time series of all available data. Instead, we have assembled our algorithms for glitch detection and removal into one Python/ObsPy toolbox. Some convenient functions for data retrieval and handling are also implemented. The package further holds MATLAB scripts to perform glitch detection and removal tasks as presented. Its link is at this site (https://pss-gitlab.math.univ-paris-diderot.fr/data-processing-wg/seisglitch). Documentation is available. Together with this code we also provide deglitched data for a selection of seismic events.