Detection, analysis and removal of glitches from InSight’s 1 seismic data from Mars

The SEIS instrument package with the three very broad-band 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 for 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 precursors. These pulses, which we term “glitches”, can be modeled as the response of the instrument to a step in acceleration, while the precursors can be modeled as the response to a simultaneous step in displacement. We attribute the glitches primarily to SEIS-

on Earth data influenced by such disturbances are often discarded especially when coinciding with 86 earthquake phase arrivals (e.g. Zahradnik & Plesinger, 2005). This obviously represents no valid 87 option for the seismic data returned from Mars and hence the correct treatment of the glitches is of 88 high importance for the scientific analyses. The present study focuses on the detection, analysis and 89 removal of glitches and extends Supplement V of Lognonné et al. (2020). 90 Glitches 91 In the literature (e.g. Iwan et al., 1985;Zahradnik & Plesinger, 2005;Vacka et al., 2015) the 92 phenomenon we are investigating here is sometimes referred to as "long-period disturbances", "accel-93 eration offsets" or even "mice", all generally describing the same type of data disturbance. Throughout 94 the present publication, however, we choose to apply the term "glitch" to these disturbances as it has 95 been established as such since their first observations in InSight's seismic data and hence been com-96 municated so to a wider audience on various occasions. Whilst we are aware that the word glitch is 97 typically associated to more general data artefacts and alike, we indeed use it here to refer to specific, 98 clearly defined disturbances in the data. A glitch (Fig. 1b,d), thus, is a particular type of transient 99 instrumental self-noise that, in the raw time series data, appears as a high amplitude, one-sided pulse 100 with a duration controlled by the seismometer's transfer function. For the VBB sensors, which have 101 76% of critical damping, glitches have a fast rise time followed by an exponential decay with a small 102 (∼9%) overshoot before almost returning to the baseline after ∼25 s. For the SP sensors, that are 103 overdamped with 110% of critical damping, glitches have a similar rise time followed by a decay before 104 almost returning to the baseline after ∼50 s. Glitches may also occur before a previous glitch has suf-105 ficiently decayed. The highest order of such "poly-glitches" we observe to date is four. Glitches (and  Many glitches, furthermore, show a high-frequency signal at their very glitch beginning that lasts 119 around 40 samples regardless of the data sampling frequency. We refer to these initial oscillations as 120 "glitch spikes". These spikes occur simultaneously with the glitch onset for both VBB and SP (Fig. 121 1b,d). Glitch spikes do not represent artifacts caused by the on-board analog or digital electronics. To 122 facilitate the analysis of glitches and help deciphering their origins, we analyse these spikes as well. To automatically detect glitches on SEIS' VBB and SP raw data, several groups (MPS,ISAE,125 This detection algorithm, implemented in Python (Rossum, 1995) and ObsPy (Krischer et al.,137 2015; Beyreuther et al., 2010), performs the following processing steps on a given period of three-138 component seismic data (components U, V, W): (i) decimate the data to two samples per second 139 (SPS), allowing all data per seismometer to be run with the same parameters and enabling faster 140 computations, (ii) deconvolve the instrument response on each component and convert to acceleration, 141 (ii) band-pass filter the acceleration data (e.g. 10-1000 s), so the steps in acceleration emerge more 142 clearly, (iv) calculate the time derivative of the filtered acceleration data so the acceleration steps 143 become impulse-like signals, and (v) on this time-derivative, trigger glitches based on a constant 144 threshold. To avoid triggering on subsequent samples also exceeding the threshold but belonging to 145 the same glitch, we introduce a window length in which no further glitch can be triggered. This 146 parameter can be thought of as glitch minimum length. We note this parameter is smaller than the 147 typical glitch length for VBB and SP, allowing our detection algorithm to detect poly-glitches. 148 A glitch simultaneously occurring on multiple components is detected on each affected component 149 but the respective start times may slightly differ. However, after modeling of the full glitch waveform 150 (Section 4) we can retrospectively establish that such glitches occur at the same time to within 151 milliseconds. This holds true for all multi-component glitches observed to date on either VBB or SP, 152 also for data with the highest available sampling frequency of 100 Hz. Therefore, we declare as glitch 153 start time the earliest time detected across the UVW-components. The list of unified glitch starts 154 contains still many false-positive triggers caused by non-glitches with a steep enough acceleration 155 change to be triggered. This is because we choose to apply a constant threshold to the time derivative 156 of the filtered acceleration, rather than a threshold based on the current seismic noise level that 157 undergoes strong diurnal changes (amplitudes varying by a factor of 100 and more) dominated by 158 meteorological influences (e.g. Lognonné et al., 2020;Banfield et al., 2020). To circumvent, we rotate 159 the gain-corrected UVW raw data of the glitch windows into the geographical reference frame (ZNE-160 components) and perform a 3-D principle component analysis (e.g. Scholz et al., 2017). Theoretically 161 a glitch is linearly polarized as the associated vector of acceleration change is not varying, however 162 slightly altered only by seismic noise. Indeed, most glitches exhibit a high linear polarization >0.9 163 which we use to discriminate against other triggered signals. The polarization analysis further allows 164 to obtain the apparent glitch azimuth and incidence angles which we use to associate glitches with 165 particular glitch sources (Section 3). Visual inspection reveals the resulting glitch onsets are usually 166 accurate to within ±1 s (e.g. green lines in Fig. 1b,d). 167 2.2 Glitch Detection by Cross-Correlation with Impulse Response Function (ISAE) 168 The principle of this MATLAB-implemented detection algorithm is cross-correlation. It performs 169 the following processing steps on a given period of three-component raw seismic data (components 170 U, V, W): (i) a synthetic glitch is constructed by convolving the poles and zeros of the transfer 171 function of the VBB and SP sensors with a step in acceleration. To increase the temporal resolution 172 to sub-sample range, we synthesise several glitches each with a different sub-sample time shift; (ii) 173 while the frequencies above 2Hz are filtered, the long period variations of the data are extracted using 174 a low-pass filter with 10 −3 normalised cutoff frequency for VBB and 0.25 × 10 −4 normalised cutoff 175 frequency for SP. These are then subtracted from the signal (and added back at the end), before (iii) 176 the synthetic glitch is cross-correlated with the data. A glitch detection is triggered for the maxima 177 of the cross-correlation function that exceed a threshold a on a given component.

178
Another step is added to prevent non-detection of glitches or false-positives, depending on the 179 correlation threshold. For that, two thresholds are chosen: threshold a and threshold b, with a ≥ b.

180
The first step presented above is done for each component, with threshold a. Then, for each component, 181 a second cross-correlation with threshold b is implemented. For the times of every maximum of cross-182 correlation exceeding threshold b, we come back to the glitches detected on the other components 183 during the first step. If a glitch had indeed been detected at that specific time on another component, 184 a new glitch is declared on the component under study. We can therefore detect small glitches with 185 low signal-to-noise ratio when a strong glitch is detected at the same time on some other component.

186
In addition, in order to be able to detect poly-glitches, a second iteration of the detection algorithm 187 is performed after the glitches from the first iteration have been removed from the data. This MATLAB based method took into account that glitch amplitudes follow a power law dis-190 tribution with many more very small glitches than larger ones (see Fig 1 in electronic supplement).

191
Therefore the strategy was to remove the largest glitches first and repeat the process on the smaller 192 ones in an iterative procedure. In this method the raw UVW VEL channel data are inspected for 193 glitches and their spikes. The instrument response to a step in acceleration was termed "Green's func-194 tion." The 20 sps data were decimated to 2 sps and each channel was tested for correlation with the 195 response function as follows. An inverse filter was designed that turned glitches into narrow Gaussians 196 with rise times equal to the glitch so that each glitch represented one peak without the overshoot.

197
This enables detection of multiple close-spaced glitches. An STA/LTA (short time average / long 198 time average) ratio was found using convolution of the data with two box car functions separated by 199 more than a glitch window. The absolute value of band-passed data was tested for peaks above the 200 STA/LTA threshold. For the first iteration the STA/LTA was set large to remove the largest glitches.

201
The Green function was correlated with the data spanning a peak and if the correlation coefficient was 202 above 0.90 the detection was registered. If multiple peaks occurred close together, multiple Green's 203 functions were fit to the data using nonlinear least squares. The data was then cleaned by removing 204 the glitches. The orocess was then repeated lowering the STA/LTA threshold=7, and the new glitches 205 removed from the data. For the last iteration the STA/LTA threshold was set to 3, i.e. lowered again 206 and the correlation threshold was also lowered to 0.8. This removed many of the small glitches. Our 207 glitch detection is applicable to SEIS' VBB and SP sensors in both low and high gain modes. Implemented in MATLAB, this glitch detection method processes mostly 2 sps continuous data 210 and is therefore focused on long period continuous signals. It first removes the aseismic signals of each 211 raw axis by subtracting the trend and the first 12 sol-harmonics (i.e., up to 1/12 sol period, about 212 0.13 mHz in frequency). Then the three axes are equalized in digital units by convolving the V and W 213 channels by the convolution ratio of the U/V and U/W transfer functions, in order to correct for the 214 gain and transfer function differences between U, V and W. Note that this process also transforms an 215 impulse response in time on V and W into an impulse response with the U transfer function. As the 216 inversion (below) is a linear one, the glitch search and deglitching can be done either on the UVW or 217 on the ZNE rotated channels, with practically no differences for the inverted glitches.

218
The glitch detection is done first by identifying all extrema in the signal and then, for all found 219 extrema, least-square testing for the occurrence of a glitch using a modeled glitch. To model a glitch, 220 we convolve a step in acceleration not only for one sample (as all other methods) but for three 221 consecutive samples. As we have equalized all components beforehand, we only use the poles and 222 zeros of the U-component for this step. Continuity of the signal is forced at the beginning and at the 223 end of the glitch window by Lagrangian multipliers. The signal is then considered a glitch when the 224 variance residual after glitch removal is less than 1-2 % of the original data squared energy over a 225 running window of 50 s, starting 5 s before the glitch maximum. To remove the glitch spikes after the 226 glitch removal, a delta impulse is then searched around the glitch time and removed if associated with 227 a 50 % variance reduction of the signal in a window of width ±3 s. Glitches and spikes amplitudes 228 are inverted on the three axes. We use these amplitudes to calculate dip, azimuth and amplitudes 229 of the spikes that we use to potentially located glitch source (Section 6.1). An average of about 170 230 glitches per sol is found for 1 % of variance residual and about 100 glitches per sol for 0.5 % of variance 231 residual. For the former case, about 40 % are detected on the three components while the other are on 232 single VBB components. As this approach is detecting the glitch through the success of the functions' 233 fit with data, glitch removal is a sub-product of the method. were made by UCLA and IPGP, and 140 by MPS and ISAE, however, the latter two detected less 240 glitches during the noise daytime. Figure 2a shows the 73 glitches that were common to all 4 groups,

241
-5-manuscript submitted to Earth and Space Science which correspond to those with the largest amplitude. Table 1 shows the number of detected glitches 242 common to pairs of groups. The non-common glitches are plotted color-coded according to each 243 group. An expanded section (Fig. 2b) reveals that the various criteria detect mutually exclusive 244 glitches as the noise level is approached. We note that the Marsquake Service (MQS, Clinton et al.,245 2018) continuously monitors InSight's seismic data to detect and catalogue seismic events (InSight 246 Marsquake Service, 2020). As part of their routine they manually seek and annotate glitches with 247 principal focus on time windows of seismic events. Our detection methods generally compare well with 248 these manual annotations both in amount and onsets of glitches, especially for larger ones. For smaller 249 annotated glitches, i. e. less than 1e −8 ms −1 in amplitude, we find that each detection method, if the 250 parameters are chosen sensitive enough, delivers satisfying results with the amount of false detections 251 only slightly increased. However, not each annotated glitch is detected as the noise level is approached 252 and the signal-to-noise ratio hence decreases. Nevertheless, our comparisons show that our algorithms 253 for glitch detection are reliable in most circumstances. Our working hypothesis is that glitches in SEIS' time series data represent sudden steps in the 256 sensed acceleration convolved with the instrument response of the respective seismometer, either VBB 257 or SP. We can use that assumption to constrain the physical mechanism that led to the glitch. When SEIS sensor assembly (including the leveling system) will therefore point in the horizontal direction.

267
This is true for both SP and VBB. Any other direction cannot be explained by a rigid motion of SEIS 268 and must be due to instrumental artifacts.

269
It is useful to recall the sign convention for accelerometers: a positive output signal corresponds  The determination of the apparent glitch azimuth and incidence angles is implemented in our  The multi-component glitches for VBB and SP are illustrated in Figure 4. Especially for VBB, 342 for which we generally detect more glitches, clear patterns emerge over the period of 2019. We discuss 343 five of these patterns in the following. 344 We observe a glitch pattern with associated acceleration change pointing towards North (blue  to the West. This configuration produces colder temperatures on the east side during the night than 401 on the west side (and the opposite during the day), with larger gradients between IF1-IF2 or IF1-IF3 402 than between IF2-IF3. This is likely the primary source of these thermal glitch patterns. We note  SP that match in azimuth were discussed already in the previous Section 3.2 (green and blue dots, 409 parts of patterns 1 and 2). As we pointed out, these glitches show incidence angles of ∼ 90 • for both 410 VBB and SP and therefore could signify the whole SEIS instrument tilting.

411
The most prominent glitch pattern in Figure 5 is the one at azimuths of ∼ times but now with an average rate of less than 2 minutes per sol. After the conjunction period, during 420 which the heaters were turned off, we observe the same as for many other glitch patterns; a more diffuse 421 signature of the glitch azimuths and incidence angles that seem to return to pre-conjunction states 422 only ∼100 sols later. Also, the onsets time now drift towards later times (red to yellow) each sol which 423 interestingly coincides with the fact that the Martian solstice occurred just after the conjunction on 424 Sol 308. For this pattern as a whole, we were able to clearly identify the critical temperature around 425 which the glitches occur. As Figure 5d,f demonstrates, the glitch onset times strikingly follow the 426 iso-temperature curve at −54 • C for both VBB and SP. In addition for VBB, there are more patterns 427 with similar behaviour for which we could find the critical temperatures; these correspond to pattern 428 3 (red and pink dots, Section 3.2). All this evidence once more supports the fact that most glitches 429 are thermally caused. Note that the temperature sensor we used here is scientific temperature sensor permanent (ever-lasting) steps in acceleration and displacement, respectively, all methods prefer to 441 correct the raw data rather than the data after conversion to physical units.

442
-9-manuscript submitted to Earth and Space Science The MPS group models a glitch waveform for each detected glitch using three parameters: an 443 amplitude scaling factor, an offset, and a linear trend parameter. An example of its glitch removal is shown in Figure 6.

465
The UCLA group carries out glitch and spike removal on 10/20 SPS data. Some glitches show 466 symmetric or asymmetric broadening relative to the glitch template, suggesting the source function 467 is more complicated than a Heaviside step in acceleration. As a first approximation, convolution with 468 a unit Gaussian or exponential decay, which adds an extra parameter, significantly improves the fit, 469 but runs the risk of over-fitting data. To minimize this effect, the approach is only applied to data 470 that show >0.9 correlation coefficient with the glitch corresponding to our acceleration step-model.

471
Glitch (sometimes broadened) and spike templates were fit to the glitches and spikes, respectively, 472 using NLSQ. Because of the delta-like shape of the spike over one or two sample intervals, the starting 473 model must find the location to within a fraction of a sample interval (e.g. 0.05 s). Glitches are easier 474 to fit than spikes, being low frequency, and requirements on the starting model are less stringent.

475
Spikes are much smaller in 2 SPS data relative to glitch sizes. Thus 2 SPS data were used to generate 476 a glitch catalog (Section 2). The starting parameters from the 2 SPS fits were then used to fit glitches 477 in the 20 SPS data and residuals were calculated. The residuals were examined for the presence of 478 a spike in the data before the glitch peak, by requiring its amplitude to be greater than 5 standard 479 deviations of the residuals after the peak. If true, an iterative forward model was run by shifting The removal algorithm of the ISAE group is basically described in Section 2.2 (glitch detection).

488
Once a glitch has been detected using cross-correlations between the model and data, the model 489 without linear trend and offset is subtracted from the data. This method is implemented for all 490 sampling frequencies available. Spike removal and deviations from the simplified acceleration step-491 model are not implemented.

492
The IPGP group inverts three consecutive acceleration step sources for the glitch which allows 493 not only to invert for multi-component glitches occurring within these 3 samples but also to invert 494 for the phase delay through finite-difference approximation of the first and second time derivative.

495
This linear approach allows the inversion to provide identical results in the U, V, W coordinates or 496 in the Z, N, E coordinates, as the rotation between the two coordinates systems is a linear relation.
the coordinate systems.

500
In the end, all the proposed deglitching methods are nevertheless based on the same idea of 501 assuming a step in acceleration and displacement to model a glitch and spike, respectively, by using 502 the instrument impulse response of either the VBB or SP seismometer. Removal differences across 503 the methods are mostly due to thresholds below which a glitch is removed or not, and by how these 504 methods attempt to fit glitches that do not fully correspond to our acceleration step-model. No general 505 rule on the thresholds can be provided as they depend on the data processing target. it is also possibly to filter them out rather than removing them from the raw data, however, small 518 artefacts depending on the exact case may remain. All these arguments combined is the reason 519 we do not provide glitch and/or spike corrected data for all available periods but instead make our 520 codes available, enabling own comparisons and removal choices to those interested. An example of 521 glitch removal showing all four methods is demonstrated in Figure 7 for two glitches occurring during 522 marsquake S0173a. 523 We lastly point out that we have discontinued our deglitching efforts using the stationary wavelet 524 transform as described in the Supplement V of Lognonné et al. (2020). Whilst this approach provided 525 promising and correct results for a fair amount of cases (as far as one can tell), there is no underlying, 526 physical model involved and the implicit data 'correction' therefore seemed too arbitrary. For many 527 cases this approach further introduced DC-offsets in the deglitched data whose amplitudes and lengths 528 depended on the length of data read (and therefore maximum decomposition level); an artifact that 529 we could never manage to fully avoid. Throughout this paper we have assumed that glitches can be understood as steps in acceleration 532 and glitch spikes as steps in displacement. This model allowed us to successfully detect, analyse 533 and remove one-and multi-component glitches for both VBB and SP. In the following we detail the 534 theoretical considerations behind this simple model.

535
Let us assume glitches are caused by a small instantaneous tilt. By instantaneous we mean that 536 the time history of the tilting is so short that it cannot be resolved with any given sampling frequency 537 available to us (maximum 100 sps). We are thus allowed to idealize any step in time by a Heaviside 538 function. Physically such short instantaneous events can for example be the result of stick-slip events.

539
The small tilt is assumed to be the result of a rotation around a horizontal axis, a. Recall that 540 the VBB is a pendulum seismometer where the (inverted) pendulum is constrained to rotate around 541 a horizontal axis, b. The sensitive direction, s, of the pendulum is perpendicular to the b axis and is 542 inclined relative to the horizontal plane by a dip angle of δ = −29.3 • . Let us also assume for simplicity 543 that all the mass of the pendulum is concentrated in its center of gravity (CoG) -which would be the 544 case for a mathematical pendulum.

545
Now we can distinguish five cases which differ by the location of the accelerometer relative to the 546 tilt axis, a:

547
(1) the two axes a and b are parallel and a passes through CoG: in this case the accelerometer 548 gets only reoriented relative to the gravity vector but the CoG stays in place.

549
-11-manuscript submitted to Earth and Space Science (2) the two axes are parallel and a does not pass through CoG but is at the same height as the 550 CoG: in this case the accelerometer gets displaced vertically and reoriented relative to the gravity 551 vector. However this reorientation is negligible because it is only a second order effect.

552
(3) the two axes are parallel and a does not pass through CoG. Furthermore a line parallel to s 553 passing through CoG intersects with a. In this case the accelerometer gets displaced vertically and 554 reoriented. However the displacement is in the direction perpendicular to the sensitive axis and hence 555 is not seen by the accelerometer. Only the reorientation is sensed.

556
(4) For all other locations of the rotation axis a for which a and b are parallel the accelerometer 557 will see both a displacement and a reorientation relative to the gravity vector.

558
(5) For the general case where a and b are not parallel the same arguments can be made but the 559 effect sensed for a given tilt angle will always be reduced relative to the case with parallel axes a and 560 b since the tilting is reduced.

561
As soon as the accelerometer gets reoriented relative to the gravity vector we expect to see the 562 response due to a step in acceleration, because the projection of the gravity vector into the sensitive 563 direction is changed. In those cases where the accelerometer gets displaced we expect to see the 564 response due to a step in displacement. The five cases then only differ in the relative size of the 565 displacement and tilting.

566
What do these signals look like? In Figure 6 we have plotted the response of the VBB sensors can use the modelled glitch and spike to remove them from the data.

573
Can these signals explain the data? As Figure 6 also demonstrates, the modeled responses have 574 been shifted in time and scaled to match the data. The fit is excellent both for the low-frequency 575 glitch and the high-frequency spike. We take this as confirmation that our simple model is capable In the following we briefly discuss other aspects of glitches and spikes that we encountered during 590 our investigations. This section shall therefore complement our understanding of glitches and detail 591 some more implications.

601
We illustrate this geometry with the glitch example of Figure 6 and recall the glitch and spike 602 characteristics in Table 2 During the night, very small but also large rotation radii are found, likely resulting from internal 612 deformation of the Evacuated Container triggered by thermal effects, as discussed previously. During 613 the day however, the rotation radii of the glitches are more stable and in the range 10-30 cm, suggesting 614 an external source and therefore rigid tilt of SEIS, likely generated by the atmospheric activity.  We also observe several glitches, circled in red, that happen at the same time as the IDA motions.

634
One of the tell-tale signs of a glitch is when we observe an offset in acceleration in the seismic 635 components. We interestingly observe that the BHE-component shows steps of the same sign for 636 both the arm loading and unloading. Two of the glitches further appear to involve the whole sensor 637 assembly as they are seen on both the VBB and SP. Other glitches seem to be limited to one or more 638 components of the VBB. This all points towards that these glitches are internally caused and only 639 triggered by the IDA movement. Attempting to remove these IDA-induced glitches show convincing 640 fits with our acceleration step-model for the BHV and BHW components, however, for the BHU panels). Nevertheless, IDA movements are limited and therefore this type of glitch does not represent 643 a major contamination of InSight's seismic data. rests undergo a large daily temperature cycle.

672
The ∼80 K peak-to-peak ambient daily temperature variations are attenuated by the different 673 thermal shields but still reach ∼15 K inside the evacuated titanium sphere hosting the three VBB 674 sensors. These temperature fluctuations inevitably lead to thermal strains and thermally induced 675 stresses at the contacts between materials with different thermal expansion coefficients. These stresses 676 will in turn lead to additional elastic deformations. Alternatively, these stresses can be relaxed by 677 a variety of irreversible mechanisms such as creep, diffusion of lattice dislocations or stick-slip along 678 mechanical contacts. While we do not know which actual stress relaxation mechanism or which 679 combination of mechanisms is at play, we attribute thermally related glitches to intermittent stress 680 relaxation events such as for example stick-slip events.

681
The question of whether external events can trigger glitches arises when we inspect marsquake 682 S0173a (Fig. 7), the VBB response to certain pressure drops, or the VBB response to ground loading 683 experiments with the scoop of the instrument deployment arm (IDA, see Fig. 8). In all these cases 684 the seismic waveforms are contaminated by a glitch. We argue that external events alone do not cause over the sol, one may suspect that triggered glitches occur already within a few seconds after an 697 arrival if following our model. We found no obvious relation (Fig. SI2-3). Whilst the number of 698 events with clear P and S arrivals is small, and a more thorough re-analysis with a larger data set may 699 be worthwhile, all our analyses combined still suggest that the timing of glitches generally has a strong 700 stochastic component next to a deterministic component. This is further supported by the frequency-701 amplitude distributions of glitches per component that seemingly follow a Gutenberg-Richter relation 702 (Fig. SI2-3), and the presence of the diurnal harmonic and all its integer multiples in a time series 703 composed of modeled glitches (Fig. SI2-4c). either North or South, that is, either close to the LSA-tether system or diametrically opposed (Fig. 4, 721 patterns 1 and 2). Whilst the picture is not fully conclusive (Section 3.2), there remains the suspicion 722 that the LSA-tether system or even the lander exert influence on SEIS and therefore promote glitch 723 production via mechanisms for which we have no unique interpretations.

724
Lastly, we mention that glitch spikes seem to largely coincide with "donks", yet another type 725 of data disturbance typically only visible on VBB and SP seismic data of 20 SPS and higher. The 726 relationship between donks and glitch spikes was not analysed within the scope of this paper but will 727 be more detailed in different publications related to non-seismic signals observed on SEIS. i.e., ∼15 K compared to a few mK. Given the harsh environments typically found on extra-terrestrial, 736 planetary bodies, it may not be easy to achieve higher thermal stability however it should be considered 737 by engineers. We can only speculate as to the exact sources of glitch production within the instrument.

738
While we have good candidates (see further above), the fact remains that InSight's seismometers, 739 especially the VBB, are complex devices consisting of many materials, joints and connections. One 740 way to approach thermal glitch reduction may therefore be to use fewer materials and thus minimise 741 potential thermal conductivity gradients, stresses and expansions. A last, ultimate step to achieve 742 thermal stability would be to completely bury the instrument and possibly even the tether but this may high-frequency spikes that occur simultaneously with the glitch onsets (Fig. 1). In this model, glitches shunt assembly / tether pushing and pulling on the SEIS instrument. We illustrate the two cases of 763 most common glitch production in Figure 9.

764
Whilst terrestrial data influenced by glitches may simply be discarded due to their difficult han-765 dling, this represents no valid option for the seismic data returned from Mars. We therefore devoted 766 much of our efforts to develop code for the glitch and spike removal (Section 4). Our algorithms have 767 proven successful in many cases for both seismometers VBB and SP (Figs. 6 and 7). Of course, 768 there remain glitches and spikes especially of smaller amplitudes that we cannot sufficiently well fit 769 and therefore confidently remove. To account for such glitches nevertheless, we have slightly deviated 770 from our step-model in acceleration to improve on their removal, i.e., we introduced fits for non-zero 771 rise times (MPS), for a combination of multiple source-functions (UCLA), and for three consecutive 772 acceleration steps of varying amplitudes (IPGP). The resulting glitch models of these adaptations still 773 produce glitch waveforms close to the ones corresponding to a zero-rise time acceleration step, allowing 774 however to fit for glitches whose responses are broader than the ones corresponding to our simplified 775 step model. As we demonstrate in Figure 10 for VBB long-period spectra to look for Phobos' tides 776 and for receiver functions of the marsquake S0173a, removing glitches following the approaches pre-777 sented here indeed allows to improve on the quality of seismic data and may hence help to accomplish 778 InSight's scientific goals.

779
As no glitch removal algorithm can warrant a perfect clean-up of all glitches and their spikes, we 780 prefer to not provide a deglitched time series of all available data. Instead, we have assembled our al-781 gorithms for glitch detection, glitch polarization analysis, and glitch removal into one Python / ObsPy 782 toolbox. Some convenient functions for data retrieval and handling are also implemented. The package 783 further holds MATLAB scripts to perform glitch detection and removal tasks as presented. Its link 784 is: https://pss-gitlab.math.univ-paris-diderot.fr/data-processing-wg/seisglitch. Docu-785 mentation is available. Together with this code we also provide deglitched data for a selection of   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, demonstrating that multi-component SP glitches occur only among the horizontal SP V and SP W components. Color code is same as in Fig. 3.
-23-manuscript submitted to Earth and Space Science Figure 5. a,d) Glitches in 2019 that occurred simultaneously on VBB and SP. Glitch azimuths agree for patterns 1 and 2 (blue and green dots, compare Fig. 4) but not for pattern 5. Color code is same as in Fig.   3; b,e) example of our polarization analysis of the same glitch for VBB and SP on 2019-07-24T18:50:01 (Sol 234). The azimuths and incidence angles for this glitch are almost identical on VBB and SP. c,f) normalised glitch amplitudes as a function of sols over local mean solar time (LMST; different detection method than in sub-plots a-d). Note how the iso-temperature curve at −54 • C (scientific temperature sensor A, channel 03.VKI) matches the glitches corresponding to pattern 5, thus supporting thermal causes for glitches of this pattern. Figure   Gray vertical lines: theoretical onsets identical for glitch and spike; a: calculated amplitudes of glitches and spikes; t: time difference between calculated glitch and spike onsets smaller than sampling period (sub-sample fitting); VR: achieved variance reduction.
-25-manuscript submitted to Earth and Space Science Figure 7. Comparison of VBB raw data at 20 SPS with the corrected data according to our four deglitching methods. The ISAE method does not correct for glitch spikes. The IPGP method only processes 2 SPS data.
Linear trends were removed for plotting purposes. The data show marsquake S0173a on 2019-05-23T02:23 (Sol 173), one of the best-quality low frequency events identified to date by the Marsquake Service (MQS, Clinton et al., 2018, catalog: InSight Marsquake Service, 2020. Vertical purple lines; P-and S-phases as identified by MQS; vertical black lines: glitches as annotated by MQS (Section 2.5). Clearly visible right after the P-phase onset is a prominent glitch. In the reconstructed ZNE-data this glitch is almost only present on the horizontal components (AZ=330 • , INC=99 • ). 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 the removal algorithms (Section 4). . We suspect such effects to be the primary reason for thermally-caused multi-component glitches such as shown in patterns 3-5 (Fig. 4). b) SEIS tilt α, corresponding to a true, rigid motion of the whole instrument. Our analysis suggests that the minority of glitches, e.g. patterns 1-2 (Fig. 4), are caused by this scenario. Note that in both cases the VBB sensors may experience a tilt and a displacement (Sections 5 and 6.1). Similar considerations apply for the SP sensors (not shown) that are mounted on the leveling system (SEIS feet) support structure (Fayon et al., 2018). This support structure is connected to the Evacuated Container containing the VBB sensors via three mounting bolts (Sections 3.2 and 6.1). The heaters are mounted to the support structure, too (not shown, Section 3.1).
For an accurate illustration of the SEIS sensor assembly, see Lognonné et al. (2019). Green lines: moving pendulum parts; P: proof mass; δ: VBB sensor dip ∼ −30 • . The tilt α is here depicted as 10 • for both cases but is in reality in the order of nano-radiant. The deglitched data (DG, ISAE method) after temperature decorrelation show reduced spectral peaks that are caused by the 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 (left) and deglitched data (right, UCLA method) and their Ps-receiver functions for marsquake S0173a. Top panels: waveform data around P-wave onset of S0173a, band-pass filtered between 0.1-0.8 Hz where most of the signal energy is located, and rotated into radial and transverse directions. Note the prominent glitch around 20 s that is still dominating the horizontal components after filtering. Gray boxes: time window used for the deconvolution in Ps-receiver function calculation shown in lower panels: the long-period contamination by the glitch becomes apparent after 8 s on the horizontal components, masks any later arrivals, and also casts doubts on the reliability of earlier phases. For example, an additional arrival near 7.3 s is now clearly visible on the radial component, a phase that is also observed in receiver functions for other marsquakes that are not contaminated by glitches (Lognonné et al., 2020, Supplement IV). with the following linear system of equations: where A represents the base transformation matrix, δ i the sensor dip of sensor i, and φ i the sensor 17 azimuth of sensor i clockwise from N. Note that sensor dips are defined as positive downwards from 18 the horizontal plane (e.g. , which is taken into account in A. To reconstruct data 19 recorded in the UVW-system into the ZNE-system, we must use the inverse operation: with A −1 the inverse matrix of A. If we now consider a glitch that occurred only on VBB U with an 21 amplitude U = 1 (V = 0, W = 0), insert those values into Equation 2, and use the following equations 22 to determine the apparent glitch azimuth defined clock-wise from N, AZ, and apparent glitch incidence 23 INC defined as the angle with respect to the Z-axis, it follows: We can calculate the inverse matrix elements (A −1 ) j1 with the known VBB sensor azimuths φ U = Thus, the apparent azimuth and incidence angles of a one-component VBB glitch will not point 28 in the direction of the sensitive direction of the affected VBB sensor. Instead, the polarization vector 29 is parallel to the vector cross-product of the remaining two components that do not show the glitch.

30
Due to the similar arrangement of all VBB's sensors (see Fig. 1a  and dips of δ U = −89.9 • , δ V = 0.0 • and δ U = 0.0 • (Fig. 1c in main paper), one finds that for SP U

42
(Z) the azimuth and incidence angles will follow one's intuition closely and be 0 • and 0 • , respectively. oriented in the horizontal plane.

53
The message from these theoretical considerations is that our glitch polarization analysis will 54 deliver azimuths and incidence angles that correctly account for the non-orthogonality of VBB and 55 SP; the vectors defined by these angles point into the only physically possible directions for a given 56 one-, two-or three-component glitch, assuming a rigid motion of SEIS. On the other hand, for the 57 interpretations of these angles, it must be born in mind that VBB incidence angles may carry counter-58 intuitive information whilst SP azimuth angles for one-component glitches will not align with the 59 respective sensor azimuths but diverge by ∼ 30 • .

60
At this stage we also note that whilst the poles and zeros of the VBB and SP seismometer 61 responses are well determined, the same does not apply fully for the generator constants (gains). In 62 the worst case they may differ up to 10% from the absolute values known by pre-mission tests. To 63 convince ourselves of the correctness of determined glitch azimuths and incidences with respect to 64 these constants we conducted a test: we took the raw data of one-and multi-component glitches of 65 different amplitudes and divided the respective components by their gains that we allowed to vary 66 each by up to ±10%. For each permutation, we then rotated into the ZNE-system and performed 67 the polarization analysis. For VBB, we find that glitch azimuths and incidences generally stay within 68 ±5 • and ±4 • , respectively. For SP, we find that glitch azimuths and incidences generally stay within 69 ±3 • and ±1 • , respectively, the latter of which is because SP multi-component glitches occur only 70 on the horizontal components. All these values are smaller than the typical errors of polarization 71 measurements and we can therefore assume the resulting glitch patterns to be reliable.

72
Let us consider a general geometry such as depicted in Figure 9 in the main paper where a cross 74 section through a VBB sensor perpendicular to its hinge is graphed. In this figure, the SEIS sensor 75 assembly is rotated around the tip of leg A by a small angle α such that the tip of leg B is raised by d·α, 76 with d being the distance between the tips of the legs. The sensitive axis of the VBB accelerometer, 77 denoted with the unit vectorσ, is inclined relative to the horizontal by the angle δ which is close to 78 -29 • , depending on the VBB sensor.

79
The force of gravity acting on the proof mass M and which the suspension spring has to counterbalance is: where g = 3.71m/s 2 is the surface gravity on Mars. After the tilting of SEIS by the angle α, the 80 projection of g onto the sensitive axes changes and it follows: The change in accelerationü produced by the tilting thus is: Since the rotation axis does not go through the center of gravity P of the proof mass M , the rotation 83 leads also to a displacement of the proof mass. In our case this displacement, y, is a small arc segment 84 of a circle with radius r = AP around the tip of leg A: y = r · α. The accelerometer only senses the 85 projection of this displacement onto its sensitive direction. If we define the unit vectorr as: the sensed displacement then becomes: 87 u = r · α · |r ×σ|.
What is the time history of this tilt and the simultaneous displacement? As we shall see, the 88 data can be very well modeled by assuming that the time dependence follows a Heavyside function, 89 that is the tilt and the displacements occur over a time interval much shorter than can be resolved 90 with the given sampling interval. In the analyzed glitches we see little to no indication for a slowly 91 progressing tilt.

92
Now we have to account for the fact that inertial accelerometers like the VBB and SP seismometers 93 in the SEIS package have a frequency dependent sensitivity to ground motion. This is described by The summed output U from the acceleration step due to the tilting at time t o and the associated 100 displacement step then becomes:

101
In this section, we provide some additional figures we have created while investigating the glitch 121 plus spike phenomenon. We will not put each figure into context but would simply like to refer to 122 their captions for understanding.  To investigate a possible triggering of glitches by seismic arrivals, we compare detected glitches with low-frequency and broadband events of qualities A-C ('A' is best quality). a) All detected glitches within one hour after the P-arrival, or the beginning of the visible signal where no clear arrival could be identified. Events with arrivals are sorted by S-P time, others by sol. Blue: P arrivals, red: S arrivals, horizontal lines: time windows of visible quake signal, stars: glitches. b) Time between glitch and the last preceding arrival (P or S). Stars: Glitches, Histogram: number of glitches in 5 min time windows. Only 6 of 72 considered glitches occur within 10 min after the last arrival. Given this small number, we do not consider the difference between the first and the second bin as significant, indicating that glitches during seismic events are not occurring significantly more than during periods of no seismic events.  Figure SI2-3: 2019 VBB glitch histograms per component, detected by the MPS method with more sensitive glitch detector settings than utilised in the main paper. We find a seemingly stable Gutenberg-Richter relation with b-values of ∼1-1.3, and roll-off glitch amplitudes of ∼1e-8 m/s (RAW data corrected for gain. The velocity response is flat for periods shorter than 16 seconds). This may indicate an underlying stochastic process behind the glitch production that, perhaps, points once more to thermal causes of glitches. Figure SI2-4: a) Cumulative contribution of glitches to the total acceleration signal. The glitches have been sorted by their variance reduction obtained from the glitch modeling. This panel shows that poorly modeled glitches (variance reduction of less than e.g. 85%) make up only a small fraction of the total acceleration signal: 25%, 25% and 18% for U, V and W respectively. b) Glitches sorted by variance reduction: For the chosen sensitivity of the MPS detector in the main paper and for the time interval Sol 70 through Sol 260, there are 13000 glitches with variance reduction less than 85% and 18000 with variance reduction greater than 85%. Taken together panels (a) and (b) show that the largest contribution in terms of signal amplitude comes from the large and well modeled glitches. In terms of signal power the contribution of the large and well modeled glitches becomes even more dominant. c) Contribution of all modeled glitches to the acceleration background for the three VBB components U, V and W. All glitches for Sols 70 through 260 for which the variance reduction in the glitch modeling stage exceeded 85% are included. A glitch corresponds to a step in acceleration at a particular time. Here we have added up in the time domain 18000 step functions, one for each glitch, with the step size corresponding to the glitch amplitude. The power spectral density of the resulting stair case like, noise free time series has been analyzed. The harmonics at integer multiples of 1 cycle/Sol are a strong indication that the glitches have a thermal cause. This analysis is a complementary method to quantify the contribution of glitches to the VBB analysis presented in figure 10a of the main paper.