A Comodulation Analysis of Atmospheric Energy Injection Into the Ground Motion at InSight, Mars

Seismic observations involve signals that can be easily masked by noise injection. For the NASA Mars lander InSight, the atmosphere is a significant noise contributor, impeding the identification of seismic events for two‐thirds of a Martian day. While the noise is below that seen at even the quietest sites on Earth, the amplitude of seismic signals on Mars is also considerably lower, requiring an understanding and quantification of environmental injection at unprecedented levels. Mars’ ground and atmosphere are a continuously coupled seismic system, and although atmospheric functions are of distinct origins, the superposition of these noise contributions is poorly understood, making separation a challenging task. We present a novel method for partitioning the observed signal into seismic and environmental contributions. Atmospheric pressure and wind fluctuations are shown to exhibit temporal cross‐frequency coupling across multiple bands, injecting noise that is neither random nor coherent. We investigate this through comodulation, quantifying the synchrony of the seismic motion, wind and pressure signals. By working in the time‐frequency domain, we discriminate between the different origins of underlying processes and determine the site's environmental sensitivity. Our method aims to create a virtual vault at InSight's landing site on Mars, shielding the seismometers with effective postprocessing in lieu of a physical vault. This allows us to describe the environmental and seismic signals over a sequence of sols, to quantify the wind and pressure injection and estimate the seismic content of possible marsquakes with a signal‐to‐noise ratio that can be quantified in terms of environmental independence. Finally, we exploit the relationship between the comodulated signals to identify their sources.


Introduction
The InSight seismic station on Mars is the first geophysical observatory ever deployed directly on the surface of another planet . It has recorded the quietest seismic traces observed on any solar-system body at periods shorter than 100 s since the Apollo mission owing to the lack of anthropogenic noise and oceanic microseismic noise on Mars . Unlike the Moon, strong motion due to local atmospheric injections dominates the data for much of the Martian day (sol). This was previously observed by the Viking lander seismic experiment (Anderson et al., 1977;Lorenz et al., 2017) where the ondeck seismometer was primarily affected by wind coupling. Prior to deployment of the seismometer, InSight recorded the impact of site noise while on-deck (Panning et al., 2020).
In light of the lessons learned from the Viking mission, InSight invested much effort in the deployment from the lander deck to the Martian surface of the Seismic Experiment for Interior Structure (SEIS) seismometer. As the mission's critical scientific instrument, SEIS was developed with the prime focus of recording ground motion caused by seismic events, offering a unique and unprecedented opportunity to identify marsquakes or meteoroid impacts. The SEIS assembly is comprised of two independent three-axis seismometers: the oblique Very Broad Band (VBB) and Short Period (SP) seismometers (Lognonné et al., 2019). SEIS measurements are complemented with the on-deck mounted Auxiliary Payload Sensor Suite (APSS, see Banfield et al., 2018;Spiga et al., 2018), which monitors environmental conditions to help quantify the impact of local atmospheric fluctuations on SEIS over both short and long timescales. APSS includes pressure, wind, and temperature sensors (TWINS). Most importantly, a robotic arm (Trebi-Ollennu et al., 2018) was included to move the SEIS assembly from the lander deck to a selected site on the ground and then cover it with the Wind and Thermal Shield (WTS). The WTS provides SEIS with at least some insulation from the wind and thermal fluctuations. As such, InSight has been able to operate as a seismic observatory on Mars and has detected several regional and distant marsquake events .
Although this careful deployment and subsequent shielding has proved significant in reducing the noise level for InSight, it has by no means completely isolated SEIS from the environment. During quiet windows, the background noise observed by SEIS has been below that seen at the quietest sites on Earth over much of its observational bandwidth. However, the seismic traces still exhibit clear correlation with environmental conditions dominated during two-thirds of each sol (mainly during the daytime) by signals attributed to atmospheric injection. The analysis and subsequent removal of this atmospheric injection is therefore required. Lander-induced vibrations, that predominantly originate from variations observed in the atmosphere, induce regolith displacements through elastic deformations and continuously couple to the seismic signal. This was predicted in prelanding assessments (Lognonné & Mosser, 1993;Mimoun et al., 2017;Murdoch et al., 2017). The observed relationships are complex and exhibit nonlinear coupling dynamics, complicating investigation of any seismic signals. In this work, we introduce a comodulation framework to investigate the continuous coupling between the atmosphere and ground acceleration, and quantify the noise contribution to the seismic signal. In particular, the observed seismic signal is modulated by injection from the atmosphere as well as vibrations transmitted through the planet, for example, by tectonic activity and meteoroid impacts. Our method is applied to the detected marsquakes to quantify their independence from environmental contamination and help discriminate the origins of observed seismic features. Finally, we show that the comodulation approach can be applied to retrieve wind speed and direction, and also the pressure fluctuations over short and long timescales. Retrieving wind speed, wind direction and pressure fluctuations from the seismic data alone could become increasingly important for the InSight mission as the solar cells begin to reduce in power output  and the APSS instrument has to be shut down to save power.

Background
Much effort on Earth has gone into understanding and minimizing the background noise in seismic observations. While this is bounded by instrument sensitivity, the installation most often provides the limiting factor in terms of site-noise injections across the recording bandwidth (Ringler et al., 2020). A significant driver of this noise is the sensitivity to the environmental parameters of temperature, pressure and wind, either through first-order forcing or second-order coupling . Such coupling could arise, for example, from the transferred energy of wind-induced motion on local trees or structures to the shallow subsurface (Johnson, Meng, et al., 2019). These have linked and overlapping contributions, and this study focuses on the pressure and wind sources.
Different mechanisms cause environmental injection into the seismic signal at different frequencies. Low frequency injections, with periods above 10 s, are relatively well understood physically; for the vertical component, the pressure dominates (Ziolkowski, 1973;Zürn & Widmer, 1995) while the horizontal component sees tilt from both pressure and wind effects (De Angelis & Bodin, 2012;Sorrells, 1971;Sorrells & Goforth, 1973;Sorrells et al., 1971). This direct forcing is accentuated by local site characteristics including the ground material and vault/installation geometry which induce extra effects based on material response and forcing on nearby objects (Dybing et al., 2019;Johnson, Meng, et al., 2019;Wolin et al., 2015).
Site contributions extend to higher frequencies and recent studies of the wind-induced noise have performed statistical analyses across arrays to determine their spatial variation (Dybing et al., 2019;Hutt et al., 2017;Johnson, Meng, et al., 2019;Lott et al., 2017). Lack of local coherence between stations has been reported (Dybing et al., 2019), showing that turbulence causes a complicated set of features which can produce tremor-like waveforms that can be confused with low signal-to-noise-ratio (SNR) local events (Johnson, Meng, et al., 2019). This can include high-frequency injections above 1 Hz (Withers et al., 1996). A theoretical approach to high-frequency injections based on the shear stress acting at the site requires a thorough knowledge of the site's characteristics (Naderyan et al., 2016). The sensitivity across a broadband range (0.1-100 Hz) has been examined in Lott et al. (2017) which found that wind speed and direction are both important, driven by the local topography. In general, the intermediate seismic bandwidth (0.1-1 Hz) is dominated by the oceanic microseism and so environmental effects are rarely observed even in moderate-noise sites (Dybing et al., 2019). In contrast, for InSight, this is the prime observational bandwidth and so our analysis requires quantifying environmental injection at unprecedented levels compared to terrestrial seismology at these frequencies.
With respect to the InSight mission, a noise model was developed prior to landing. Its aim was to describe and predict the noise contributions from different sources including those excited by the wind and pressure . This model identifies that the major injections will consist of: 1. Instrument self-noise 2. Temperature/thermoelastic sensitivity 3. External magnetic fields 4. Ground tilt from both local effects (e.g., vortex passage ) and regional atmospheric effects (e.g., gravity waves ) 5. Lander motion The ground tilt and lander motion are analogous to the low-and high-frequency injection pathways, with a transition at about 0.1 Hz, as discussed in the Earth literature. At longer periods, this tilt is predominantly the response of the ground to direct forcing from the pressure field (Kenda et al., 2017Lognonné & Mosser, 1993) and so a coherency analysis can be performed leading to a decorrelation . Tilt can also be caused by wind forcing on the lander and the WTS. Lander injection was predicted to be a significant contributor  and as Mars lacks ocean-generated microseism, these injections will be observed across the full bandwidth. This injection is the analog of the coupling of wind from nearby structures and plants on Earth (Johnson, Meng, et al., 2019) and causes vertical motion in addition to the tilt seen in the horizontal components. This has been analyzed for InSight in terms of the drag and lift forces on the lander in Murdoch et al. (2017). Teanby et al. (2017) and Myhill et al. (2018) produced an analysis of how such forcing would couple through the regolith, using analog studies in Iceland to inform the contribution to the noise model. Furthermore, the lander has several resonances which are excited by the wind . As a result, the contributions are expected to be multiple and complex.
Understanding these atmospheric contributions is further complicated when only one seismic station is available; spurious signals cannot be rejected through local array correlations as is commonly performed on Earth. On the other hand, as there is lower background ambient seismicity on Mars, the noise injection can be examined across the full bandwidth. This means that approaches to tackle the wind and pressure injections should be separable from the seismic signal at lower amplitudes, potentially leading to better techniques for studies on Earth.
In this work, we introduce a novel approach in which we analyze the sensitivity of the seismometer response to injections with an environmental source across the full bandwidth. In doing so, we are able to determine the independence of seismic events and other observed features. This approach effectively creates a "virtual vault," where noise factors are attenuated by postprocessing of the seismic output rather than physically attenuated by the seismic installation. Conversely, the method can further exploit the quantified seismic and pressure sensitivity, and be applied in reverse to obtain estimates for the atmospheric variables of wind speed and wind direction.

Data
To help constrain the atmospheric injection into the seismic signal, InSight measures multiple environmental parameters continuously at a high sampling rate and sensitivity: wind speed and direction, temperature, pressure, and the vector magnetic field, using APSS . The TWINS sensors measure winds facing in opposite directions in the east and west directions, and stand at slightly different heights (of less than 10 cm) due to InSight's tilt in its landed position on Homestead hollow . A detailed study of the characteristics and meteorological phenomena observed from InSight's atmospheric science instrument data are described in Banfield et al. (2020).
Data are recorded continuously on board the lander at 1 sps (samples per second) for the wind and at 20 sps for pressure with the higher rate selected to capture transient excursions caused by passing vortices. Data are downlinked to Earth at selectable frequencies depending on the available bandwidth: 0.1 sps or 1 sps for the wind sensor and between 2 sps and 20 sps for the pressure sensor. Seismic data are recorded continuously and returned at 20 sps for the VBB. The SP records continuously at 100 sps, but also downlinks at 20 sps. Downlink requests for specific events of interest allow small time windows of pressure, wind and SP data to be transmitted at the higher recorded rates.
This study focuses on the continuous data (InSight Mars SEIS Data Service, 2019b). This allows for longer durations of uninterrupted data, permitting a deeper investigation into the atmospheric coupling measured by SEIS. Examples of some of the longest duration of uninterrupted 100 sps data are also investigated to demonstrate the broadband characteristics of atmospheric injection into the system and any variability in the seasonal evolution.

The Comodulation Framework
This work explores the injection of wind and pressure into the seismic signal from 0.01 to 50 Hz, the bandwidth of observed seismic events (InSight Marsquake Service, 2020). We adopt a data-driven analysis, the comodulation framework, to determine how the power in the environmental and seismic signals is related without needing to consider the fundamental physics driving such relationships.
This particularly builds upon the approaches discussed in Dybing et al. (2019), Johnson, Meng, et al. (2019), andLott et al. (2017) which use a statistical approach. Analysis of related signals through comodulation has notably been applied to understanding monaural and binaural processes of the human auditory system, where the perceptual comparison of sound by energy content reflects both the time and frequency information in the signal (Buus, 1985;J. W. Hall et al., 1984;Richards, 1987;Verhey et al., 2003). Further comodulation interactions were noted in electrical field oscillations in the brain during cortical cooperativity and spatial memory tasks (Bruns & Eckhorn, 2004;Shirvalkar et al., 2010). In this case, we are examining the power comodulation between signals from three different sensors: pressure, wind, and seismic.
To introduce this method, we first consider a qualitative investigation of sols 237-239 (at a solar longitude of L s = 59°), shown in Figure 1. This figure shows spectrograms for the three seismic axes, vertical (Z), north (N/S), and east (E/W), followed by the spectrogram of the pressure and the wind speed time series. The weather during this period was typical of northern spring conditions at the InSight landing site in Elysium CHARALAMBOUS ET AL.

10.1029/2020JE006538
Planitia . The ambient (large-scale) wind direction is from the south-east in the daytime with excursions from the south-west at night (see InSight weather observations in Banfield et al., 2020 andprelanding predictions in Spiga et al., 2018). Multiple dust-devil-like convective vortices occur through the daytime, identified by an abrupt drop in the atmospheric pressure with broadband injection and ground tilt seen on SEIS. The maximum pressure drop during this 3-sol period occurred on sol 239 with a magnitude of 4.7 Pa and can be observed in Figures 1a and 1b as a broadband impulse in the pressure and seismic spectrograms at 11.41 Local Mean Solar Time (LMST). This pressure drop is indicated by the star symbol in Figure 1b.
Three regimes can be identified in the meteorological diurnal cycle at the InSight landing site across the seasons : intense convective turbulence in the daytime between 08:00-17:00 LMST; a quiet period after sunset of low wind speed below the wind sensor's detection threshold, when Mars is quiet enough to reveal the instrument self-noise  and the majority of the identified seismic events are detected (17:00-19:00 LMST) ; and, following a transition in wind direction, a nighttime regime with gravity waves and small-amplitude shear-driven turbulence (19:00-07:00 LMST). Although transition times between these three regimes vary on a seasonal timescale, opportunistic windows for seismic detection emerge during the low noise periods. One example of such a window is indicated by the vertical gray dashed lines in Figure 1. The quiet periods can also be distinguished by the emergence of a continuously excited 2.4 Hz mode, predominantly on the vertical component with an origin that is as yet unknown . A 15-min high-frequency event exciting this 2.4 Hz resonance is observed in the final quiet window, at approximately 18:52 LMST on sol 239, with the time of occurrence indicated by a triangle in the third panel of Figure 1a.
The wind speed magnitude and power fluctuations observed in the spectrograms of pressure and ground acceleration in Figure 1 are strikingly similar across a wide range of frequencies. Despite this, the coherence between the raw time series is often low . This indicates a breakdown in the phase relationship between the signal amplitudes. However, the similarity of these spectrograms and the wind speed magnitude is reflected in the comodulation between the atmospheric and seismic signal powers.
The seismic signal power will be driven by the known effects previously outlined, for example, the effect of wind on local features, tilting, and vibrations. We can consider these effects by constructing a model where the measured atmospheric conditions, pressure, and wind are considered as inputs with the seismic ground acceleration as the output, that is, the response to forcing from various physical injection pathways. These include: 1. Direct forcing on the ground from the atmosphere 2. Vibrations of the lander coupled through the regolith 3. Forcing on the WTS coupled to the SEIS assembly through local ground deformations and direct forcing acting on the umbilical tether connecting the deployed instrumentation to the lander body The response of each of these pathways is complex and the proposed comodulation framework is invoked to capture the multiple pathways and describe the full effects of the system. Hence, we consider these pathways in terms of energy transfer, rather than signal amplitude. This energy injection is predominantly driven by dynamic pressure P , which is given by: where ρ is the air density and U is the wind speed, which injects energy through various paths that include the lander, WTS, and ground itself via the drag force F which is given by:

 F SCP
where S is the surface area and C the drag coefficient. These relationships highlight the many sensitivities of the system with the drag coefficients dependent on the Reynolds number which in turn is affected by the wind properties itself. The Reynolds number characterizes the flow which can be described as orderly laminar flow or chaotic turbulent flow, or an intermediate point between the two. Furthermore, the surface topography will also affect the wind properties. Although wind speed is the major variable, the relationship with pressure shown in Figure 1 indicates this cannot be ignored in these bandwidths as instrumental sensitivity to wind, turbulence, or the introduction of additional forcing pathways all vary with pressure.
This approach of relying on the signals themselves to quantify the injection is in contrast to the more theoretical studies in this issue Kenda et al., 2020;Murdoch et al., 2021) which can be applied within certain frequency bands to allow the extraction of near-surface elastic properties using the seismic injection from, for example, convective vortices. Studies on the Earth have often looked at the impact of only isolated variables, such as wind speed (Dybing et al., 2019), pressure (Zürn & Widmer, 1995), and temperature (Doody et al., 2017), considering that the interrelations between these variables may yield benefits. The comodulation presented in this work allows us to more generally quantify environmentally injected seismic noise across the full bandwidth.

Envelope Calculation and Comparison
The seismic acceleration and pressure spectrogram intensities are both correlated with the wind speed, as shown in Figure 1. To quantify this relationship, we examine the root mean square (RMS) envelopes, e[t], which yield an estimate of the signal power. In this section, we derive how these envelopes are obtained and then compared. For a given time series x[n] (n = 1, 2…, N), the spectrogram is calculated as a set of power-spectral-density (PSD) slices X[k, t] for a window length W len and overlap OL win , where k denotes the frequency bin and t the time step at which the window is taken. The envelope is then calculated as: where P and Q designate the frequency bins that equate to the desired bandwidth of study f min −f max .

Envelope Regression Through Moment Matching
The envelopes e Z [t], e N [t], e E [t], e P [t], and e W [t] are the signals that allow the correlation between the wind, pressure, and seismic measurements to be examined through a regression. We apply the method of moments as a statistical basis for this regression, a tool that has enjoyed a wide variety of applications as a common approach used for estimation, from econometrics (A. R. Hall, 2005;Hansen, 1982), to aerosol dynamics (McGraw, 1997).
The method of moments compares two time series in terms of their statistical moments, in this case with the first moment being the mean and the second central moment being the variance about the mean. In this approach, we assume one distribution to be the reference distribution, with the other being the empirical. The moments are therefore the primary parameters in determining a forecasting equation that transforms the empirical distribution, from its sample moments, to the corresponding population moments of the reference distribution (Hansen, 1982;Pearson, 1894 [t] up to second-order statistics, the mean of e Y is first set to zero and the variance is normalized to one. The variance is then scaled to that of e X and the mean added. Only one theoretical moment equation is required for each unknown parameter describing the time series. This is solved by the ensuing set of equations in terms of the moments arising from the time series of interest. An estimator is then obtained by substituting the moments evaluated from the sample to the theoretical moments of the reference time series. This process is summarized as: (2) where the mean operator is:  Figure 1c. The pressure and vertical ground acceleration envelopes are scaled and offset using the matched moments technique to obtain a regression to the wind speed allowing an estimate of the equivalent wind speed from the total power in each of these signals on the basis of wind-induced coupling. The correlation between signals in power as opposed to the usual amplitude is shown in Figure 1d, where the averaged correlation between the wind speed and seismic envelope is maintained throughout the sol. The correlation indicates synchronous fluctuations in the signal envelopes and therefore indicates the atmosphere is the primary noise injector masking the seismic signal.

Diurnal and Seasonal Sensitivity to Atmospheric Conditions
In this section, we describe the relationships between the seismic and atmospheric variables in terms of the comodulation at the daily and seasonal timescales. The analysis quantifies the distinct regimes and their relationships within a Martian sol. We achieve this through moment matching, allowing comparison of the degree of coupling between the atmospheric and seismic signals. These regimes are then analyzed over a range of bandwidths and the evolution of the relationships considered over the course of the mission. As the environmental components of pressure and wind are the driving factors, their relationship is then discussed. Finally, we show how the wind data can be predicted from the seismic and pressure signals using moment matching.

Diurnal Variation
We first examine at the diurnal timescale: 1. The ground acceleration response to the atmospheric variation over the course of a sol 2. The frequency dependence of this variation 3. The difference between the vertical and horizontal acceleration response for both 1 and 2 In order to explain the observations, we continue with the description of the relationships between the seismic and environmental variables over the diurnal cycle.

Time Dependence of Diurnal Ground Acceleration
Figure 2a plots the pressure and temperature variations for sols 237-239, with the vertical lines indicating sunrise and sunset in each sol. The wind speed is shown in Figure 2b. Figure 2c shows the vertical seismic acceleration envelope (in the 0.1-8 Hz band, calculated using windows of 50 s) against the wind speed for this 3-sol period. Three regimes of the wind-induced ground response are identified which coincide with the typical regimes of the diurnal cycle encountered at various seasons .
During the Martian sunset, a quiet period emerges when the wind drops below a threshold of 2.4 m s −1 , with the acceleration response confined within the red tail in the distribution of Figure 2c. We quantify this threshold as the point below which there is no correlation with the seismic signal, as we further discuss in Section 5.2. During this period, the acceleration response reaches the noise floor of the VBB sensor and so no relationship can be established between the wind speed and acceleration envelope. Impulse excursions are seen in the seismic acceleration envelope at stationary wind speeds, attributed to short-duration nonatmospheric sources, most often due to glitches (Scholz et al., 2020). The variance in this period is at its lowest for the ground acceleration, in line with this being the calmest daily interval. However, due to the unreliability of wind retrieval below this wind-speed threshold, the variance in the wind measurements is high.
Following the quiet period, in blue, the wind speed gradually rises from 3 to 6 m s −1 during the nighttime until dawn. This has been attributed to nocturnal low-level jet , with the ground acceleration appearing to vary much more strongly from low amplitude wind-speed variations. This is quantified by a steep power-law relationship with an exponent of a ≈ 4 ( Figure 2d). Generally, the nighttime regime is highly stable on Mars with a strong near-surface inversion, allowing only for shear-driven weak turbulence to develop to drive the wind. In contrast, the variance for acceleration is high due to increased sensitivity in this regime. This period is characterized by long-period fluctuations in the pressure, related to gravity waves and bores, while the temperature profile is stable as the surface is cooler (Figure 2a; Banfield et al., 2020).
The intense convective turbulence regime during daytime (dark gray) follows a shallower power law with an exponent of a ≈ 2 ( Figure 2d). The transition to this turbulent regime occurs at a mean ambient wind ∼6 m s −1 . The daytime regime is characterized by the highest wind-speed fluctuations in the Martian sol, which in turn generate the highest moment values in the distribution, for both mean and variance. Rapid fluctuations of both temperature and pressure produce the highest short-term variance in this period.
By separating these established power-law regimes, the relationships are quantifiable by a gain and offset and hence the moment-matching technique from Equation 2 can be applied in a piece-wise approach to provide a regression. In Figure 2b, we have transformed the acceleration domain to that of the wind speed, our reference. The moment-matching method captures the power-law variation by matching the logarithmic distribution in each regime using a single scaling and offset for the second-order moments.
In reality, these power laws do not abruptly shift but follow short transitions between each regime. The rapid rise of temperature during the Martian morning imposes steep positive gradients on both the atmospheric and surface temperature (Figure 2a), while fast cooling in the early Martian afternoon imposes an inverse temperature gradient, with such temperature swings also measured by previous missions Davy et al., 2010;Hess et al., 1977;Mueller et al., 2020;Schofield et al., 1997;Sutton et al., 1978). The transition from the nighttime to intense turbulence daytime regime takes place after sunrise, a period during which these steep atmospheric temperature gradients induce a distinct short-duration regime. Similarly, the reverse occurs before sunset and these effects could play a dominant role in the transitional behavior. A detailed breakdown of these regimes and their transitions with fitted power laws is illustrated in

Frequency Dependence of the Ground Acceleration Response
This analysis can be expanded to consider several bandwidths and include pressure, as shown in Figure 3. The envelopes for the vertical acceleration response and pressure are computed in 0.01-0.1 Hz, 0.1-0.5 Hz, 0.5-1 Hz, and 1-4 Hz bandwidths and plotted against the wind speed. The pressure envelope also correlates with the seismic-acceleration envelope. The relationships have lower dispersion than the wind speed and there are similar, although less distinct, separate regimes. At the lowest frequencies, 0.01-0.1 Hz, the correlation between the seismic acceleration envelope and pressure envelope/wind speed is much weaker, and restricted to higher signal powers. This may be caused by injections not captured by the environmental sensors, such as instrument glitches. This additional noncorrelated injection can be clearly seen recurring in the VBB Z spectrogram of Figure 3a below 0.1 Hz, with a weak energy injection during the quiet period and a stronger injection occurring from midnight and extending beyond sunrise. As the frequency increases, the instrumental noise raises the seismic signal threshold and a steeper power law represents a greater sensitivity of the seismic signal to environmental injection. In contrast, pressure and wind maintain a more consistent linear relationship, with an exponent of a ≈ 1 above 0.1 Hz. Below a threshold of 2.4 m s −1 , the response of both the pressure and seismic signals is very weak. As discussed in Section 4.1, the driving factor for the environmental injection is the dynamic pressure. While we have shown the pressure envelope also correlates, it is the wind speed that allows the direct determination of dynamic pressure and so we will limit our analysis of the seismic injection to the wind speed before returning to this aspect in Section 5.2 where we will discuss the comodulation factors between all three signals.
We can also consider the full bandwidth over the range of frequencies in which marsquake events have been observed by SEIS (0.1 Hz to the antialiasing filter cut-off of 40 Hz). A similar correlation between the seismic and atmospheric variable is seen, as shown in Figure 4 for both VBB and SP data for sols 361-362.
In panel (c), the acceleration response envelopes (frequencies of 0.1-1 Hz, 1-8 Hz, and 8-40 Hz) are here shown for all three components; the east, north and vertical, plotted against the wind speed. These sols take place in the early summer (solar longitude of L s = 115°), in contrast to mid spring depicted in Figure 2.
Here, the typical three-regime cycle of the weather is repeated, however, the ground acceleration indicates changes in its response. This is further addressed under the seasonal variation in Section 5.1.2. Each weather regime is well described by a power-law relationship over increasing bandwidths, similarly to as observed above in Figure 3. There is a pronounced increase in each slope with increasing frequency, indicating a higher sensitivity under the same weather fluctuations. Nevertheless, all bandwidths still share a windspeed threshold, bounding the acceleration envelope to the ambient noise limit of that frequency range. The threshold here for these two sols 361-362 has increased to ∼2.8 m s −1 , slightly above the threshold observed for sol 237-239 at ∼2.4 m s −1 .
To track how frequency relationships change with wind speed, panel (d) of Figure 4 shows the average power spectrum of the three-component seismic signal as a function of the wind speed from 0.1 to 10 Hz on the VBB and from 10 to 40 Hz on the SP. The dominant feature in this representation is the concentration of energy in mechanical modes evident as horizontal features at their resonant frequencies. These prominent modes are artificial and associated with the lander body, solar panels, instrument deployment arm, load shunt assembly, and tether, coexciting with surges in atmospheric activity. Such features are also observed on Earth, for example, resonant modes in ocean-bottom seismometers attributed to head-buoy cables strumming from fluid motion, even at moderate current velocities (Stähler et al., 2018). Here, a further characteristic of the wind speed threshold is apparent: the abrupt attenuation of modes below 10 Hz at the 2.8 m s −1 wind-speed threshold. The high-frequency bandwidth (8-40 Hz) exhibits increased scatter in the acceleration compared to the lower frequencies in all three regimes of the weather's diurnal pattern, as illustrated by Figure 4c. This is due to the persistence of multiple superimposed resonant modes and features below the wind-speed threshold at higher frequencies, suggesting an excitation mechanism other than wind or pressure. The variance below the threshold is even further amplified by the frequent occurrence of a distinct class of spikes, called "donks," impulse-like high-frequency (>8 Hz) occurrences likely from elastothermal relaxations that excite the mechanical resonances as seen in Figure 4d Lognonné et al., 2020). These very short-duration excursions cannot be attributed by the energy observed in either atmospheric variable.

Comparison of Vertical and Horizontal Ground Acceleration Response
The horizontal east and north components follow similar power laws to the vertical for all the frequency bandwidths examined in Figure 4c. The transition between the regimes occurs at the same wind speed for all three components, with the threshold at ∼2.8 m s −1 distinguishing the quiet period for sols 361-362. A further corner point at 4 m s −1 seen in Figure 4d represents the nighttime regime in our exemplar sols up until sunrise. During this period, the 2.4 Hz resonance is excited below a boundary set at 4 m s −1 particularly CHARALAMBOUS ET AL.
10.1029/2020JE006538 13 of 35 on the vertical component. This excitation disappears for higher wind speeds along with the 1 Hz tick noise-electrical crosstalk between the seismic channels and the temperature house-keeping channels )-observed on the N/S and Z. Above 10 Hz, it can be seen that there are persistent modes even below the wind thresholds, with their intensity increasing with wind speed for all three components. For example, there is a particularly strong resonance at 25 Hz.
However, at the lowest frequencies (<0.3 Hz), the seismic power increases with wind speed only on the horizontal components, likely through a tilt injection as described in Murdoch et al. (2017). Even during the more turbulent parts of the day, the higher-and lower-frequency injections on the horizontal seismic components maintain a relatively undisturbed frequency gap, only beginning to merge at ∼0.3 Hz during the stronger wind speeds. Owing to the injection by separate mechanistic processes, this frequency noise isolation suggests a possible bandwidth for event detection in the horizontal components even during more turbulent periods, something that is not applicable to the vertical. No events have been recorded to date within this frequency band during the daytime.

Seasonal Changes
The diurnal comodulation relationships are driven by the cycle of convective daytime and shear-driven nighttime turbulence. The Martian seasons trigger a change in this cycle. Figure 5 compares the wind speed and envelope correlations for the (previously examined) bandwidths 0.1-1 Hz, 1-8 Hz, and 8-40 Hz calculated with 50 s windows on data from sols 98-99, 290-291, and 361-362 corresponding to the Martian late winter, late spring and early summer sol, respectively (in the northern hemisphere). The selection of these sols is based upon availability of 100 sps seismic data. It can be immediately observed that the correlation relationships between wind speed and acceleration response envelopes do not follow a single trajectory, but rather exhibit seasonality with shifting patterns in the expected three-regime diurnal cycle of the weather.
The three cases of Figure 5 offer a good sample of the different environmental conditions in the first 400 sols of operation ( Figure 7 in Spiga et al., 2020). During the late-winter sols 98-99, close to spring equinox, InSight experiences high surface temperatures and particularly low ambient wind speed. Conversely, during the late-spring sols 290-291, close to summer solstice, surface temperatures reach a seasonal low while ambient wind speed is strong. The early summer period for sols 361-362 features rising surface temperature again, with ambient wind speed remaining strong. Of all those three cases, sols 98-99 are characterized by the strongest wind gustiness (Figure 9 in Spiga et al., 2020).
Sols 98-99 are turbulent throughout, the mean ambient wind speed is low, and the relationship is maintained predominantly within a single regime described by a single power law. Because the surface temperature is at its highest during these two sols (Mueller et al., 2020) and wind speed is low , this could lead to purely buoyancy-driven convection throughout the sol. While sols 98-99 fall mainly within a single regime, there is considerable variance. Furthermore, the noise floor is not reached due to the absence of a calm period. The high variance across this period is seen by different sensitivities to the ground acceleration for the same wind speeds and neighboring local times of occurrence. This increased sensitivity is observed particularly in the nighttime up to dawn.
For sols 290-291, the acceleration response trajectory during nighttime is split into at least three paths. During the distinct inertial oscillations caused by the low-level jet-induced weak turbulence, two of the pathways are observed: the first emerges around midnight and a further one at 4 a.m. The latter transitions to a flat region before re-emerging shortly after sunrise from ∼6 to 9 a.m. LMST, characterized by increased seismic injection compared to other sols, appearing as an outlier with a clear trajectory. Although it is well matched to the steeper power law as observed in other sols, this divergence reveals a prominent feature. Upon closer inspection, this behavior matches a reversal of the wind direction during the weakly turbulent, near laminar flow of the low-level jet stream. This regime matches to a power-law exponent of 2 for 0.1-1 Hz, and further continues into the daytime turbulence. The quiet period is entered at the lowest wind-speed threshold close to 2 m s −1 , in comparison to the other two sols shown in Figure 5 which occur at higher wind-speed thresholds. The flat slope of the quiet period is more strongly evident at the two higher bandwidths, with the "donks" apparent in a cluster from 8 to 40 Hz, before the re-emergence of comodulation after midnight for a period of time, returning to the flat region for a short period before the injection rises again. For the 0.1-1 Hz band, the whole sequence almost lies on a single power law aside from a branch of higher sensitivity during the morning.
On sols 361-362, the three regimes are present similar to sols 237-239 with clearly defined trajectories. The nighttime regime here can be seen gradually climbing up along a steep power-law trajectory. After sunrise, the morning turbulence injects the highest energy but follows a shallower power law with a ≈ 2 for the 0.1-1 Hz band, with this slowly shifting in its intercept until sunset. At the transition to the quiet period, this returns to a steeper power law with lower injection after noon, likely due to an atmospheric density shift, giving the greatest separation of the daytime branches. It then promptly enters the regime with the wind speed below the threshold, identified here at 2.8 m s −1 . The lack of injection during the quiet period is distinct at all frequencies with the highest wind-speed threshold of the three periods.
For the latter sols 290-291 and 361-362, there is an increase in both ambient wind speed and turbulence which is observed in the seismic signal. The regimes are more clearly separated with more variable power-law relationships. The daytime regime preserves the power-law relationship with an exponent close to 2. The power law begins during the strongest ground deformations early in the day, while slowly shifting down to lower responses under similar wind magnitudes in the afternoon. This negative shift in acceleration response could be related to the shift in the air density through the day, triggered by steep temperature gradients at the onset of sunrise and sunset. Denser air in the nighttime to early morning would exert a higher dynamic pressure on the surrounding elements inducing a stronger lander-coupled motion for the same wind speeds experienced through the day. This can be seen as the high variability in the acceleration envelope during the nighttime under relatively low wind speeds, forming a steep power law relationship. In contrast, less dense air caused by daytime heating will weaken wind-induced lander motion. Since the evolution through each regime essentially maintains a single power-law exponent only with a shift in the intercept, this could provide an alternative way of determining the air density variation at InSight.
In summary, these changes in the diurnal variation previously discussed in Section 5.1.1 are driven by the changing wind profile across the seasons, shown in Figure 5. For sols 98-99, the wind speed is generally low with indistinct regimes and a constant, low level of turbulence as measured by its high wind gustiness . For late spring into early summer, the turbulent daytime increases in length and strength.
As a result, it can be concluded that the seismic signal's sensitivity to the environment depends on specific weather conditions such as atmospheric pressure, wind direction, turbulence intensity/gustiness, air and ground temperature. It should be added that these conditions also affect the TWINS sensor calibration. TWINS is calibrated to retrieve the Reynolds number at the sensor position and shows a lower measuring threshold for Reynolds numbers of approximately 100. This translates to a seasonally varying wind-speed retrieval threshold value at sunset given the pressure and temperature evolution. The measured wind-speed threshold would vary at sunset from 2 to 2.8 m s −1 for sols included in this study. Complete seasonal coverage by InSight over 2 Martian years will help confirm the trends observed in this study. Figure 6 shows the correlation coefficient between the environmental and seismic signals, over the full bandwidth for wind and as a function of frequency for the pressure, calculated as a correlation matrix. The broadband noise injection induced by the atmosphere is correlated with the wind (top row) and within and across frequencies for the pressure and acceleration envelopes. The correlation coefficient matrix can be regarded as an estimate of the cross-frequency power comodulation. This matrix is constructed with the correlation coefficients of the spectral power variation in each overlapping frequency band of logarithmic spacing, between 0.01 and 10 Hz for the VBB Z and 0.01-4 Hz for the pressure. The wind is already in the appropriate domain for comparison, and the correlation to wind against the other two variables is shown at the top row. This establishes that cross-frequency power variations are driven by the atmosphere. Whereas comodulation breaks down below 0.1 Hz in the ground acceleration due to multiple glitches in these long period windows, the pressure to wind relationship is maintained.

The Relationship Between Pressure and Wind
These atmosphere-seismic relationships were illustrated in detail in Section 5.1.1, where Figure 3 showed that the pressure envelope exhibits a similar correlation to the seismic acceleration envelope as wind speed. In Figure 3b, the wind speed was plotted against the pressure envelope for multiple bands of increasing frequency, 0.01-0.1 Hz, 0.1-0.5 Hz, 0.5-1 Hz, and 1-4 Hz. It can be seen that for >0.5 Hz, above a 2.4 m s −1 threshold for sols 237-239, these relationships correlate with a single linear relationship a ≈ 1, while the exponent a becomes increasingly larger at lower frequencies. The increased spread observed at lower frequencies is also reflected by the slightly lower correlation coefficient reported by the top row of Figure 6. We investigated this relationship across multiple seasons and found a similar correlation.
As a result, we can see that the pressure envelope provides similar information to the wind speed above the seasonally variable wind speed threshold. The pressure sensor is designed to observe static pressure up to a sampling rate of 20 sps and so it is shielded by an inlet designed to transmit only the static component . In Banfield et al. (2020), it was suggested that the pressure sensor may be sensitive to wind at higher frequencies due to reduced efficiency of the inlet shielding mixing a proportion of the dynamic pressure into the measurement which would then be indistinguishable from the high-frequency static pressure. Moreover, static pressure and wind speed will be correlated through turbulence inducing static-pressure variations . A precise uncoupling of these factors is beyond the scope of this work.
Through both instrumental constraints and turbulence, joint modulation of the pressure and wind data can be expected to reflect the dynamic pressure generating forcing on the lander or the WTS which couples through the regolith to the SEIS assembly. In Section 4.1, we outlined how the dynamic pressure drives the observed seismic response through drag forcing. As discussed, these equations are in turn dependent on different factors (e.g., wind Reynolds number, atmospheric-density changes due to seasonal and diurnal CHARALAMBOUS ET AL.

10.1029/2020JE006538
17 of 35 The distance between each pair of points represents the frequency band at which the pressure or VBB Z signals were band passed and envelopes extracted. The calculation for the correlation coefficient was subsequently applied to this extracted envelope. (Bottom) Envelope correlation matrix for sols 237-239 between the vertical ground acceleration and pressure for frequencies from 0.01 Hz up to Nyquist for each sensor's data rate acquisition, spaced in frequency intervals of half an octave. The mean correlation coefficient across the VBB Z frequency is shown above the matrix, while for the pressure frequency the mean values are shown to its right. VBB, Very Broad Band.
pressure and temperature changes, surface area via a wind-direction dependence and regolith response) which would give rise to an array of sensitivities for different atmospheric conditions, including mechanisms dependent on the detailed geometry of the lander such as vortex shedding.
We now discuss the wind speed threshold below which SEIS acceleration no longer clearly correlates with the wind speed and often reaches the noise floor in several bandwidths. This is the quiet period regime during which most seismic events have been identified. In Lognonné et al. (2020), it was stated that this wind speed corresponds to a drop in Reynolds number, decreasing the resultant injection. In addition, wind speed retrieval also becomes more difficult and so exact values should be interpreted with caution . Although the wind speed can vary below this threshold at this period, its fluctuations cannot be reliably quantified. Figure 7 visualizes the threshold for sols 361-362 by comparing the VBB vertical acceleration component to the pressure-frequency histogram using pressure as a proxy for wind over an increased bandwidth due to their observed comodulation. The threshold is indicated by the dividing gray dashed line at 2.8 m s −1 for these two early summer sols, a cutoff point that differentiates the response in the behavior for both pressure and VBB at all frequencies. It is noticeable that the lander resonances also diminish in power at this point, falling below the instrument noise floor with an abruptness that is not seen in the ambient broadband signal. This could be associated with the collapse of the daytime turbulence which provides excitation of the lander modes. This figure further illustrates the joint broadband nature of the environmental injection for energy across the 0.1-10 Hz bandwidth with the energy across all bandwidths correlated between the pressure and seismic signals.
The overlain transparent gray lines on the 2-D histograms of Figures 7b and 7c are the logarithms of the total power showing the comodulation of the seismic signal power and wind strength during the 2-sol observation time. The seismic signal reflects the rising acceleration power with frequency and the excitation of mechanical modes observed at >1 Hz which concentrate excess power, as indicated by the mean PSD in Figures 7b and 7c. For pressure, the signal rises at low frequencies, with no modes present. Notice the co-occurrence of the wind speed threshold, flattening both integrals at 2.8 m s −1 . Above this threshold, the increase in the summed power behaves differently for the two signals, rising linearly for the pressure in line with the relationships in Figure 3. In contrast, the increase above this threshold for the vertical acceleration is a steep power law with an exponent close to a = 4 (see Figure 5) up to ∼7 m s −1 . Both signals capture a breakpoint at this wind speed of ∼7 m s −1 , a value which here signifies the transition from the nighttime low-level jet to the turbulent daytime after sunrise, to the transition after sunset following the end of turbulence. Above this breakpoint, the seismic acceleration transitions into a power law with an exponent closer to a ≈ 2.

Predicting the Wind Speed and Direction
Earlier in Figure 6, it was identified that the ground acceleration and pressure could act as a proxy for wind speed due to the cross-frequency coupling. We exploit this again here in Figure 8. The envelope is extracted in 1 Hz narrow frequency bands, starting from 0.1 Hz up to 8 Hz for the VBB and 4 Hz for pressure to avoid antialiasing effects. The square root of the envelope of each incremental band is moment matched to the global moments of the wind over the 3-sol period for sols 237-239. Due to the consistent broadband characteristics of the noise induced by the atmosphere, both pressure and acceleration indicate that any frequency band is directly coupled to any other frequency band between pressure and VBB, at all times. This can be also realized by the good overlap of these scaled relationships, providing predictions with minor differences. Therefore, the energy content of both pressure and ground acceleration exhibits the same behavior and are good predictors of the wind. This is shown by the close linearity in the performance of the predictions against the measured wind speeds in the middle plots of Figures 8a and 8b.
On the right-hand side of Figures 8a and 8b, the evolution of the mean and variance in incremental frequency bands of 0.2 Hz is shown for both vertical ground acceleration and pressure. The mean and variance of the pressure drops with increasing frequency, while it increases for VBB. This is in line with the expected noise curves of the two systems Lognonné et al., 2019). The triangles in the rightmost panel of Figure 8a show synchronized covariations between the mean and variance which occur at the resonant modes about 4, 6, and 7 Hz, the tick noise at 1 Hz and a small perturbation of the 2.4 Hz ambient mode of the VBB Z seismic acceleration. The reader is also referred to Figure 1 for the visualization of the diurnal mode variations in frequency, as shown by the spectrograms of the VBB components. The diurnal temperature variation causes a shift in the frequency of a subset of the mechanical modes. Therefore, the variance will be higher within the frequency range encompassing these lander modes since they vary diurnally. On the contrary, the persistence of a constant amplitude tick noise in the quiet and intermediate regimes identified earlier in Section 5.1.1, promotes an inverse relationship with less variance and a higher mean. These modes therefore provide a good representation of the lander's diurnally evolving state and can also be used as a proxy of the excitation induced by the atmosphere.
We have also successfully established a prediction of the measured wind direction by TWINS from the seismic energies of the horizontal components. The measured wind direction by TWINS is determined by a combination of the measurements of the sensors and computational fluid dynamics of the lander-induced CHARALAMBOUS ET AL.

10.1029/2020JE006538
19 of 35 disturbances to the wind flow . Our approach exploits the acceleration induced by ground deformations from the lander to estimate the direction of forcing on the lander from the wind. We have estimated this by moment matching with a single mean and variance from the inverse tangent function of the horizontal acceleration envelopes of the VBB components E and N, band passed between 0.5 and 1 Hz and shown in Figure 9. The agreement between the VBB prediction and the actual wind direction partly breaks down during the steep temperature gradient just before sunset (see Figure 2) and through the late afternoon's quiet period. This corresponds to the lack of environmental injection which induces mechanical vibrations and falls steeply below the wind-speed threshold (Figure 7). This period also occurs at the windspeed threshold observed by both the pressure and VBB sensors. Our method therefore provides an independent confirmation of the atmospheric diurnal cycle, synoptic (day-to-day) variations and the seasonal cycle of wind direction and wind speed retrieved by TWINS. These predictions offer an alternative means of redundancy in maintaining observations of the atmospheric variables, both for short and longer timescales.

Discussion on Observed Environmental Sensitivity
In this section, we have examined the comodulation of the seismic ground motion with respect to the pressure and wind data to determine the level of environmental injection. The seismic acceleration envelope correlates with the pressure envelope and wind speed and follows a set of power-law regimes which vary seasonally. Furthermore, we have analyzed the correlation between the pressure and wind data in order to understand the core factors of the observed comodulation. Here, we will summarize the underlying mechanisms for the observed relationships.
CHARALAMBOUS ET AL.
10.1029/2020JE006538 20 of 35 As already stated, the main injection pathway is through the dynamic pressure, proportional to U 2 , which induces a forcing described by the drag equation . These forces are applied predominantly to the lander but also to other site components such as the WTS and the tether. This forcing then produces ground deformation which couples to the SEIS assembly . Such forcing can also be directly applied to the regolith producing a tilt and requires consideration of the wind direction and air density. On top of this direct wind-induced forcing, other coupling mechanisms such as injection from a pressure field advected by the wind (Sorrells, 1971) are possible. At different times, different mechanisms may dominate.
Over the course of a sol, three regimes are generally observed in terms of the relationship between the seismic and environmental signals: (i) a shallow power law with exponent near 2 during the daytime convective turbulence, (ii) a steep power law with exponent near 4 during the nocturnal low-level jet from nighttime through morning, and (iii) a quiet period where the wind speed drops below a threshold (in the Reynolds number) where no seismic response is observed. These regimes have been observed to vary across seasons and are weather dependent. For example, we saw that only the steeper power law is observed on sols 98-99. For the other sols, the daytime regime is maintained in the shallower power law, with an intercept value that shifts negatively, in line with changes in atmospheric density. At least one example, sols 361-362, indicates a correlation to the wind direction at sunrise, where a prominent divergent branch in the relationship merges the transition from nighttime to daytime regimes.
The comodulation approach has enabled this quantification without a full development of a theoretical model for each injection pathway. As stated, these pathways are numerous and have complex dependencies. For example, the drag coefficient depends on the Reynolds number which, in turn, depends on wind speed, air viscosity, and density. This dependence can account for several diurnal and seasonal patterns. The slow decrease in sensitivity of the daytime turbulence regime can be attributed to such a factor. On a seasonal scale, the air density on Mars can change by 30% over a year, with a diurnal density varying by 50%, between daytime and nighttime. This not only causes a change in sensitivity but may also cause the observed variation in the cutoff threshold which changes between 2 and 2.8 m s −1 . Different atmospheric properties and weather patterns do induce some of the variation observed. On top of the variation in the initial forcing, they may also influence the coupling mechanisms. The InSight lander and deployed SEIS are supported CHARALAMBOUS ET AL.
10.1029/2020JE006538 21 of 35 by a poorly sorted, thin layer of unconsolidated sand that is underlain by a cemented fine-grain matrix, called duricrust, of several centimeters thickness . Beneath this layer, the regolith consists of buried rocks from impact ejecta with a distribution that increases in size with depth (Charalambous, 2014). The near-surface stratigraphy of the regolith and the lander itself are subject to steep temperature gradients. As a result, their thermoelastic properties are affected, as seen in the diurnal variation of the lander resonances. On top of this, radiative temperature effects on the lander may also impact TWINS winds retrieval. Further analysis complemented by surface brightness temperature sensing (Mueller et al., 2020) and air temperatures may provide a deeper insight to the properties of the Martian regolith deformation and lander dynamics. We consider that the comodulation approach has short-circuited the requirement for this full model but posit that these theoretical aspects are encompassed in the observed variation. The precise quantification of this physics is the subject of further work.
Our comodulation method can be further exploited to provide the resolution which TWINS lacks during the passage of convective vortices. Charalambous et al. (2020) identified wind retrieval issues due to perturbations and die saturation from high-vorticity conditions during convective vortex encounters. Our method can provide continuous wind predictions and therefore provide estimates of wind-peak speeds to fill in the missing wind data during such turbulent conditions. This may allow estimations of high-frequency wind speeds that have never been resolved before, to better identify the sudden peaks in the wind speed and infer the fluid threshold for the onset of particle detachment from the surface, for both sand motion and dust lifting . This will provide an insight into the long-standing paradox of aeolian transportation and the mechanistic relationships in the coupling between atmospheric processes and the surface of Mars.

Environmental Injection Analysis for Seismic Events
The InSight mission includes the Marsquake Service (MQS) which identifies events as part of mission operations, nearly 400 in the current catalog V2 Giardini et al., 2020;InSight Marsquake Service, 2020). The understanding of the injection of environmental sources into seismic data is key for the identification and analysis of seismic events. This is especially true in light of the Viking mission (Anderson et al., 1977;Lorenz et al., 2017), and the key Level 1 science goals of this mission are dependent on the identification of seismic events in the downlinked data (Banerdt et al., 2013;Lognonné et al., 2019).
In this section, we apply our comodulation analysis to the task of quantifying the energy injection from wind and pressure into events; that is, assessing the independence of the seismic events from the environmental noise. To do so, we first show how the moment-matching estimation is adapted to only take conditions local to the event into account and thus provide a suitable regression between seismic and environmental envelopes. We subsequently propose two SNR metrics to apply to these matched envelopes in the energy domain. The initial SNR and quality metrics so far implemented by the MQS are based on the amplitude of the seismic signal and its separation from the signal background. On the other hand, the proposed comodulation-based SNR incorporates the independence of the event from environmental sources providing a new dimension for analysis. We next apply these methods to several marsquakes, showing how the proposed approach can be tailored to an event's characteristics. We also apply these methods to the entire event catalog for comparison.

Moving Moment Matching for Envelope Prediction
In order to implement our comodulation approach to provide an SNR metric, we must first develop a regression between the envelopes obtained in Section 4.2. The method of moments has already been introduced for this task in Section 4.2.1, and it has been used to compare the wind speed, seismic, and pressure envelopes within distinct power-law regimes in Section 5.1.1. The method of moments is an estimation method which matches the population moments (up to second-order in our case) of two time series. This is a regression, where the moments of a time series are matched to the reference. To analyze seismic events, we moment match to the conditions local to the event. This is because the variance and mean of wind speed, pressure, and seismic acceleration are not constant throughout each regime, that is, each envelope time series exhibits unconditional heteroskedasticity, with a time-varying variance (Bollerslev, 1986). For forecasting equations, for example, in equity markets (Engle & Sokalska, 2012), heteroskedasticity imposes a serious challenge due to unpredictable variance measures (Bollerslev et al., 1992).
In our case, heteroskedasticity can be expected mainly from the identified changes in atmospheric conditions in Section 5.4, for example, air density variation throughout the convective turbulence daytime regime or wind direction reversals. To this end, a finer-grained approach is required to assess environmental injection into potential marsquakes which also vary over a smaller scale than these regimes and have typical durations ranging from 5 to 30 min. Therefore, the moment matching method is extended to operate over a sliding window, enabling an adaptive regression where the estimation is locally linear, akin to a manifold. This proposed moving-moment-matching procedure captures the variation in the mean and variance over time in order to create a smooth prediction which adequately represents the current conditions as they evolve through the Martian day.
To calculate the prediction, consider the moving average filter defined as the operation: which outputs the mean value of the input e[t] over K time steps before and L following time steps. Next we define the moving-variance operation: which similarly outputs the variance of the input e[t] over K time steps before and L following time steps.
Using these operations, the moment matching in Equation 2 is then done on a sliding basis, where the matched time series

An SNR Metric for Seismic Events
The regression between the seismic and environmental signal envelopes is used to estimate the independence of the observed seismic signal from the weather. This allows us to obtain an SNR for the events in the seismic catalog.
The first proposed SNR metric to compare two events is the ratio of the instantaneous energy, given as: where e[t] and The peak value of SNR2[t] during the event then yields a second metric that averages over the event and thus avoids instantaneous features to better assess the overall power injection. However, an averaging may CHARALAMBOUS ET AL.

10.1029/2020JE006538
24 of 35 Figure 10. A diagram illustrating the steps of the environmental independence SNR algorithm. The pressure and seismic acceleration time series are first transformed to a spectrogram and the envelopes obtained as in Section 4.2. The mean Ex and standard deviation σ are extracted from the seismic acceleration envelope, pressure envelope, and wind speed in order to perform a regression through moment matching. The ratio of the resulting matched envelopes then yields an SNR. SNR, signal-to-noise-ratio.
bias the result if there is an extremely large divergence. For a given event, this is reported for both regressions to wind and pressure and the full process is summarized by Algorithm 2.

An Application of the SNR Analysis to Recorded Marsquakes
In this section, we apply our proposed SNR metrics to a set of marsquakes identified so far. In doing so, we show how our comodulation analysis identifies the divergence in the seismic acceleration envelope from the wind speed and pressure envelope, which is then quantified by the proposed SNR metrics. This is extremely useful in light of the fact that envelope alignment and analysis has so far been particularly fruitful for the interpretation of what InSight has observed in terms of the structure of Mars Lognonné et al., 2020).
The events identified by the MQS have been categorized into the following types: 1. Low frequency (LF)-energy below 2.4 Hz 2. Broadband (BB)-energy mostly below 2.4 Hz with some slight excitation above 3. 2.4 Hz, and events with energy solely around the ambient resonance which is proposed to be of a potentially geophysical origin  4. High frequency (HF)-energy centered around and above 2.4 Hz, considered to be larger versions of the 2.4 Hz events 5. Very high frequency (VF)-energy increasing with frequency in the horizontal components often exciting the 2.4 Hz resonance For more details of these events, see Clinton et al. (2020) and Giardini et al. (2020). Figure 11 shows the spectrogram, matched envelopes and SNR metrics for the events S0173a, S0235b, S0185a, and S0128a, obtained from Algorithm 2 implemented with the following parameters: 1. Event bandwidth f min −f max in Table 1 2. Window length W len = 50 s and overlap OL win = 90% between spectrogram slices 3. Spectrogram slice PSD number of averages N av = 2 4. Length of moving moment matching K MM = 1,000 s and L MM = 0 for longer period LF/BB events and K MM = and L MM = 0 500s for the HF/VF/2.4 Hz events 5. Number of standard deviations to exclude from moment calculation σ = 5 6. Length of SNR averaging K SNR = 500 s and L SNR = 500 s These events are a good subset to demonstrate the applicability of comodulation to event analysis. The aim is to highlight how to apply the proposed technique to understand the event energy divergence and the various pitfalls. The S0173a and S0235b events are the two category A (highest quality) events observed so far Giardini et al., 2020). From the examination of the matched envelopes, both show a clear divergence of the seismic energy from the wind and pressure. The power of the moment matching for event analysis is well demonstrated by S0173a as it occurs during a noisy period and so the envelopes track well until the event where the SNR metrics are able to be obtained. The transient spikes in the SNR1 values are due to VBB sensor glitches (Scholz et al., 2020). On the other hand, the S0235b event highlights the data dependency issues of this metric. At the beginning of the time window demonstrated at approximately 11:20, a broadband high amplitude signal indicates a mechanical recentering process of one of the VBB components, with a clear divergence from the atmospheric noise. About 10 min prior to the event onset, the wind speed drops very low and this divergence is captured by both the wind and pressure SNRs. In fact, the wind SNR2 is seen to become biased for the surrounding samples, leading to an overestimation. Here, the wind speed has actually dropped below the wind speed threshold of 2.4 m s −1 discussed in Section 5.1.
For a high quality determination, the parameters K MM , L MM , K SNR , and L SNR must be adjusted. To show this adjustment consider event S0128a in the fourth column, a VF event. For such higher-frequency events, we CHARALAMBOUS ET AL.
10.1029/2020JE006538 26 of 35 Figure 11. The derivation of the SNR metrics for a selection of quality A LF (S0173a) and BB (S0235b) and quality B LF (S0185a) and VF (S0128a) events InSight Marsquake Service, 2020). The top row panels show the vertical (Z) spectrogram of each event (except S0128a, showing north [N]), the second row the seismic envelope (black, event period taken from InSight Marsquake Service (2020) is highlighted with the thick black line) and movingmoment-matched wind (cyan) and pressure (red), third row the SNR1 metric of Equation 7 for wind and pressure and the fourth row the SNR2 metric of Equation 8. The black stars refer to the peak picked for the SNR metrics. Note that in many of these events the wind speed drops below the threshold discussed in Section 5.1. SNR, signal-to-noise-ratio.
choose a shorter moment-matching window of K MM = 500 s in order to capture the faster variation in the bandwidth, owing to its higher sensitivity observed in Section 5.1. Note that this envelope is calculated on the north component, as for the higher frequency events the horizontal power is greater than vertical . The selection of the window length is critical in our SNR estimation; too short a duration and the algorithm risks overfitting, whereas a long enough window risks contaminating the statistics of different regimes observed throughout the Martian sol introducing back heteroskedasticity. Therefore, these selected window lengths provide the best SNR calculations without any risk of overfitting or regime contamination. Now consider the event S0185a in the third column of Figure 11. This is a relatively low amplitude signal, however, with our analysis its separation is clear and it receives a relatively high SNR. This indicates that there is little corruption of the data by environmental injection and so any estimated parameters from this event would be relatively reliable with respect to its strength. This shows the power of our approach to indicate a level of confidence in derived interpretations.
The proposed SNR metrics are based on the measured data and therefore they rely on good data quality in all data. For example, consider the above discussion for S0235b. These data issues include: Glitches in either the seismic, wind speed, or pressure data. These corrupt both moments, the mean and variance, in the statistics within an event and inject power to the signal, leading to an incorrect SNR determination. These glitches are prevalent in S0173a in Figure 11. Drop out/missing data, which renders an assessment impossible. Very low amplitudes (i.e., low wind speed), which causes an instability in obtaining a result (e.g., S0235b in Figure 11). This often occurs when the wind speed drops below the threshold and so interpretation of these values should be cautious. In addition, wind retrieval issues often result in no data points during which the algorithm depends on interpolated values. Figure 12 shows the moving moment-matched envelopes for a selection of key LF/BB and HF/VF events as identified in Clinton et al. (2020) and Giardini et al. (2020) with the SNRs given in Table 1. The event bandwidths f min and f max , manually optimized for each event, are also provided in Table 1.
CHARALAMBOUS ET AL. Abbreviations: MQS, Marsquake Service; SNR, signal-to-noise-ratio.   Shown in Figure 12 The key value of this measure is that it indicates the independence from atmospheric injection as opposed to simply the strength of the seismic signal. On top of the SNR value, the moving moment-matching regression yields the ability to assess correlation of the envelopes. For the largest events, while they are already clear, this analysis can show the degree to which portion of an event is independent. For smaller events, this metric may more appropriately represent their significance. For example, event S0189a is a low-energy event with a higher SNR from our analysis than from the amplitude-based MQS SNR. The divergence can be best seen as the envelope does not correlate with the pressure envelope or wind speed. On the other hand, CHARALAMBOUS ET AL.
10.1029/2020JE006538 28 of 35 the interpretation of an event which exhibits some correlation with the pressure envelope or wind speed must take into account a potential injection, for example S0167a in Figure 12.

An Operational SNR for the Event Catalog
In this section, we discuss the implementation of our proposed SNR metric for the entire MQS catalog. In this case, a consistent set of parameters must be used as the SNR must allow equivalent comparisons between events. This is also the case for the SNR already implemented in the MQS catalog, calculated based on the ratio of the power (obtained from the PSD) in the event bandwidth during and just prior to the event . To this end, the Algorithm 2 was applied to the entire catalog with the parameters used in the previous section and the frequency bands: (i) f min = 0.2 Hz and f max = 0.5 Hz for LF and BB events and (ii) f min = 2.2 Hz and f max = 2.6 Hz for the HF, VF and 2.4 Hz events. These are consistent with the MQS approach to aid comparability. For the low frequency group of events (LF and BB), there is a strong relationship between the SNR2 and MQS SNR for both atmospheric variables with power-law exponents just under 1. This shows that the low amplitude events are actually generally independent in this bandwidth and, while still significant, the SNR2 for the larger events is slightly lower owing to the averaging over the event of our method compared to the PSD-based SNR. At lower MQS SNRs, the matched-moment values are relatively higher, with the divergence from environmental injection more apparent than in the seismic signal alone, thus providing a higher degree of confidence.
In comparison, the higher frequency group shows a different relationship. The power-law exponent is around 2 for the 2.4 Hz, VF and HF groups against both atmospheric variables, showing that the larger events tend to have a higher SNR2 than that derived from the PSD and the lower MQS SNR events have a lower SNR2 and are more difficult to distinguish from the environment than low frequency events. This suggests that this bandwidth is more susceptible to environmental injection. This is in line with Section 5.1 which showed an increased environmental sensitivity with frequency.
The SNR2 derived from the pressure is often more robust than from wind. For example, in Figure 13, the pressure-based plots show a tighter clustering than the wind based. This is due to the limitations of TWINS inherent to wind sensing. At wind speeds below the threshold of 2.8 m s −1 , the TWINS measurement can be inaccurate . Moreover, the higher sample rate of the pressure sensor data is useful. As discussed in Section 5.2, we use the pressure envelope here as a partial proxy for the dynamic pressure.
High-frequency events have larger horizontal than vertical energy, owing to their scattering nature . However, the 2.4 Hz resonance is vertically polarized and so we do not see this directly in the SNR2s obtained with the fixed parameters for the catalog. This feature can be seen in the bandwidth-specific analysis in the case of S0128a where we have chosen to display the SNR2 from the North component.
This implementation of the proposed SNR2 metric provides an automated system to assess events in the catalog InSight Marsquake Service, 2020) in terms of independence from environmental injection. As this is an automatic operational approach, some SNR2s may not be truly representative owing to the unpreventable data quality issues outlined above. Moreover, the implementation requires a set of fixed parameters as inputs for Algorithm 2. These parameters are outlined in Section 6.2.1 with the bandwidth f min = 0.2 Hz to f max = 0.5 Hz for the lower frequency group and f min = 2.2 Hz to f max = 2.6 Hz for the higher frequency events. However, in general, this implementation has been reasonably robust and therefore acts to provide a necessary reference for catalog users to identify suitable events for their analysis and to understand their distribution in terms of environmental independence. Furthermore, a methodology is in place where this analysis can be fine-tuned for a specific event. This is demonstrated in the sequel, Section 6.3. The SNRs for the entire catalog are reported in the Supporting Information (SI).

Discrimination Between Atmospheric and Other Nonseismic Source Injections
We have proposed a method for event analysis based on a moving-moment-matching regression. This has been shown for event analysis and implemented as an automatic algorithm to obtain an SNR. In this section, we employ this method more generally to identify episodes of loss in correlation in the energy transferred by wind or pressure to that of the observed seismic signal. Figure 14 illustrates the seismic envelope against the matched pressure energy and wind speed over four distinct sols and shows various prominent features observed so far at the InSight landing site. Here we demonstrate how identification of a divergence in the moment-matched atmospheric and seismic envelopes (calculated here as the magnitude of the three-dimensional ZNE vector of the acceleration response as   2 2 2 ( ) Z N E ) allows us to discriminate different potential sources of noise contamination to the combined seismic signal from all directions. Figure 14a illustrates this approach for the quality B, VF event on sol 128. The arrival of this particular VF event is seen with the seismic envelope diverging from the expected match, indicating that the signal source was not due to environmental injection. The seismic envelope during the event is significantly above the acceleration-equivalent environmental injection, suggesting a nonsite source for the seismic signal, uncontaminated by wind or pressure energy. The robotic-arm excursion can be seen following the event promptly, performing a maneuver from which the Instrument Deployment Camera (IDC) achieves the tau pose, in order to perform imaging for variations above the horizon in the atmospheric optical depth (Trebi-Ollennu et al., 2018).
A BB event from sol 133 can be seen diverging in Figure 14b shortly after a broadband glitch, with both of these features indicating a nonatmospheric origin. After the occurrence of some further VBB glitches later on, a "pressure burst" diverges from the seismic and wind match. Pressure bursts are a rare phenomenon of nighttime occurrences exclusively observed on the pressure measurements with a source mechanism that is not yet understood . Figure 14c demonstrates the discrimination between features of analogous spectral signatures from sol 253. A broadband impulse resulting from the robotic-arm motion performing a push with the scoop onto the regolith is observed at 19:50 UTC. The signature resembles the "glitch-like" spectral character of a weak convective vortex occurring at 20:10 observed by SEIS as a tilt in the ground response (Kenda et al., 2017;Murdoch et al., 2021). Notice how the broadband injection from the lander vibration induced by the robot-arm motion is not observed in the atmospheric energy, in contrast to the vortex that is matched by the envelope forecasting equation. Another resemblance between spectral features can be observed from the vibrations induced by a subsequent robotic-arm motion during which the arm retracts further away from Figure 14. Discriminating origins of atmospheric injection from other sources at the InSight landing site. (a) Sol 128-VF event and robotic-arm motion switching to an atmospheric tau-measuring position. This is the power in the frequency band 1-8 Hz of the VBB ENZ signal. (b) Sol 133-glitches, a BB event with dominant energy below 1 Hz and a "pressure burst," a poorly understood excess energy injection for which SEIS is known to be insensitive to Banfield et al. (2020). This is the power in the frequency band 0.1-1 Hz of the VBB ENZ signal. (c) Sol 253-robotic-arm motion with the scoop pushed onto the regolith. A pressure bump is seen later matching to the SEIS envelope, however, notice the similarity in the spectral character of the robotic arm motion 5 min before. A convective vortex at 20:10 appears "glitch-like"-notice the resemblance of the vortex's signature to the HF spectral character of the regolith push at 19:50. This is the power in the frequency band 1-8 Hz of the VBB ENZ signal. (d) Sol 334-a VF event amid a forest of donks. This is the power in the frequency band 5-40 Hz of the SP ENZ signal. the surface at 20:01. The duration and high-frequency spectral characteristics are not dissimilar to the pressure-energy injection 3 min later, marked as ΔPressure.
Donks, prominent impulse-like high-frequency features exciting lander resonances and described earlier as occurring intermittently at all time regimes of the day are seen on Figure 14d for sol 334 during the sunset, apparent at frequencies above 8 Hz from the SP ZNE acceleration response magnitude. Amid this forest of donks, the S0334a VF event occurs with its envelope diverging from the atmospheric match.

Discussion and Conclusion
The auxiliary sensors of InSight were part of the mission's design to help identify the effect of the atmosphere on the prime seismic payload of SEIS (Banerdt et al., 2013). However, the mapping between each sensor's data is not entirely straightforward for two major underlying physical reasons. First, the coupling between the atmosphere and the lander is a result of the sum of the atmospheric forcing on all the surfaces of the lander and the regolith which will vary according to the wind direction and the level of turbulence Myhill et al., 2018;Teanby et al., 2017); InSight's wind sensor's can only provide at most a two-point (at different locations and height than SEIS) measurement. Second, the transmission of the forces through the regolith to SEIS will depend on the multilayer regolith  properties that will in turn depend on the environment in the near-surface region. The quantification of the injection can therefore be seen as a hidden-variable problem, breaking the 1-1 mapping that might be hoped for between the environmental and seismic signals.
Two possible approaches to reveal these variables are through modeling and analog testing. The first is challenging due to the geometrical and time-resolution requirements to allow a convergence on the dynamics. Analog testing would also be a major undertaking, requiring a representative model of the lander, regolith and atmosphere under Mars conditions. Furthermore, we are looking at environmental-injection levels on Mars which are below the lowest observed background on Earth over much of our observational bandwidth: the deployment of SEIS was successful in physically attenuating environmental injection to levels better than can be observed in the quietest vault on Earth. Hence, it is impossible to provide analog testing that would elucidate injections at the levels we have observed with InSight, and assumptions of linearity would be required to extrapolate from Earth testing to the mission environment.
The approach outlined in this work is instead to use the data itself to reveal the effect of these hidden variables as a time-dependent coupling between the environmental and seismic signals. We use the observed comodulation of the environmental and seismic signals to estimate the power transferred from the environment to SEIS. The method of moments applied here allows this coupling to be quantified within related regions even in the absence of an injective function which would provide a phase relationship between the signals to allow the estimation of a transfer function.
In terms of the qualitative relationship revealed, the environmental injection is generally broad band, except at the resonances of the lander, with the environmental signal variance injecting across much of the bandwidth of the seismic signals. Only during the quietest periods is the sensitivity of seismic detection set by the noise floor of the VBB and SP sensors of SEIS. The relationship between the signals varies with both the diurnal and seasonal cycle, with varying power laws identifiable as a function of frequency. The diurnal covariance is cyclic, with three seasonally evolving identifiable regimes of which the quiet evening period is of prime interest for seismic event detection on the planet. This work has quantified the relationship by matching the first and second moments of the environmental and seismic signal variance in a time-advanced window longer than the observed duration of possible seismic events. Such a matching allows an estimated attribution of signal in these events between environmental and tectonic sources and the derivation of a corresponding environmental SNR. These SNRs show that while the lower frequency events generally have a high enough SNR consistent with a tectonic source, at higher frequencies above 1 Hz many events show significant environmental injection.
The limitations of this approach are two-fold. First, during the quietest periods, the environmental sensors themselves are at the threshold for reliable quantification. It is therefore quite possible that environmental injection could still be occurring, but can only be detected on SEIS signals themselves. An extension of the existing approach would be to use the broadband characteristics of the injection, and its enhancement at the known resonances of the lander to give an internal attribution; an excess energy at such a resonance comodulating with an environmentally undetectable signal would then allow an estimate of the partition of the signal energy between tectonic and environmental sources. The ability of the seismic signal to predict the wind speed and direction suggests the viability of such an approach.
While the physical attenuation of any environmental injection into a seismic signal is preferable, mission constraints will always limit the ability to provide such attenuation comparable to a terrestrial deployment. However, we show here that environmental sensors can provide a quantification of environmental injection, even if the pathway for such an injection is only partially understood. Sensors will usually levy a lighter resource burden on a mission than deployment mechanisms and in addition provide science return in their own right. Hence this work provides a motivation for including sufficient environmental sensors in future mission payloads to provide critical attenuation. The sensors should be chosen to measure the largest expected environmental injections, and have sufficient performance to quantify the injection down to the seismic sensors' noise floor. However, even without such dedicated sensors, we have shown it is possible to use different bandwidths of seismic sensors for seismic and environmental attribution, providing an internal approach to quantify such injection, and thus effectively allowing SEIS itself to build its own virtual vault.