Empirical H/V spectral ratios at the InSight landing site and implications for the martian subsurface structure

SUMMARY

The horizontal-to-vertical (H/V) spectral ratio inversion is a traditional technique for deriving the local subsurface structure on Earth. We calculated the H/V from the ambient vibrations at different wind levels at the InSight landing site, on Mars, and also computed the H/V from the S-wave coda of the martian seismic events (marsquakes). Different H/V curves were obtained for different wind periods and from the marsquakes. From the ambient vibrations, the recordings during low-wind periods are close to the instrument self-noise level. During high-wind periods, the seismic recordings are highly contaminated by the interaction of the lander with the wind and the martian ground. Therefore, these recordings are less favourable for traditional H/V analysis. Instead, the recordings of the S-wave coda of marsquakes were preferred to derive the characteristic H/V curve of this site between 0.4 and 10 Hz. The final H/V curve presents a characteristic trough at 2.4 Hz and a strong peak at 8 Hz. Using a full diffuse wavefield approach as the forward computation and the Neighbourhood Algorithm as the sampling technique, we invert for the 1-D shear wave velocity structure at the InSight landing site. Based on our inversion results, we propose a strong site effect at the InSight site to be due to the presence of a shallow high-velocity layer (SHVL) over low-velocity units. The SHVL is likely placed below a layer of coarse blocky ejecta and can be associated with Early Amazonian basaltic lava flows. The units below the SHVL have lower velocities, possibly related to a Late Hesperian or Early Amazonian epoch with a different magmatic regime and/or a greater impact rate and more extensive weathering. An extremely weak buried low velocity layer (bLVL) between these lava flows explains the data around the 2.4 Hz trough, whereas a more competent bLVL would not generate this latter feature. These subsurface models are in good agreement with results from hammering experiment and compliance measurements at the InSight landing site. Finally, this site effect is revealed only by seismic events data and explains the larger horizontal than vertical ground motion recorded for certain type of marsquakes.

Coda waves, Seismic noise, Site effects, Martian seismology, Marsquakes

1 INTRODUCTION

On 2018 November 26, NASA's InSight mission landed on Elysium Planitia, Mars, and deployed a set of geophysical instrumentation aimed at investigating the planet's inner structure and dynamics. Along with other instruments, SEIS (Seismic Experiment for Interior Structure, Lognonné et al. 2019) monitors the martian ground motion using one six-channels seismological station, which was deployed on the ground, around of 2 m away from the InSight lander.

The InSight landing site is located at 4.5024 N/135.6234E inside a degraded impact crater, the so-called Homestead hollow, in Elysium Planitia (Golombek et al. 2017, 2020c). The surficial geology of this site was studied before the mission (Warner et al. 2017), using mainly orbital imagery and analysis of rocky ejecta craters, and also assessed after landing (Golombek et al. 2020b; Warner et al. 2022). According to the pre-landing observations, a model proposed for the shallow subsurface structure at the InSight landing site consists of a shallow fine regolith layer, a second layer of coarse ejecta, a deeper layer of fractured basalt followed by a layer of more pristine basalt and, finally, a deep layer of possible weakly bonded sediments located at ∼200 m depth below the lander (Knapmeyer-Endrun et al. 2017; Pan et al. 2020). Knowing the shallow subsurface layers of a site and their elastic properties is relevant for understanding the recorded ground-motion amplitudes. In this regard, the amplification of the ground motion in a specific frequency range is affected by the waves propagating through soil and soft rock layers, commonly referred to as a site effect (Borcherdt 1970; Anderson et al. 1986). The analysis of analogue samples has allowed studying the mechanical properties of the martian regolith at the InSight landing site, suggesting that the top layer is made of soft material (Delage et al. 2017; Morgan et al. 2018). Therefore, ground-motion amplification due to a site effect might be expected at the InSight landing site.

On Earth, a classic technique used to evaluate ground-motion amplification is the horizontal-to-vertical spectral ratio (H/V or HVSR), which is defined as the ratio between the Fourier amplitude spectra of the horizontal and the vertical components of the seismic recordings. This technique, first introduced by Nogoshi (1971) and popularized by Nakamura (1989), allows the estimation of the fundamental resonance frequency and ground-motion amplification using seismic recordings of either ambient vibrations or earthquake motions. Even though this method has become popular in recent decades, the discussion on its strengths and limitations is still due to the absence of a clear theoretical basis (Bonnefoy-Claudet et al. 2008). In this context, Sánchez-Sesma et al. (2011) have proposed that the ambient vibrations can be represented as a diffuse wavefield containing all types of body (P and S) and surface waves (Love and Rayleigh). This diffuse field approach (DFA) has a strong theoretical background based on the principle of equipartition of energy (Weaver 1982; Margerin 2017). Under this assumption, the autocorrelation in the frequency domain is proportional to the imaginary part of the Green's function (GF), so the latter can be directly linked to the H/V curve. Several studies have applied this theory to explain the features observed on the H/V curve and have demonstrated it to be a robust approach (e.g. Lontsi et al. 2015; Piña-Flores et al. 2016; Spica et al. 2018; Bora et al. 2020).

The diffuse field regime can also be extended to the seismic event recordings, as shown by Margerin et al. (2009). Recently, Kawase et al. (2018) obtained similar H/V curves from ambient vibrations, the S-wave window and the late coda window, especially up to their first peak frequency. Thus, assuming the wavefield in the coda is equipartitioned, the DFA can be applied to model the H/V curve from this data window.

Since the landing, some attempts have been made to reveal the subsurface structure at the InSight landing site. First, the SEIS recordings of the hammering sessions of the HP³ (Heat Flow and Physical Properties Package) experiment were used to derive the seismic wave velocities in the first tens of centimetres below the ground surface (Brinkman et al. 2019; Lognonné et al. 2020). A re-evaluation of this data has determined the P-wave velocity to be 119 ± 29 m s⁻¹ and the S-wave velocity to be 61 ± 20 m s⁻¹ in the first ∼40 cm of the shallowest layer (Brinkman et al. 2022). Kenda et al. (2020) derived a subsurface model based on ground-motion compliance (the surface response to atmospheric pressure loading) measurements, which was extended by Onodera (2022), who proposed a model with a single discontinuity at ∼1 m depth in the first ∼75 m of the subsurface below the lander, corresponding to the martian regolith over a layer of coarse blocky material. In addition, Hobiger et al. (2021) proposed a subsurface model to explain the surface wave ellipticity features that they observed between 1 and 5 Hz. They propose a buried low velocity layer (bLVL) sandwiched between older and younger lava flows to explain a trough in the ellipticity curve. The proposed model spans down to ∼200 m depth but it lacks resolution in the first 20 m below the surface. A model with better resolution at shallow depths can be obtained by evaluating the H/V curve over a broader range of frequencies, so in this study we used the frequencies between 0.4 and 10 Hz.

After more than two years of operation, SEIS has recorded persistent ambient vibrations during different periods of time. In addition, more than 900 martian seismic events (hereafter, Marsquakes) have been identified by the Marsquakes Service (MQS, Clinton et al. 2018). In this work, we compute and compare the H/V ratios from both ambient vibrations and the coda waves of marsquakes. Then, based on the DFA, we invert the H/V curve to find the most representative shallow subsurface models.

2 THE INSIGHT DATA

2.1 Characteristics of the InSight data

The SEIS instrument records the martian ground motion at the InSight landing site with two sensors: a tri-axial very broad band (VBB) seismometer, with a higher sensitivity in the low-frequency (LF) range (especially 0.01–5 Hz), and a tri-axial short period (SP) seismometer, more sensitive to higher frequencies (>5 Hz), which together cover the frequency range between 0.01 to 50 Hz (Lognonné et al. 2019). Also, the TWINS (Temperature and Winds for InSight) provide wind and air temperature measurements at the InSight landing site to help in understanding the seismic recordings (Spiga et al. 2018; Banfield et al. 2019). All these sensors transmit data at different sampling rates depending on the requirements of the team and operational restrictions. In this sense, the most continuously available seismic data has been recorded by VBB at 20 sps (channel 02.BH). The availability of seismic data at 100 sps, either from VBB (00.HH) or SP (65.EH), is sparse and less continuous, mainly due to operational restrictions.

2.1.1 The ambient vibrations

Even though the 20 sps data would be enough to retrieve the H/V curve for the desired frequency range, the 100 sps data yields a clearer perspective of the whole spectral curve. Throughout the mission, simultaneous recordings by both VBB and SP at 100 sps spanning longer than one martian Sol are available only for two time periods: from Sol 422 to Sol 423, during the summer on the northern hemisphere, and from Sol 762 to Sol 763, during the northern winter. We analysed the seismic data recorded by both VBB and SP at 100 sps on Sol 422 (Fig. 1), which is a good representation of the period between spring and fall on the northern hemisphere, the least noisy time of the martian year. In contrast, data recorded on Sol 762 is representative of the period between northern fall and spring (including winter), when strong winds generate high noise levels throughout the Sol (Dahmen et al. 2021b).

Figure 1.

Spectrogram of martian ground motion recorded by VBB and SP at 100 sps on Sol 422. (a) Wind speed (grey circles) and air temperature (black lines) as a function of time. Data from Sols 422 and 423 are used to cover TWINS data gaps. (b) spectrogram of the vertical component of VBB acceleration at 100 sps. The median of horizontal (grey) and vertical (black) PSDs of VBB 100 sps data are computed between (c) 00:00–03:00 LMST, (d) 08:00–11:00 LMST and (e) 20:00–23:00 LMST. Plots (f)–(i) are analogous figures for SP 100 sps data.

For each instrument, we pre-processed the data by correcting for the instrument response, high-pass filtering at 0.01 Hz and rotating from the sensor orientations (U, V, W) to the geographic coordinate system (Z, N, E). When both VBB and SP data are simultaneously available, the spectral acceleration from a 120 s time window with 50 per cent overlap was computed for each component. The spectrograms of the vertical component of both VBB and SP are shown in Figs 1(b) and (f) along with the wind speed and atmospheric temperature data provided by TWINS (Fig. 1a).

The main characteristics of the ambient vibrations can be observed on Fig. 1. As shown by Giardini et al. (2020), the background noise level is highly affected by the wind and temperature conditions at the InSight landing site. The noisiest period corresponds with the highest wind speeds, which are recorded roughly between 08:00 LMST (Local Mean Solar Time) and 16:00 LMST. In contrast, the evening (between 16:00–24:00 LMST) is the quietest period of the day, where the lowest wind speeds are recorded (< 3 m s⁻¹). The very narrow peaks observed at frequencies below 10 Hz, especially during the evening and night, correspond to the so-called tick noise at 1 Hz and its overtones, which is an artefact generated by cross-talk in the measurement system. Further details are provided by Ceylan et al. (2021) and Zweifel et al. (2021).

Special attention must also be paid to the lander-related resonances, a set of major distinctive spectral peaks whose origin can be attributed to the lander system (Schimmel et al. 2021; Dahmen et al. 2021b). They have been studied in detail in the frequency range between 1 and 9 Hz, where they are characterized by a temperature-dependent peak frequency, a wind-sensitive amplitude, a predominantly horizontal polarization and are clearly excited at 1.6, 3.3, 4.1, 6.8 and 8.6 Hz (Ceylan et al. 2021; Dahmen et al. 2021b). Similar patterns are observed at higher frequencies (>10 Hz), so additional lander-related resonances are likely also present in this spectral range. Furthermore, Hurst et al. (2021) showed there are a number of resonances at about 2.86, 5.3, 9.5, 12, 14 and 23–28 Hz related to the Load Shunt Assembly (LSA), which provides mechanical separation between SEIS and the tether connecting to the lander in order to reduce the effect of lander perturbations on the seismometer.

Another remarkable feature of the ambient vibrations is the strong resonance at 2.4 Hz, which is more visible during the quiet times of the northern summer, roughly between 16:00–24:00 LMST. Dahmen et al. (2021b) investigated this mode and have shown its behaviour is contrary to the lander-related modes, as it is not excited by the steady winds during the mornings and it never disappears during quiet conditions (Fig. 1). Hobiger et al. (2021) proposed this mode is due to the response of a bLVL excited by Rayleigh waves generated by sources located at regional distances.

Besides the lander resonances, other features such as one-sided long-period pulses, termed glitches, can be observed in the data (Scholz et al. 2020). The description of other particular features in the seismic data recorded by SEIS are provided by Ceylan et al. (2021) and Kim et al. (2021).

2.1.2 The seismic events

As of 2021 September 30, the MQS has reported more than 900 potential seismic events (Clinton et al. 2021, see Data availability), whose origin and nature are under investigation, but a classification based on their spectral content can be made. These events are classified into two main groups: the LF family, characterized by long-period energy below 1 Hz and the high-frequency (HF) family, characterized by the excitation of frequencies mainly above 1 Hz. The LF family events can be divided into the LF events, which do not excite the 2.4 Hz mode, and the broad-band (BB) events, which do excite this mode (Clinton et al. 2021). On the other hand, in the HF family, we can distinguish three different groups of events: (1) the 2.4 Hz events, which excite the 2.4 Hz mode and frequencies up to 4 Hz; (2) the HF events, exciting the 2.4 Hz mode and higher frequencies up to 10 Hz and (3) the very high frequency (VF) events, which show a strong excitation of the horizontal components at frequencies above 5 Hz, up to as high as 35 Hz (Clinton et al. 2021).

Along with this classification based on the spectral energy content, the seismic events also have assigned a quality index based on the signal strength and the ability to identify and interpret the phase arrivals (Clinton et al. 2021). Therefore, the quality index can be either A (high), B (medium), C (low) or D (suspicious). As we make use of the coda of the marsquakes to compute the H/V curve, only events with quality A and B were used. In total, 139 marsquakes were used to compute the H/V curve, from which the 2.4 Hz events are the most abundant, as shown in Table 1. In the same way, the 100 sps channels were preferred whenever they were available, otherwise the 20 sps channels were used. Further details on the events used for this analysis, including the time windows for the H/V computation, can be found in the Supplementary Informtion.

Table 1.

Summary of martian seismic events recorded by SEIS and used in the H/V analysis, as reported by MQS through 2021 September 30. The marsquakes are differentiated by event type as classified by MQS: LF (low frequency), BB (broad band), 2.4HZ, HF (high frequency) and VF (very high frequency). The two main columns indicate: the number of events identified by the MQS per quality (A or B), and the number of events recorded either by VBB or SP at 20 or 100 sps.

	Events			Recordings
	QA	QB	Total	VBB@20sps	VBB@100sps	SP@20sps	SP@100sps
LF	4	9	13	13	2	3	5
BB	2	6	8	8	2	2	2
2.4 Hz	0	46	46	46	7	27	20
HF	0	51	51	51	12	27	22
VF	0	21	21	21	10	7	9
Total	6	133	139	139	33	66	58

The VF events are of special interest for this study as they contain energy well above 5 Hz, which allows mapping of shallower depths at the InSight site. Even though quality A VF events have not been reported by the MQS, there are 21 VF events with quality B (see Table 1). As shown in Table 1, VBB data at 100 sps are available for 10 VF events whereas SP data at 100 sps are available for nine events, with just one of them being recorded simultaneously by both SP and VBB.

3 H/V CURVES AT THE INSIGHT LANDING SITE

3.1 H/V analysis

3.1.1 The diffuse field approach

Based on Sánchez-Sesma et al. (2011), assuming a diffuse wavefield, the imaginary part of the GF between two sites is proportional to the average cross correlations of the corresponding displacements:

$$\begin{eqnarray*} \langle {u}_i\left( {{{{\bf x}}}_{\rm{A}},{\rm{\omega }}} \right)u_j^*\left( {{{{\bf x}}}_{\rm{B}},{\rm{\omega }}} \right)\rangle \propto {\rm{Im}}\left[ {{G}_{ij}\left( {{{{\bf x}}}_{\rm{B}},{{{\bf x}}}_{\rm{A}};{\rm{\omega }}} \right)} \right], \end{eqnarray*}$$

(1)

with |$\omega \ = \ 2\pi f$| the angular frequency, |${u}_i( {{{{\bf x}}}_{\rm{A}},\omega } )$| is the displacement field component in i direction at point |${{{\bf x}}}_{\rm{A}}$|⁠, the GF |${G}_{ij}( {{{{\bf x}}}_{\rm{B}},{{{\bf x}}}_{\rm{A}};\omega } )$| is defined as the displacement in the direction i at point |${{{\bf x}}}_{\rm{B}}$| due to the application of a unit harmonic point force in the direction j at point |${{{\bf x}}}_{\rm{A}}$|⁠. The asterisk (*) corresponds to the complex conjugate. If the sites |${{{\bf x}}}_{\rm{A}}$| and |${{{\bf x}}}_{\rm{B}}$| are the same, as well as the components i and j, then the average of the autocorrelation corresponds to the directional energy density |${E}_i( {{{\bf x}},\omega } )$| at point x (Sánchez-Sesma et al. 2008; Perton et al. 2009), such that:

$$\begin{eqnarray*} {E}_i\ \left( {{{\bf x}},\omega } \right) = {\rm{\ }}\rho {\omega }^2 \left\langle \|{u}_i\left( {{{\bf x}},\omega } \right){\|}^2 \right\rangle \propto {\rm{Im}}\left[ {{G}_{ii}\left( {{{\bf x}},{{\bf x}};\omega } \right)} \right], \end{eqnarray*}$$

(2)

where ρ is the mass density. Thus, the H/V spectral ratio in a diffuse field, defined as the square root of the ratio between the horizontal and the vertical energies (Arai & Tokimatsu 2004), is equivalent to the square root of the ratio between the imaginary parts of the corresponding GFs, that is,

$$\begin{eqnarray*} \mathrm{ HVSR}\ = {\rm{\ }}\frac{H}{V}\left( {{{\bf x}},\omega } \right) = \sqrt {\frac{{\langle \|{u}_1\left( {{{\bf x}},\omega } \right){\|}^2\rangle + \langle \|{u}_2\left( {{{\bf x}},\omega } \right){\|}^2\rangle }}{{\langle \|{u}_3\left( {{{\bf x}},\omega } \right){\|}^2\rangle }}} \ , \end{eqnarray*}$$

(3)

$$\begin{eqnarray*} = \sqrt {\frac{{{\rm{Im}}\left[ {{G}_{11}\left( {{{\bf x}},{{\bf x}};\omega } \right)} \right] + {\rm{Im}}\left[ {{G}_{22}\left( {{{\bf x}},{{\bf x}};\omega } \right)} \right]}}{{{\rm{Im}}\left[ {{G}_{33}\left( {{{\bf x}},{{\bf x}};\omega } \right)} \right]}}} \ , \end{eqnarray*}$$

(4)

where i = 1, 2 correspond to the horizontal components and i = 3 is the vertical component.

We computed the H/V spectral ratios for both the ambient seismic noise data (nHV) and the seismic events (eHV) following eq. (3), thus adhering to the diffuse field approach. Our methodology is slightly different from the classical average of H/V spectral curves as we compute the final H/V from the average spectrum of each component. Because of this, to have the same spectral contribution from each time window, we normalized the vertical and horizontals spectra by the total power of the respective time window. The normalized ground motion corresponds to:

$$\begin{eqnarray*} {\hat{u}}_i\ \left( {{{\bf x}},\omega } \right) = \frac{{{u}_i\left( {{{\bf x}},\omega } \right)}}{{\sqrt {\mathop \sum \nolimits_{i = 1}^3 {{\left| {{u}_i\left( {{{\bf x}},{\rm{\Delta }}\omega } \right)} \right|}}^2} }}\ , \end{eqnarray*}$$

(5)

where the denominator is the normalization term, given by the total power in a certain bandwidth ∆ω. For simplicity, we used ∆ω = 50 Hz to normalize by the power of the total spectrum. This way, the final H/V is computed as in eq. (3) but replacing |${u}_i$| by the normalized ground motion |${\hat{u}}_i$|⁠.

For the computation of nHV, we used the data recorded during Sol 422 and Sol 423, which were divided into three different groups of data representing three different periods of time along a typical summer-time Sol (Fig. 1). These are a high-noise period (between 08:00–11:00 LMST), a low-noise period (20:00–23:00 LMST) and a moderate-noise period (00:00–03:00 LMST). For each instrument, SP and VBB, the H/V was computed using time windows of 120 s length with an overlap of 50 per cent. In order to compare the performance of VBB against SP, we used the same 120 s data windows such that the segment is skipped when there is no data available for one of the sensors.

The H/V of the marsquakes data was obtained in a similar way but based on the phase picks provided by the MQS. For events with quality A and B, the MQS provides the times of the P- and S-wave onsets and also includes the times of the characteristic background seismic noise before the event (Clinton et al. 2021). According to this catalogue, we computed the H/V of the S-wave coda, which has been shown to be depolarized and the energy to be equipartitioned (Menina et al. 2021; van Driel et al. 2021). The duration of this coda window is variable and is provided by the MQS such that it is a good representation of the S-wave coda and free of non-seismic perturbations (e.g. glitches). The H/V computation was performed in a similar way as for the noise data but using consecutive 20 s data windows with an overlap of 50 per cent.

For each event and sensor, the vertical spectral signal-to-noise ratio (SNR_Z) was computed using the median power of the vertical component of the coda and the pre-event noise windows. For each sensor, a single eHV curve is derived as the median of all the segments with SNR_Z larger than 2 dB. The DFA-HV curve is expected to be less smoothed than the classical H/V approach, which typically implements a Konno–Ohmachi smoothing window (Konno & Ohmachi 1998). Because of this, a visual inspection and manual smoothing following the median of the final eHV is performed in order to also discard any contamination due to lander resonances.

3.1.2 Polarization attributes

In order to have a better assessment of the spectral features of the H/V curves, we also computed the polarization attributes of both the seismic noise and the seismic events using a similar approach as presented by Park et al. (1987), which is based on the work of Samson (1983). For each data window, the calculation of the polarization features relies on the eigen-decomposition of the spectral covariance matrix (S) of the three-component ground motion data. The largest eigenvalue and the corresponding eigenvector of S represent the polarization of the seismic signal (either noise or event coda) during the given time window for each frequency. The directional attributes of this polarization are of particular interest: the horizontal incidence angle θ_H, measured in degrees clockwise from North, and the vertical incidence angle θ_V, restricted to lie between 0° (vertical) and 90° (horizontal).

We also computed the degree of polarization (dop), as proposed by Samson & Olson (1980), which is a measure of how well the seismic signal can be described by fewer than three degrees of freedom. For a three-component ground motion record, the degree of polarization is defined as

$$\begin{eqnarray*} dop\ = \frac{{3{\rm{tr}}\left( {{{{\bf S}}}^2} \right) - {{\left[ {{\rm{tr}}\left( {{\bf S}} \right)} \right]}}^2}}{{2{{\left[ {{\rm{tr}}\left( {{\bf S}} \right)} \right]}}^2}}\ , \end{eqnarray*}$$

(6)

where tr is the trace and S is the spectral covariance matrix of the ground motion record. This parameter ranges from 0, when the signal is depolarized, to 1, when a single non-zero eigenvalue exists and the signal is completely polarized in the direction of the corresponding eigenvector. High values of dop correspond to a polarized signal with a particle motion that can be either linear or elliptical. As lander resonances have been observed to be highly polarized (e.g. Schimmel et al. 2021; Dahmen et al. 2021b), the polarization attributes will be used for a better assessment and understanding of the resulting H/V features.

3.2 Ambient noise H/V (nHV)

From the seismic data simultaneously recorded by VBB and SP during Sols 422 and 423, a total of 1576 segments of 120 s length were used for our analysis. From these data, 215 segments correspond to the low-noise period (between 20:00 and 23:00 LMST), 202 segments to the moderate-noise period (00:00–03:00 LMST) and 216 segments to the high-noise period (08:00–11:00 LMST). The spectral (vertical, horizontal and H/V spectra) and polarization features (dop, θ_H and θ_V) are illustrated in Fig. 2.

Figure 2.

Spectral acceleration and polarization features during three different time periods of Sols 422 and 423 (northern summer). (a)–(f) correspond to the vertical spectral acceleration, the horizontal spectral acceleration, the HVSR, the degree of polarization (dop) computed from VBB data, the horizontal (θ_H) and vertical (θ_V) incidence angle, respectively, of the time period between 20:00 and 23:00 LMST. Analogously, plots (g)–(l) correspond to the time period between 00:00 and 03:00 LMST (moderate noise) whereas plots (m)–(r) are the respective features during the time period between 08:00 and 11:00 LMST (high noise). The plots (a)–(c), (g)–(i) and (m)–(o) show the median of the VBB (black) and SP (grey) data. Dashed red curve corresponds to the theoretical self-noise of VBB.

The spectral curves from both VBB and SP data (Fig. 2) reveal clear amplitude differences in the vertical components, at different noise periods, but an overall good match between the VBB and SP horizontal components.

3.2.1 Low-noise period

During the low-noise period (2a–f), the tick noise is a set of clear and narrow peaks observed at frequencies below 10 Hz, especially in the SP data. When no smoothing is applied, their presence in the H/V curve can be clearly identified and discarded. There is an overall good match between the SP and VBB horizontal (VBBH) components for frequencies above 6 up to 21 Hz, which, in the vertical component is only observed between 14 and 20 Hz. Between 6 and 14 Hz, the amplitudes recorded by SP are slightly lower than those recorded by VBB, while below 6 Hz, the VBB recordings have lower amplitudes than SP. These amplitude differences below 14 Hz have been previously observed by Dahmen et al. (2021b), who interpreted the VBB spectral amplitudes between 6 and 14 Hz as its spectral noise floor, close to the theoretical instrument noise (Lognonné et al. 2019; Stutzmann et al. 2021), as shown in Fig. 2(a). The same is observed for the SP amplitudes below 6 Hz, where the instrument noise level is reached. Besides the 2.4 Hz trough clearly visible from the VBB data, the resulting H/V curves are flat for frequencies below 10 Hz.

Above 20 Hz, the recorded VBBH amplitudes are slightly higher than those from the SP (SPH), whereas the opposite is true for the vertical components (VBBZ and SPZ). As the horizontal amplitudes are higher than the vertical ones for both sensors (see Fig. 1), a broad peak in the H/V curve is observed, with higher amplitudes for the VBB sensor, particularly in the frequency range between 20 and 30 Hz. In this frequency range, singular H/V peaks related to highly polarized ground motion (dop > 0.7) are revealed (Fig. 1d). The same is observed for other peaks at 10–12 Hz and around 15 Hz. All these peaks have a nearly horizontal incidence (θ_V > 60°), with different incidence angles θ_H. Although the low-noise period is related to mainly low-wind velocities, high-wind bursts may also occur, which would explain the high polarization and horizontal incidence of the peaks observed at high frequencies (>10 Hz).

3.2.2 Moderate- and high-noise period

During the moderate- and high-noise period (Figs 2g–l and m–r), the ambient vibrations level increases due to stronger winds and the aforementioned high frequency peaks become even more distinguishable from the background noise. The overall spectral amplitudes are higher than in the low-noise period and, above 10 Hz, VBBZ has a persistent lower amplitude than SPZ. A variety of singular peaks are clearly observed on the H/V curve, including the lander modes in the ranges 3–3.3 Hz, 3.8–4.1 Hz and around 6.8 Hz, which reveal a characteristic high and nearly horizontal polarization. All the singular peaks at higher frequencies (e.g. ∼11, ∼15, ∼20 and ∼25 Hz) are even more polarized than in the low-noise period, with a predominantly horizontal rather than vertical incidence angle θ_V, similar to the 4.1 Hz resonance, a well-known lander resonance (Dahmen et al. 2021b). Even more notable than in the low-noise period, striking differences between VBB and SP are observed for the H/V curves at frequencies above 10 Hz.

3.3 Marsquakes H/V (eHV)

The eHV curve was computed from the coda of 139 marsquakes. As we are interested in the computation of the H/V for actual seismic coda signal, we used the corresponding SNR_Z to choose the spectral segments with seismic energy above the pre-event seismic noise. Each type of marsquake excites different frequencies (Fig. 3), so using a pre-defined frequency range for different event types is meaningless, because there might be frequencies that are not excited by some seismic events. Thus, we combine all the H/V curves sections with SNR_Z > 2 dB.

Figure 3.

Spectral features of different classes of marsquakes. From top to bottom: LF (S0173a), BB (S0235b), 2.4 Hz (S0372a), HF (S0992a) and VF (S0794a and S0756a) class. Data from VBB at 20 sps (02.BH) illustrate the characteristics of the LF, BB, 2.4 Hz and HF classes. The VF class is represented by two events independently recorded by VBB (00.HH) and SP (65.EH) at 100 sps. For each event, the SNR_Z is plotted as a function of time, where a characteristic noise and S-wave coda window are indicated as grey and black boxes, respectively. In the second column, the vertical (solid) and horizontal (dashed) spectral power of the S-wave coda signal is compared to the characteristic noise window. The H/V spectral ratio computed under the DFA approach is represented by the black solid line. The shaded rectangular areas indicate the frequencies covered by the different types of events.

In order to check the SP and VBB differences revealed by the noise data, we also show two strong VF events independently recorded by VBB (event S0794a) and SP (event S0756a). Both events carry seismic energy above 1 Hz, but especially above 10 Hz (Fig. 3). They share H/V similarities between 1 and 10 Hz, but there are striking differences above 10 Hz. The H/V obtained from the VBB recording presents a plateau-like shape followed by a strong peak at around 25 Hz, whereas the SP recording reveals a decay towards higher frequencies.

The origin of the VF events, which shape the H/V curve above 5 Hz, remains unknown, so the observed waveforms and different H/V curves might be influenced by the environmental conditions. Nevertheless, during the occurrence of both VF events, no major environmental changes related to the first arrival of the seismic energy or to the proposed S-wave arrival are observed (Fig. 4). We thus conclude that the H/V curves of the VF events are unrelated to the environmental conditions. Also, as we are using the seismic coda recorded at great distances, we exclude the source as a determining feature.

Figure 4.

Environmental conditions during two large VF events: S0756a (left) and S0794a (right). From top to bottom: air temperature, wind speed, wind direction, vertical acceleration and horizontal acceleration. The vertical (black) and horizontal (grey) ground-motion recordings are bandpass filtered between 6 and 9 Hz. The vertical yellow line indicates the S-wave arrival time for each event.

We computed the median H/V curve for both SP and VBB sensors, separately (Fig. 5). As the DFA-HV computation leads to unsmoothed curves, we applied a manual smoothing by following the resulting median H/V curve, as shown in Fig. 5, which allowed us to mitigate the effect of lander resonances such as the 4.1 Hz. Due to the lack of events with enough energy above the noise, the eHV curve for the SP sensor is unreliable for frequencies below 2 Hz. Therefore, given the target frequency range between 0.4 and 10 Hz, the preferred eHV curve is obtained from the VBB recordings.

Figure 5.

H/V curves derived from marsquakes recorded by VBB (left) and SP (right) with vertical SNR larger than 2 dB. The curves are coloured corresponding to different types of marsquakes. The median and smoothed median H/V are denoted by the continuous grey and dashed black lines. The grey areas illustrate the smoothed median ± 15 per cent uncertainty range. The bottom plots indicate the amount of marsquakes with SNR_Z > 2 dB per frequency (N_f). Separate H/V curves for each type of marsquakes are also plotted with the corresponding colour.

Although there are only a few events with energy well above the ambient noise for frequencies above 10 Hz, the H/V differences between SP and VBB are still notable, similarly to what was derived from the noise data. Thus, we infer this phenomenon is most likely related to the instruments themselves. The reason behind this SP-VBB difference is still unclear but seems to be related to the amplification of the vertical components above 10 Hz. When compared to the horizontals, the VBBZ seems to deamplify the ground motion in the frequency range roughly between 10–35 Hz (see Fig. 1), whereas the SP sensor reveals similar amplitudes between SPZ and SPH, especially for frequencies above 5 Hz, where the instrument response is flat (Lognonné et al. 2019). The SP sensor was designed to record ground motion at HFs, so it would be preferable here. Nevertheless, as the VBB and SP sensors are placed in different spots on the SEIS Leveling system (LVL, Lognonné et al. 2019), a possible reason for the SP recording larger vertical ground motions is that it could be exposed to larger vertical motions due to rotation of SEIS (Fayon et al. 2018). Due to the persistent incompatibility of VBB and SP recordings at high frequencies (>10 Hz), both from noise or marsquakes data, we dismissed using the H/V curves at high frequencies for any further analysis.

Finally, for frequencies below 10 Hz, the VBB H/V curve is chosen to be the representative eHV curve for the InSight landing site, which has a characteristic trough around 2.4 Hz and a strong peak around 8 Hz.

3.4 nHV versus eHV

Clear differences are observed between the VBB nHV and eHV curves, as shown in Fig. 6. On Earth, despite some differences, a generally good agreement between the nHV and eHV curves is typically observed (e.g. Parolai et al. 2004; Pilz et al. 2009; Rivet et al. 2015). Napolitano et al. (2018) suggest the differences between nHV and eHV curves can be related to an effect of the local topography, however, the InSight landing site is mostly flat without complex topography in the surrounding. Hence, the difference between the recorded eHV and nHV curves can be hardly related to a topographical effect. Rather, the reason for the difference between nHV and eHV is different for the quiet and noisy periods.

Figure 6.

H/V curves for three different contexts at the InSight landing site as recorded by VBB (black) and SP (grey). (a) Low noise: mainly recorded during the martian evening and night, with wind speeds lower than ∼ 2.8 m s⁻¹, (b) high noise: typically recorded during the martian day-time, when high-wind speeds (> 5 m s⁻¹) are recorded and (c) events: collection of early S-wave coda of the best recorded marsquakes.

In the low-noise period, the nHV curve is similar to the eHV curve at frequencies below ∼4 Hz, including the trough around 2.4 Hz. Above 4 Hz, the eHV presents a strong peak around 8 Hz that is not visible on the low-noise nHV curve. We consider this difference might be caused by VBB most likely recording the instrumental self-noise at frequencies between 4 and 14 Hz. These observations suggest that, at least during the quietest period of a martian Sol, there is either an absence of energy coming from seismic waves propagating through the InSight landing site or the corresponding seismic energy is lower than the instrument self-noise of the SEIS instruments. Thus, due to the proximity and shape similarity of VBBZ to the theoretical self-noise, we suspect there is a risk during this time period of measuring the corresponding self-noises and thus the obtained H/V curve might be not representative of the subsurface structure, as has been experienced in some Earth cases (e.g. van Ginkel et al. 2020).

During the quietest period of a martian Sol, the potential sources of noise at the InSight landing site are minimized due to the low-wind speeds (Mimoun et al. 2017; Murdoch et al. 2017), so the instrument self-noise is likely recorded either on the horizontal or vertical components. Velocity variations in the frequency range between 6 and 9 Hz have been observed especially during the quiet periods. These variations are evidenced only on the horizontal (north and east) but not on the vertical component and can be understood as the thermo-elastic perturbation of the martian regolith sensed by a seismic wavefield (Compaire et al. 2022). Consequently, even if the horizontal ground motion recorded during the quiet period is generated by the propagating seismic wavefield, it is unreliable to retrieve any information related to the subsurface structure from the H/V curve, because the vertical component is likely dominated by the instrument self-noise.

During moderate- and high-noise periods, the differences between the eHV and the nHV curves are caused by different factors. The high frequency signals recorded during these windy times, especially in the high-noise period, are above the instrumental noise, so they can be interpreted as actual ground-motion signals. The strong winds excite frequencies around 3.3, 4.1 and 6.9 Hz, which coincide with the lander resonances investigated by Dahmen et al. (2021b). Stronger excitation of frequencies around 11, 15, 20 and 25 Hz, which likely have a mechanical origin due to the interaction between the wind and the lander, are also observed. In the high-noise period, many of the narrow peaks are also related to the lander resonances, for example, 3.0, 3.8 and 5.8 Hz. These are the same lander resonances observed during the moderate-noise period, but their eigenfrequencies are shifted because of the daily variation due to temperature changes (Dahmen et al. 2021b). A similar behaviour is observed at higher frequencies, where the frequencies around 11, 15, 20 and 25 Hz are strongly excited and shifted.

The excitation of these non-natural resonances contaminates the H/V curves due to their mechanical origin, and thus do not represent the local subsurface structure. In addition, the broad-band shape of the moderate- and high-noise nHV curves can be understood as the mechanical noise generated by the interaction of strong martian winds with the InSight lander and the WTS, as proposed by Murdoch et al. (2017). In this case, even at high frequencies, the horizontal ground response to the wind is expected to be greater than the vertical one due to propagation of shear stress (Naderyan et al. 2016; Lott et al. 2017). This broad-band nHV shape is not directly related to the subsurface structure. Indeed, Mucciarelli et al. (2005) showed that high wind speeds enhance the horizontal components, likely explained by flow rotation near the soil and the instrument, leading to unreliable H/V peaks at high frequencies.

In general, the influence of the lander is highly visible on the nHV curves during moderate- and high-noise periods. Between 1 and 10 Hz, even though these nHV curves also show an increasing towards higher frequencies (Fig. 6b), its shape is different (less steep) than the preferred eHV curve (Fig. 6c). For example, the H/V amplitude at 8 Hz is lower on the high-noise nHV (∼2) than the eHV (∼8). Also, the effect of Sorrels propagative pressure waves is identifiable as a broad trough below 1 Hz (Kenda et al. 2020; Stutzmann et al. 2021). Because of this high contamination by external factors, we discarded these noisy periods from the H/V analysis.

Thus, the VBB eHV curve is preferred for the inversion for the subsurface structure instead of any of the nHV curves. In particular, the VBB eHV is suitable for the H/V inversion under the DFA because of the high diffusivity and strong multiple-scattering revealed for the martian crust (Lognonné et al. 2020; Menina et al. 2021; van Driel et al. 2021). These phenomena lead to a diffuse seismic coda, where a combination of scattered body and surface waves is expected and therefore the full wavefield DFA is a good representation.

4 H/V INVERSION

4.1 H/V inversion using DFA

The eHV curve (hereafter, the final H/V and defined as fHV) is chosen to be a good representation of the H/V response of the subsurface structure at the InSight landing site. We inverted this fHV curve for a 1-D shear wave velocity model (V_s) using the diffuse wavefield approach. The inversion scheme can be divided into two parts: the theoretical H/V forward modelling and the nonlinear inversion procedure.

4.1.1 Forward modelling of H/V

Given a 1-D layered subsurface model (including P- and S-wave velocities, density and thickness of each layer), the theoretical H/V spectral curve under the DFA can be computed from the corresponding GF, following Eq. (4). In this work, the GFs are calculated using the technique developed by Margerin (2009), which consists on the summation of the generalized eigenfunctions of the layered elastic medium. This technique is implemented via an open-source code that has been already successfully applied to terrestrial data (e.g. Margerin et al. 2009; Galluzzo et al. 2015).

The code calculates the vertical and horizontal spectral ground-motion for each P, S, Love and Rayleigh waves individually, so the contribution of body and surface waves can be separately derived. The H/V curve for the full wavefield, including the contribution of all the types of waves, can be also obtained. Given the diffuse nature of the seismic coda used for the eHV curve, instead of using one type of wave propagation mode such as fundamental mode Rayleigh waves as in Hobiger et al. (2021), we make use of the full wavefield H/V curve as the forward theoretical prediction to explain the interpreted signal of the eHV curve at the Insight landing site.

4.1.2 Inversion algorithm and parametrization

In order to find the best subsurface models explaining the observed fHV curve, we employed the Neighbourhood Algorithm (NA), which was developed by Sambridge (1999) and adjusted for this case by Wathelet (2008). This technique was implemented via the dinver toolbox included in the Geopsy package (Wathelet et al. 2020). The NA is suitable for this type of inversion due to the intrinsic non-uniqueness problem of the H/V curves, that is, without further information or constraints, the same H/V curve can be explained by many different subsurface models (e.g. Piña-Flores et al. 2016). The NA provides an efficient exploration of the parameters space to retrieve different subsurface models fitting the H/V curve.

The NA makes use of the misfit between observed and theoretical H/V curves to explore the area of the parameters space where the models with low misfit are expected. The misfit function is defined as follows:

$$\begin{eqnarray*} {\phi }_{{\rm{hv}}}\ \left( {{\bf m}} \right) = \sqrt {\frac{1}{N}\mathop \sum \limits_{i = 1}^N \frac{{{{\left( {\log \left( {{{\frac{H}{V}}}_{{\rm{m}}i}} \right) - \log \left( {{{\frac{H}{V}}}_{{\rm{d}}i}} \right)} \right)}}^2}}{{\sigma _{{\rm{d}}i}^2}}} \ \end{eqnarray*}$$

(7)

where m is the evaluated model, N is the number of frequency points of the observed (⁠|$\frac{H}{{{V}_{\rm{d}}}}$|⁠) and full-wavefield modelled (⁠|$\frac{H}{{{V}_{\rm{m}}}}$|⁠) curves. |${\sigma }_{{\rm{di}}}$| is the standard deviation of the |$\frac{H}{{{V}_{\rm{d}}}}$| curve at a given frequency. Following Lontsi et al. (2015), we implemented a logarithmic misfit to avoid overweighting strong peaks in the H/V curve. After the parameters space is explored with the NA, for a given parameters setting, several best models with the misfit around the global minimum are selected as the possible models explaining the observed fHV curve.

We inverted for the shear wave velocity structure only, including the corresponding layer thickness, as it has been demonstrated to be the main parameter affecting the H/V curve (e.g. Malischewsky & Scherbaum 2004; Bonnefoy-Claudet et al. 2006). For the inversion, the parameters space explored by the NA is bounded by the prior knowledge based on the pre-landing model. This pre-landing model, based on the size range of rocky ejecta craters and exposures in nearby fracture zones, suggests a transition from the shallow regolith to a coarse ejecta layer, followed by a fractured basaltic layer and a deeper pristine basaltic unit (Golombek et al. 2017; Knapmeyer-Endrun et al. 2017; Warner et al. 2017). Based on the size range of the rocky-ejecta craters and ghost craters, previous mapping suggests that the basalt layer has a thickness of ∼170 m thick near the landing area (Warner et al. 2017; Golombek et al. 2018; Pan et al. 2020) and are underlain by sedimentary rocks likely of Noachian age (Pan et al. 2020).

Thus, based on the pre-landing model for the first 50 m, the parameters space is defined for a model with four layers over a half-space (4LOH). The results of the inversion are affected by this constraint such that a different number of layers will likely lead to different solutions. Table 2 presents the value ranges of the parameter spaces explored in the different inversion runs. Further details on these values are explained in the following.

Table 2.

Range of values defining the parameter spaces to be explored during the inversion of the fHV curve. The values in parenthesis correspond to the maximum shear-wave velocity for the small (S) and medium (M) parameter spaces, respectively. Density and V_p/V_s are fixed for each layer.

	Bottom depth (m)	VS (m s−1)	Density (kg m−3-)	Vp/Vs
Layer 1	0–10	Top: 40–80 Bottom: 70–300	1350	1.87
Layer 2	3–30	100–2000 (1000, 1500)	2000	1.73
Layer 3	5–60	100–2000 (1000, 1500)	2900	1.73
Layer 4	10–100	100–2000 (1000, 1500)	2900	1.73
Half-space	–	200–2000 (1000, 1500)	3300	1.73

The inversion of the H/V curve alone, without further observations, has been demonstrated to be heavily affected by the non-uniqueness phenomenon (e.g. Piña-Flores et al. 2016). Therefore, constraining the parameters space is relevant for obtaining the proper model. So far, the best new constraints of the inversion are the shear wave velocities obtained from the hammering experiment (Lognonné et al. 2020; Brinkman et al. 2022), determined to be lower than 80 m s⁻¹ at the surface, and the thickness of the shallow regolith layer, expected to be no deeper than 10 m (Golombek et al. 2020b,a). Due to the expected compaction, the S-wave velocity of the first layer is modelled as a power law increase. Because the H/V inversion is not very sensitive to variations in density, this parameter is fixed for each layer, as well as the V_p/V_s ratio. This way, the P-wave velocity is derived from the inverted S-wave velocity. In particular, for the shallowest layer, we used V_p/V_s = 1.87, as reported by Brinkman et al. (2022), and ρ = 1350 kg m⁻³, following the pre-landing assessment of this site (Morgan et al. 2018). This V_p/V_s lies in the range 1.68–2.56 for dry unconsolidated material, as shown by investigations on Earth samples and also derived for the Lunar regolith (e.g. Uyanik 2010; Dal Moro 2015). The density and V_p/V_s of the deeper layers are set to standard values of volcanic regions on Earth (e.g. Lesage et al. 2018).

Regarding the deeper layers, in volcanic areas on Earth, the shear wave velocities of basaltic layers have been shown to range from 500 to 2000 m s⁻¹ at shallow depths (Delage et al. 2017). On the other hand, an unexpected low-velocity layer should be placed at larger depths to explain the 2.4 Hz resonance (Hobiger et al. 2021). We explored a large parameters space with shear wave velocities between 100 and 2000 m s⁻¹ extending down to 100 m. In particular, we performed different inversion runs with different maximum shear wave velocities of 1000 (small), 1500 (medium) and 2000 m s⁻¹ (large model).

As shown in Fig. 6, our target fHV has a narrow trough around 2.4 Hz and a steep left flank leading to a strong peak at around 8 Hz. No clear peaks or troughs are observed at low frequencies between 0.4 and 2 Hz. For each frequency, a common uncertainty range of 15 per cent of the H/V amplitude was set.

For each 4LOH parameters space presented in Table 2 (S1, M1, L1), a total of ∼45 000 models were explored using the NA, for which three different inversion runs were performed in order to have a better search of the models fitting the fHV curve. For each parameters space, the best fitting models with a misfit lower than 0.3 were chosen, from which 200 random models are plotted in Fig. 7. This way, three different subsets of models that fit the fHV curve well can be identified. The best misfit of these models is obtained from the M1 inversion runs (∼0.242). However, an overall good fitting of the fHV curve is obtained for models with misfit values lower than 0.3 and thus all these models can be a good representation of the subsurface structure.

Figure 7.

Sets of models obtained from the inversion of the fHV curve, between 0.4 and 10 Hz, using different parameter spaces. 200 best random models are plotted for each of the three parameters space explored: (a) S1 (small, orange), (b) M1 (medium, green) and (c) L1 (large, blue). For each parameters space, the models with darkest colours indicate a best (lower) misfit. (d) Synthetic full wavefield H/V response of the different models drawn in subfigures (a)–(c). Two representative models from each parameters space are plotted following the corresponding colour pattern. Besides, the constrained MAP and ML models from Hobiger et al. (2021) and the model derived from ground compliance measurements (Onodera 2022), together with their corresponding synthetic full-wavefield H/V response, are plotted as magenta, yellow and red dashed lines, respectively.

5 DISCUSSION

5.1 Implications for the subsurface structure

5.1.1 Models from fHV inversion

Some similarities and differences can be highlighted from the comparison of the S1, M1 and L1 models (the SML1 models). They all present a characteristic shallow high velocity layer (SHVL) followed by a bLVL with a sharp contrast between the SHVL and the bLVL. The latter bLVL has been proposed to be the cause of the 2.4 Hz resonance due to seismic waves trapped in this low-velocity layer (Hobiger et al. 2021). Against the pre-landing subsurface models, the bLVL in the resulting models of this work, as well as Hobiger et al. (2021), is an unexpected feature of the subsurface structure. A similar sharp contrast is expected between the SHVL and the subsurface structure above, which can be understood as the mechanism to generate the 8 Hz peak on the eHV curve. The bLVL is followed by a faster half-space, which has lower velocities than the SHVL and therefore represents a less sharp discontinuity.

The bLVL of the L1 models is in good agreement with the bLVL proposed by Hobiger et al. (2021, hereafter, HOB) to explain the 2.4 Hz trough. Nevertheless, in our case, the units below the bLVL have relative lower velocities than the SHVL whereas, in the HOB models, the bLVL is followed by a layer with velocities even higher than the SHVL, which generates an H/V peak between 1 and 2 Hz that is not observed from the SEIS data, as shown in Fig. 7(d).

The models arrive at different velocities and thicknesses of each layer, which is a demonstration of the non-uniqueness problem of the H/V inversion. In general, depending on the thickness and velocity of the SHVL, the bLVL will be placed at shallower or larger depths with a thickness between 15 and 25 m and a velocity ranging between 100 and 500 m s⁻¹. In any case, a generic subsurface pattern of a bLVL below a SHVL is required to simultaneously generate a H/V peak around 8 Hz and trough around 2.4 Hz.

Our results suggest the features of the fHV curve can be explained by the full wavefield rather than by a single type of waves. In particular, the 2.4 Hz resonance is hardly derived from the fundamental Rayleigh mode as computed under the DFA (see Fig. S1, Supporting Information). Instead, this resonance is closely coinciding with a P-wave resonance and modelled as a resonance of the full wavefield.

For the L1 models, the bLVL coincides with what has been proposed to be a geological unit sandwiched between Early Hesperian to Early Amazonian basalt flows, which may have a sedimentary origin (Hobiger et al. 2021). The SHVL then corresponds to the Early Amazonian basaltic lava flows (∼1.7 Ga, Warner et al. 2017; Wilson et al. 2019). Meanwhile, the half-space below the bLVL must be either a layer within the Early Amazonian or the Early Hesperian basalt lava flows (Warner et al. 2017; Hobiger et al. 2021).

Crater statistics suggest significant erosion and landscape degradation between an Early Hesperian volcanic surface (3.6 Ga) and the early Amazonian effusive volcanism (1.7 Ga), which is estimated to be ∼140 m thick (Warner et al. 2022). Furthermore, throughout the dichotomy boundary between the southern Noachian highlands and the northern plains, there are Noachian through Hesperian transition units that indicate active erosion and deposition of sedimentary materials (Tanaka et al. 2014; Pan et al. 2017, 2020). To the south of InSight lander, some of the Amazonian–Hesperian transition units are sedimentary deposits (Pan et al. 2020) as well as the Medusae Fossae Formation (Tanaka et al. 2014). In addition, alluvial activity has occurred further south in Gale crater during this period (Grant et al. 2014; Grant & Wilson 2019). Given this evidence for aqueous activity in the geological record, it is possible that there was both erosion of the Early Hesperian surface and deposition of sedimentary deposits prior to the deposition of the Early Amazonian basalt flows.

If Early Hesperian lava flows are beneath the bLVL, these rocks would be expected to be more fractured by impact cratering due to the higher cratering rate of the Hesperian (e.g. Hartmann & Neukum 2001) or longer exposure to the surface; and to be more altered by aqueous activity than the Early Amazonian volcanics. The expected greater physical and chemical weathering of the basalt flows in Early Hesperian presents a reasonable explanation for the relatively lower seismic velocities of the half-space beneath bLVL compared to the SHVL.

In the case of the S1 and M1 models, the bLVL is located at shallower depths (< 40 m) and the corresponding SHVL right above it is a thin and fast layer. Although we cannot completely discard these models, the thin SHVL and the extremely low velocities of the bLVL, considering the effects of the overburden pressure, seem to be unlikely and disfavour these models. Nevertheless, if the SHVL had higher velocities, its thickness would also increase and therefore the bLVL would be placed at larger depths. In any case, the thickness of the bLVL would range between 15–25 m, which is less than the ∼50 m thickness proposed by Hobiger et al. (2021).

There is a lack of resolution for the structure above the SHVL because the higher frequencies were excluded from the analysis. Nevertheless, the inversion of compliance measurements suggests the regolith layer is most likely about 1 m thick and everything below is associated with coarse blocky ejecta material (Onodera 2022). This shallow interface would generate an H/V peak between 25 and 30 Hz, which is out of the frequency range here investigated. A H/V peak at these frequencies is hardly observed from the data due to the strong excitation of a highly polarized resonance at this frequency range (see Figs 1 and 2), which might be related to the LSA (Hurst et al. 2021). Anyway, this shallow discontinuity between the regolith and the coarse ejecta layer cannot be discarded. Therefore, we can infer the top of the SHVL is likely below the coarse ejecta layer and does not correspond to the regolith interface.

5.1.2 Alternative model without 2.4 Hz trough

The inversions above were performed assuming the origin of the 2.4 Hz resonance is related to the subsurface structure. In fact, this has not yet been proven and the mechanism to generate this resonance remains unclear. An interesting observation is that this trough is clearly observed from seismic noise data and also from seismic events data, which is in contradiction to previous observations on Earth, where the effect of a bLVL on the H/V curve has been observed to be less visible in the noise data than the events (Di Giacomo et al. 2005; Panzera et al. 2015).

Because of this, we performed a second round of fHV inversions using the same parametrization but simultaneously fitting the frequency ranges between 0.4–2 and 3–10 Hz, so the 2.4 Hz is ignored. Once again, three parameter spaces are explored and are analogously referred as S2, M2 and L2 models (SML2) and illustrated in Fig. 8 in the same way as for the SML1 models (Fig. 7).

Figure 8.

Sets of models obtained from the inversion of the fHV curve in two different frequency ranges, between 0.4–2 and 3–10 Hz, using different parameter spaces. 200 best random models are plotted for each of the three parameters space explored. The colour patterns and other features are analogous to Fig. 7.

From these new inversion runs, the non-uniqueness problem is still apparent, so it is difficult to provide one single model for the subsurface structure. One can see the S2, M2 and L2 models are analogous to the SML1 models. However, some interesting features can be highlighted. In both cases, ignoring the 2.4 Hz trough or not, the SHVL is required for modelling the strong 8 Hz peak and low velocities below it are needed. Nevertheless, in the SML2 models, the velocity contrast between the SHVL and the underlying bLVL is less sharp than those derived for the SML1 models. Thus, the presence of a bLVL is likely, even when the 2.4 Hz trough is ignored. The difference between generating a 2.4 Hz resonance or not will depend on the sharpness of the contrast between the SHVL and the bLVL.

Consequently, the geological interpretation of the SML2 models is similar to the interpretation of the SML1 models. The main difference lies in the properties of the bLVL such that a weaker bLVL will generate the 2.4 Hz resonance, due to the sharp contrast with respect to the SHVL, whereas a more consolidated bLVL will properly model the fHV curve without the 2.4 Hz trough.

Lower velocities are required immediately below the SHVL to model the fHV curve, with or without the 2.4 Hz trough. In this sense, Fig. 9(a) shows that models without a clear bLVL (lightest colours) will also generate an H/V peak around 8 Hz, but with a larger amplitude and a poor fitting of the left-flank decay. Thus, in order to fit both peak and decay, besides the SHVL, we require the presence of a bLVL, which becomes an essential feature to explain the whole fHV curve. Furthermore, as shown in Fig. 9(b), if the layer below the bLVL has higher relative velocities, then a secondary peak at lower frequencies, which is not observed from the data, would be generated. Therefore, general low velocities are required below the SHVL.

Figure 9.

Synthetic modelling of H/V using different velocity profiles below the SHVL. (a) Different overall velocities below the SHVL. (b) Fixed bLVL velocity and varying velocity of the half-space below.

The lower velocities below the SHVL can be explained by the presence of a sedimentary layer and also differences in the physical weathering process between the Early Hesperian and Early Amazonian lavas. This difference between the shallow Amazonian and deeper Hesperian layers is always required, either for the SML1 or SML2 models, as it allows for the characteristic steep left flank of the fHV curve.

Another relevant aspect is the different velocity ranges explored in the parameter spaces of the SML1 and SML2 models. On Earth, Lesage et al. (2018) studied the velocity profiles for different volcanic areas, including andesitic and basaltic deposits, and they found that shear-wave velocities lower than 500 m s⁻¹ are typically observed in the first 50 m depth. In a similar way, Panzera et al. (2015) compared a stratigraphic sequence from borehole data with the velocity model derived from H/V inversion for a site near Mt. Etna. They associated shear-wave velocities between 500–1000 m s⁻¹ with the corresponding basaltic layer underneath. These regions on Earth share some similarities with western Elysium Planitia (Golombek et al. 2017, 2020b). Under this context, the velocities obtained for models S1 and S2 are well within the expected range of shallow velocities in a volcanic region. Assuming these lower velocities, the thickness of the SHVL would be no greater than 15 m, which falls in the reasonable range of basaltic lava flow thicknesses. Nevertheless, shear-wave velocities around 2000 m s⁻¹ have been obtained from laboratory experiments for basaltic rock samples under the influence of confining pressure (e.g. Vinciguerra et al. 2005). Also, differences between the shear wave velocities derived from indirect seismic and lab samples have been revealed (Lesage et al. 2018). Therefore, an SHVL with higher velocities is still possible, which would be associated with a thicker layer.

5.1.3 Compatibility with compliance observations

Besides the H/V analysis, the ground compliance is an independent observation that can help to understand the first tens of metres underneath the InSight landing site. In this regard, Onodera (2022) computed the ground compliance using the SEIS and pressure data for a large set of convective vortices occurring in the nearby region. The inversion under a Bayesian approach of this data set allowed them to provide a most likely model consisting of a thin shallow regolith (∼1 m thick) followed by a coarse ejecta layer without further discontinuities below (see Figs 7, 8 and 10), at least down to ∼75 m.

Figure 10.

Characteristic models from set of models S1, M1 and L1 (top) and S2, M2 and L2 (bottom) and their synthetic ground compliance. (a) Shear wave velocity structure for characteristic models from the SML1 (top) and SML2 models (bottom). The model proposed by Onodera (2022) is plotted with dashed red line. The respective synthetic compliance response and the measured ground compliance (black curve and dark grey area) are plotted for different wind levels: (b) 0–6 m s⁻¹, (c) 6–12 m s⁻¹ and (d) 12–18 m s⁻¹.

Our proposed models differ from the models obtained from the inversion of compliance data as we are unable to properly resolve the first few metres of the subsurface structure at the InSight landing site, whereas the compliance inversion is highly sensitive to the shallow depths close to the surface. Nevertheless, our models are constrained by the velocities derived from the hammering experiment, which are in good agreement with those from the compliance inversion.

In order to assess the validity of our proposed models in terms of the observed ground compliance, among the best-fitting models, we chose one model representing each S1, M1, L1 and S2, M2, L2 model subsets, and computed their synthetic ground compliance response following the same approach as in Onodera (2022). The characteristic models and their synthetic compliance are plotted in Fig. 10.

The synthetic ground compliance modelling shows that our proposed models generate a compliance response that fits well within the error range of the measurements. Because of this, the rejection of models based on the compliance measurements can be hardly carried out. Particularly, the presence of the SHVL and the bLVL cannot be discarded. It is important to note that, given the Bayesian approach, the models obtained from the inversion of the compliance data are the most likely models given a set of a priori conditions. Nevertheless, models with unexpected rigidity such a high- or a low-velocity layer can also explain the observations, as shown in Fig. 10 and Fig. S2 (see Supporting Information).

5.2 The relevance of the marsquakes data

One of the advantages of using ambient vibrations data for the H/V analysis is certainly the abundance of recordings. On the contrary, using seismic events on Mars is less efficient not only due to the low seismic activity (Giardini et al. 2020) but also because the events are required to excite a wide frequency range.

Even though more than 900 events have been reported by the MQS, only the VF events have excited the resonance frequency at 8 Hz and just a few of them have been properly recorded by SEIS. The scarcity of these events is not a problem because we make use of their S-wave coda, which is produced by the multiple scattering of seismic waves. The multiple scattering phenomenon is expected to be the same for events occurring at different distances and/or azimuths, as long as they carry enough seismic energy above 5 Hz (Aki & Chouet 1975). Thus, future events exciting the high-frequency part of the spectrum are expected to generate the same H/V curve from the S-wave coda, which can be interpreted as the interference of multiple reverberations due to local conversions in the subsurface structure underneath, which generates surface-wave like ground motion (rather than recording direct Rayleigh or Love waves). Therefore, the full wavefield approach is an appropriate technique for modelling the corresponding H/V.

Although the hypocentral distance can be estimated by the S_g—P_g time, determining the backazimuth of the VF events has been not possible yet due to the depolarization of the first arrival (van Driel et al. 2021). Therefore, the epicentre locations of these events have not been obtained yet and the high-frequency propagation mechanism is still unclear (Clinton et al. 2021; van Driel et al. 2021). At least, according to the meteorological data, these events seem to be unrelated to atmospheric or environmental perturbations. A possible mechanism explaining the propagation of such high frequencies for very long distances (i.e. large S_g—P_g time) might be similar to the propagation of T-phases on Earth through the SOFAR channel (Walker et al. 1992; Okal 2008), but travelling through a thick waveguide with a really thin scattering layer inside so the energy is not leaked out. Further investigation of this mechanism is out of the scope of this work.

An interesting feature of the VF events is the almost negligible vertical ground-motion in the frequency range between 6 and 9 Hz. This observation can be explained by a strong local site effect at the InSight landing site due to the subsurface structure underneath. Given an input ground motion, the horizontal ground motion is largely enhanced while the vertical one is weakly augmented. In fact, the horizontal ground motion is about 7–8 times larger than the vertical one (see Figs 5 and 6) due to the site effect and so the amplitudes may easily surpass the horizontal ambient noise. On the other hand, the apparent absence of vertical ground motion can be understood as it having lower amplitudes than the ambient noise. This is also revealed by the analysis of the average ground acceleration around 8 Hz for all the events with quality A, B and C from the MQS catalogue, as shown in Fig. 11. In this case, the data shows that, given a horizontal ground motion slightly above the noise level of the quietest period of the Sol, the expected vertical ground-motion lies below the median vertical noise level obtained from the same quiet period, when most of the marsquakes are recorded. Hence, as the vertical noise level is higher than the expected vertical ground motion, there is no seismic signal on the vertical component associated with such event.

Figure 11.

Vertical and horizontal ground acceleration around 8 Hz for all the marsquakes with quality A, B and C reported by the MQS. The red dashed lines represent the median ground acceleration recorded during the quietest period of Sols 422 and 423 for both vertical and horizontal components. Different markers correspond to different type of marsquakes. The dotted–dashed blue line is a linear regression between the horizontal and the vertical acceleration, using events with SNR_Z > 2 dB (related to the H/V amplitude at 8 Hz). The filling pattern of the marker indicates SNR > 2 dB for the horizontal (lower half), vertical (upper half) or both components (full filled).

6 CONCLUSIONS

We used the VBB and SP seismic data recorded by SEIS instruments at the InSight landing site to investigate the H/V spectral ratios at frequencies between 0.4 and 10 Hz. The H/V was computed from both ambient vibrations (nHV) and seismic events data (eHV).

The H/V curves obtained from VBB and SP data are similar for frequencies below 10 Hz whereas, above this frequency, the VBB H/V curves (both nHV and eHV) show higher amplitudes. This difference is related to a different amplification of the vertical component at high frequencies. However, the origin of this difference is still unclear, so we suggest ignoring the frequencies above 10 Hz for any further analysis of the H/V curve.

Different nHV curves are observed during different periods of the martian day (Sol) due to the daily variation of the atmospheric conditions (mainly wind speed and temperature). These nHV curves are unreliable and do not represent the subsurface structure beneath InSight because the recorded seismic signals either (i) correspond to seismic signal close to the instrumental noise during the quietest period of the Sol, usually the martian evening when low-wind velocities are recorded (< 3 m s⁻¹), or (ii) correspond to non-natural lander mechanical noise, including the strong excitation of lander resonances, during the noisier periods of the martian Sol, usually the daytime when high-wind speeds are recorded. Unlike most places on Earth, using the ambient vibrations to reveal the shallow structure through the H/V technique is unfeasible at the InSight landing site on Mars.

The eHV curves were obtained from the S-wave coda of all recorded martian seismic events with energy well above the background noise. The resulting eHV curve shares similarities and differences with the nHV curves, particularly during the low-noise period. They both present a trough at 2.4 Hz and similar shape for frequencies below ∼4 Hz but a steep left flank leading to a peak at around 8 Hz is observed on the eHV curve. The inversion of the fHV (the eHV from VBB data) under the full diffuse field approach results in three different subsets of models, all of them suggesting the presence of a SHVL followed by a bLVL. Because the origin of the 2.4 Hz resonance remains unclear, the fHV curve was also inverted in the frequency ranges between 0.4–2 and 3–10 Hz. The resulting models share the presence of a high-velocity layer emplaced over layers with lower velocities and also a bLVL. In any case, the features of the H/V curve can be properly explained by the effect of the full wavefield and not only one type of waves.

We propose the emplacement of a high velocity layer at shallow depths (SHVL), likely an effusive Early Amazonian basaltic lava flow, to explain the strong excitation of the 8 Hz resonance during the VF events. We are unable to resolve the subsurface structure above the SHVL, but the pre-landing models and inversion of compliance data suggest the presence of a thin regolith followed by a coarse ejecta layer, which overlay the SHVL and are compatible with the fHV. Lower velocities below the SHVL are required to obtain the observed steep left flank of the 8 Hz peak. These deeper low velocities may be related to a more physically and chemically weathered Early Hesperian or Early Amazonian basaltic unit. A bLVL, likely corresponding to a sedimentary layer interrupting these two basaltic units, is needed to explain the fHV curve with or without the 2.4 Hz trough. The 2.4 Hz resonance is generated as long as there is a sharp contrast between the bLVL and the SHVL.

SUPPORTING INFORMATION

Supplementary data are available at GJI online.

Figure S1. Theoretical H/V curves for two exemplary models obtained from the SML1 inversion results. Left: H/V curves for an exemplary S1 model. Right: H/V curves for an exemplary L1 model. The H/V curves are computed for: the full wavefield (blue), body waves (orange), surface waves (green), P waves (red), S waves (magenta) and the fundamental mode of the Rayleigh waves (brown).

Figure S2. Discrete subsurface structure models including an SHVL (black lines) picked from the Bayesian inversion performed by Onodera (2022). The density maps correspond to the output probability density distributions of the Young modulus from surface to 200 m depth, where red and blue colours show high and low probabilities, respectively.

Figure S3. Spectral density plot of H/V curves of the S-wave coda of marsquakes data recorded by VBB (left-hand panel) and SP (right-hand panel). Lighter colours indicate a higher spectral density.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

The authors acknowledge National Aeronautics and Space Administration (NASA), Centre National D’études Spatiales (CNES), their partner agencies and institutions (United Kingdom Space Agency [UKSA], Swiss Space Office [SSO], Deutsches Zentrum für Luft-und Raumfahrt [DLR], Jet Propulsion Laboratory [JPL], Institut du Physique du Globe de Paris [IPGP]—Centre National de la Recherche Scientifique-École Normale Supérieure [CNRS], Eldgenössische Technische Hochschule Zürich [ETHZ], Imperial college [IC], Max Planck Institute for Solar System Research [MPS-MPG]), and the flight operations team at JPL, SEIS on Mars Operation Center (SISMOC), Mars SEIS Data Service (MSDS), Incorporated Research Institutions for Seismology–Data Management Center (IRIS-DMC) and Planetary Data System (PDS) for providing SEED Seismic Experiment for Interior Structure (SEIS) data. We acknowledge funding from (1) Swiss State Secretariat for Education, Research and Innovation (SEFRI project ‘Marsquake Service-Preparatory Phase’), (2) ETH Research grant ETH-0617–02, and (3) ETH + 02 19–1: Planet MARS. The research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). This is the InSight contribution number 268.

DATA AVAILABILITY

The Interior Exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) seismic event catalogue version 9 (InSight Marsquake Service 2022), the waveform data and station metadata are available from the Institut du Physique du Globe de Paris (IPGP) Datacenter and Incorporated Research Institutions for Seismology Data Management Center (IRIS-DMC, InSight Mars SEIS Data Service 2019b). Seismic waveforms are also available from the National Aeronautics and Space Administration Planetary Data System (NASA PDS), available at https://pds.nasa.gov/(last accessed April 2022, InSight Mars SEIS Data Service 2019a). The channel location and codes follow an adapted version of the Standard for the Exchange of Earthquake Data (SEED) naming convention. The seismic catalogue (InSight Marsquake Service 2022) provides information on all detected events, including the event classification based on frequency content, phase picks, backazimuth estimates and event quality estimates. The code for the computation of the theoretical H/V under the diffuse field approach is available upon request from L. Margerin. Data processing and plotting has been done with Obspy (Beyreuther et al. 2010), Numpy (Harris et al. 2020) and Matplotlib (Hunter 2007). The authors declare no conflicts of interest to disclose.