Rayleigh Wave Ellipticity Modeling and Inversion for Shallow Structure at the Proposed InSight Landing Site in Elysium Planitia, Mars

The SEIS (Seismic Experiment for Interior Structure) instrument onboard the InSight mission will be the ﬁrst seismometer directly deployed on the surface of Mars. From studies on the Earth and the Moon, it is well known that site ampliﬁcation in low-velocity sediments on top of more competent rocks has a strong inﬂuence on seismic signals, but can also be used to constrain the subsurface structure. Here we simulate ambient vibration waveﬁelds in a model of the shallow sub-surface at the InSight landing site in Elysium Planitia and demonstrate how the high-frequency Rayleigh wave ellipticity can be extracted from these data and inverted for shallow structure. We ﬁnd that, depending on model parameters, higher mode ellipticity information can be extracted from single-station data, which signif-icantly reduces uncertainties in inversion. Though the data are most sensitive to properties of the upper-most layer and show a strong trade-off between layer depth and velocity, it is possible to estimate the velocity and thickness of the sub-regolith layer by using reasonable constraints on regolith properties. Model parameters are best constrained if either higher mode data can be used or additional constraints on regolith properties from seismic analysis of the hammer strokes of InSight’s heat ﬂow probe HP 3 are available. In addition, the Rayleigh wave ellipticity can distinguish between models with a constant regolith velocity and models with a velocity increase in the regolith, information which is difﬁcult to obtain otherwise.


Introduction
Propagation through a soft soil layer can significantly amplify ground motion amplitudes, specifically on the horizontal components, resulting in strong site effects which may considerable increase earthquake damage (e.g. Borchert 1970;Anderson et al. 1986;Sánchez-Sesma and Crouse 2015). Selective amplification of the horizontal components' amplitudes has also been observed in data recorded on the Moon by the Apollo seismic lunar network (Lammlein et al. 1974;Nakamura et al. 1975) and attributed to resonances in the surficial layer of lunar regolith. A similar effect can be expected from the regolith at the proposed landing site for NASA's InSight mission. This mission will for the first time place a three-component broad-band seismometer and colocated three-component short-period seismometer, the SEIS (Seismic Experiment for Internal Structure) instrument, on the surface of Mars in Elysium Planitia in November 2018 (Banerdt et al. 2013). A study of the expected site amplification is not only important to understand how and in which frequency range it will affect the recorded seismograms, but also because the observed amplification can help to constrain the elastic properties of the regolith at the landing site. As all seismic waves recorded by SEIS pass through the surficial regolith layer, understanding its properties will reduce the uncertainty associated with other seismic measurements. In addition, the elastic properties of the Martian soil are of profound interest for future robotic and human exploration missions.
The ambient vibration horizontal-to-vertical spectral amplitude ratio (H/V) is a common tool for estimating site effects and soil properties with a single station (Nakamura 1989). The resulting H/V curve often shows a prominent frequency peak that provides a good proxy for the fundamental resonance frequency of the site (e.g. Lachet and Bard 1994;Lermo and Chávez-García 1994;Malischewsky and Scherbaum 2004;Bonnefoy-Claudet et al. 2008). Thus, microzonation in densely populated, earthquake-prone regions often makes use of H/V measurements to map variations in resonance frequencies (e.g. Panou et al. 2005;Bragato et al. 2007;Souriau et al. 2007;Bonnefoy-Claudet et al. 2009;Picozzi et al. 2009;Poggi et al. 2012). Inverting the H/V curve for shallow sub-surface structure requires some understanding of which part of the noise wavefield is responsible for this effect, though, and different theories have been put forward to that end: Nakamura (2000Nakamura ( , 2008 explains the H/V peak by SH-wave resonances in the soft surface layer, whereas a number of other authors consider the H/V curves as measurements of the frequency-dependent Rayleigh wave ellipticity (e.g. Lachet and Bard 1994;Lermo and Chávez-García 1994;Konno and Ohmachi 1998;Fäh et al. 2001;Bonnefoy-Claudet et al. 2006). Bonnefoy-Claudet et al. (2008) show that, for a variety of structural models, the H/V peak frequency provides a good estimate of the theoretical 1D resonance frequency, regardless of the contribution of different wave types to the wavefield. In these simulations, surface waves are found to dominate the wavefield for high to moderate impedance contrasts between sediment and bedrock and surficial sources (Bonnefoy-Claudet et al. 2006. However, both simulations (Bonnefoy-Claudet et al. 2008) and actual measurements (Okada 2003;Köhler et al. 2007;Endrun et al. , 2011Poggi et al. 2012) indicate that Love waves may contribute significantly to the measured H/V curves, and their contribution may also vary with frequency and time. Accordingly, recent studies focus on either extracting the Rayleigh wave ellipticity from the ambient vibration wavefield before inversion, or modeling of the complete noise wavefield using either diffuse field theory García-Jerez et al. 2013;Kawase et al. 2015;Lontsi et al. 2015) or stochastic fields (Lunedei and Albarello 2015). As discussed in the overview by Hobiger et al. (2012), Rayleigh wave ellipticity can be estimated from ambient noise recordings by both array and single station methods. Single station methods are either based on time-frequency analysis using a continuous wavelet transform (Fäh et al. 2001(Fäh et al. , 2009Poggi et al. 2012), or on the random decrement technique (Hobiger et al. 2009(Hobiger et al. , 2013Bard et al. 2010;Garofalo et al. 2016;Gouveia et al. 2016). Hobiger et al. (2013) investigate which part of the Rayleigh wave ellipticity curve contains relevant information on soil structure and is thus most useful in an inversion. They find that for curves with a strong singularity, the right flank of the ellipticity peak together with the peak frequency, which might be constrained by including the left flank, is the most informative part. This is consistent with observations by Fäh et al. (2001) that H/V curves are most stable and dominated by Rayleigh wave ellipticity in the frequency band between the fundamental resonance peak and the first minimum, where they are determined by the layering of the sediments. However, Scherbaum et al. (2003) have shown that the inversion of Rayleigh wave ellipticity alone is subject to strong trade-offs between layer thickness and average layer velocity. Better results have been obtained when combining ellipticity curves with other information, e.g. stratigraphic layering (Mundepi et al. 2015) or sediment velocities (Yamanaka et al. 1994;Satoh et al. 2001;Arai and Tokimatsu 2008) from borehole logging, or surface wave dispersion measured actively or passively with surface arrays Tokimatsu 2005, 2008;García-Jerez et al. 2007;Picozzi and Albarello 2007;Hobiger et al. 2013;Dal Moro 2015). In the later case, the inclusion of ellipticity information can significantly improve estimates of bedrock depth (Garofalo et al. 2016;Gouveia et al. 2016) and velocity . Besides, Lontsi et al. (2015) recently found that the inversion trade-offs can be resolved through the additional use of H/V curves measured at one or several borehole receivers at depth.
H/V curves for the Apollo lunar seismic data have been determined from the coda waves of shallow source events, due to the lack of a continuous ambient background wavefield strong enough to be observable by the Apollo seismometers (Lognonné et al. 2009), and interpreted in terms of Rayleigh wave ellipticity by Mark and Sutton (1975), Nakamura et al. (1975) and Horvath et al. (1980). Regolith thickness has then been obtained by using Pwave velocity results of the active seismic experiments as prior information. More recently, Dal Moro (2015) inverted H/V curves for two Apollo sites in combination with dispersion measurements from the co-located lunar active seismic experiments, considering both Love and Rayleigh wave contributions. On Mars, InSight is expected to observe a micro-seismic background wavefield caused by atmospheric sources, mainly in the form of surface waves as the sources interact with the surface of the planet. This background wavefield can thus be used to extract the high-frequency Rayleigh wave ellipticity and invert it for shallow subsurface structure.
Here, we simulate this application of Rayleigh wave ellipticity measurements and inversion of the ellipticity curve to InSight SEIS data from Mars. Based on a priori information on the landing site geology and laboratory measurements of seismic velocities in regolith analogue material, we build a plausible model reaching to ∼ 50 m depth and generate a noise wavefield that consists of fundamental and higher mode Rayleigh and Love waves. We show how the Rayleigh wave ellipticity can be extracted from this synthetic dataset and inverted for ground structure using the Conditional Neighbourhood Algorithm (Sambridge 1999;Wathelet 2008). We discuss how reasonable variations in the elastic parameters may alter the obtained ellipticity curve, how site amplification will influence the recorded seismograms, and how a combination of ellipticity information with other data can best be used to constrain properties of the shallow subsurface at the InSight landing site.

Construction of the Model
The geology and shallow subsurface structure of the InSight landing site was determined by mapping and analyses described in more detail by . The landing site in western Elysium Planitia is located on a surface mapped as Early Hesperian transition unit (Tanaka et al. 2014) and is most likely volcanic based on: 1) the presence of rocks in the ejecta of fresh craters ∼ 0.4-20 km diameter arguing for a strong competent layer ∼ 4-200 m deep and weak material above and beneath (e.g. , 2) exposures of strong, jointed bedrock overlain by ∼ 10 m of relatively fine grained regolith in nearby Hephaestus Fossae , 3) platy and smooth lava flows mapped in 6 m/pixel CTX images south of the landing site (V. Ansan Mangold, written comm.), and 4) the presence of wrinkle ridges, which have been interpreted to be fault-propagation folds, in which slip on thrust faults at depth are accommodated by asymmetric folding in strong, but weakly bonded layered material (i.e., basalt flows) near the surface (e.g. Mueller and Golombek 2004). The thermophysical properties of the landing site indicates the surface materials are composed of cohesionless, very fine to medium sand (particle sizes of 40 μm to 400 μm, average particle size ∼ 170 μm) or very low cohesion (< few kPa) soils to a depth of at least several tens of centimeters with surficial dust less than 1-2 mm. Highresolution images of the landing site and surrounding areas show surface terrains that are dominantly formed by impact and eolian processes . The sand grains are likely equant to rounded by saltation as they are exposed to surface winds by repeated impacts (e.g. McGlynn et al. 2011).
The landing ellipse, sized 130 km by 27 km, is located on smooth, flat terrain that generally has very low rock abundance . Most rocks at the landing site are concentrated around rocky ejecta craters larger than 30-200 m diameter, but not around similarly fresh smaller craters Warner et al. 2016a). Because ejecta is sourced from shallow depths, ∼ 0.1 times the diameter of the crater (Melosh 1989), the onset diameter of rocky ejecta craters has been used to map the thickness of the broken up regolith. Results indicate a regolith that is 2.4-17 m thick at the landing site (Warner et al. 2016a,b), that grades into large blocky ejecta over strong intact basalts. This is also consistent with regolith thickness estimates based on morphometric properties of concentric craters (Warner et al. 2016a) and SHARAD radar analysis suggesting low-density surface material overlying more intact rock within 10-20 depth of the surface . Because craters larger than 2 km do not have rocky ejecta, material below the basalts at ∼ 200 m depth is likely weakly bonded sediments. An exposed escarpment of nearby Hephaestus Fossae (Fig. 1) shows this near surface structure with ∼ 5 m thick, fine grained regolith, that grades into coarse, blocky ejecta with meter to ten-meter scale boulders that overlies strong, jointed bedrock. The grading of finer grained regolith into coarser, blocky ejecta is exactly what would be expected for a surface impacted by craters with a steeply dipping negative power law distribution in which smaller impacts vastly outnumber larger impacts that would excavate more deeply beneath the surface (e.g. Shoemaker and Morris 1969;Hartmann et al. 2001;Wilcox et al. 2005). Fragmentation theory in which the particle size distribution is described by a negative binomial function (Charalambous 2014) was applied to the InSight landing site using rock abundance and cratering size-frequency measurements to derive a synthesized regolith with a relatively small component of particles > 10 cm (Charalambous et al. 2011;Charalambous and Pike 2014;).

Fig. 1
Shallow structure nearby the InSight landing site. HiRISE image PSP_002359_2020 of a portion of the Hephaestus Fossae in southern Utopia Planitia at 21.9 • N, 122.0 • E showing ∼ 4-10 m thick, fine grained regolith, that grades into coarse, blocky ejecta that overlies strong, jointed bedrock. Image shows a steep escarpment with talus on the steep slope below As a result, for our modeling, we use a baseline model with an intermediate regolith thickness of 10 m (Fig. 2a). We discuss the influence that a different regolith thickness within the estimated range will have on the results in Sect. 3.3.
We proceeded as follows to translate this subsurface structure into a seismic velocity model ( Fig. 2b and Table 1): The regolith velocity is based on laboratory experiments with three regolith simulants, for which compressional (v P ) and shear (v S ) wave velocities were determined under various confining pressures corresponding to lithostatic stresses at 0-30 m depth . For all regolith simulants, a power-law increase of velocities with depth was observed , as is also common for terrestrial soils (e.g. Faust 1951;Prasad et al. 2004). The velocities used here are based on the results for two sands (Mojave sand and Eifelsand), which are rather similar, as these sands are closer in particle size to the expected regolith in Elysium Planitia than the third, rather fine-grained, silty simulant tested. The low velocities and low v P /v S ratios obtained also agree with laboratory data on dry quartz sands (Prasad et al. 2004) as well as terrestrial in situ measurements on shallow unconsolidated sands (e.g. Bachrach et al. 1998). A velocity increase from the surface to the value corresponding to the maximum depth of the regolith layer was implemented in the model (Fig. 2), spanning 20 layers, and the unit mass density as used in the lab tests assumed.
Below the sandy regolith, somewhat more blocky ejecta are expected, based on less frequent larger impact that would eject material from deeper levels. Velocities in this layer are based on field measurements in an analogue environment on Earth, lava flows in the Californian Mojave Desert (Wells et al. 1985). The stratigraphy of the Cima volcanic field consists of a thin layer of tephra and eolian material on top of a so-called rubble zone of basaltic clasts, grading into highly fractured basaltic flow rock. P-wave velocities of the dif-    Table 1 ferent units have been determined along seismic profiles. The shallowest layer in this area consists of silt, so it cannot be compared to the sandy regolith at the InSight landing site. However, the rubble zone beneath is considered equivalent to the coarse ejecta which have a similar thickness as the sandy regolith , whereas the upper-most part of the basalt flows in Elysium Planitia is also expected to exhibit some crack damage that will reduce the seismic velocities compared to pristine basalt (Vinciguerra et al. 2005). The velocity model tries to mimic these variations and includes gradational changes between the different layers. Based on a HiRISE image of a steep exposed portion of Hephaestus Fossae, southern Utopia Planitia , gradient layers have a thickness of 1 m between the regolith and the coarse ejecta and between the coarse ejecta and the fractured basalt, whereas the change from fractured to unfractured basalt extends over a larger depth range. S-wave velocities in these layers were derived from the P-wave velocities by assuming v P /v S decreasing from 1.9 to 1.8 with depth.
In addition to (mainly) shear wave velocities, the Q factor has a non-negligible influence on H/V curves: By modeling, Lunedei and Albarello (2009) showed that damping has a significant effect on H/V peak amplitudes, and concluded that Q values, which are otherwise difficult to obtain, might be derived from H/V curves. For the Moon, high Q values in the upper few 100 m of the lunar subsurface strongly influence the H/V peak amplitude and are essential in obtaining a good fit to the measured data (Dal Moro 2015). The unusually high Q values observed on the Moon (Nakamura and Koyama 1982) are caused by the extremely dry rocks from which even thin layers of adsorbed water have been removed by strong outgassing under vacuum conditions (Tittmann 1977;Tittmann et al. 1979). In the Martian crust, a comparable evacuation of trapped fluids is prevented by atmospheric pressure. Accordingly, Q is predicted to be larger by at most a factor of two compared to Earth (Lognonné and Mosser 1993). The above studies only consider rocks at larger depths, though, and not the properties of surficial soils. Any liquid or frozen surface water would not be in equilibrium in the equatorial regions of Mars and thus quickly sublimate, and during planetary protection review, it was confirmed that the InSight landing site does not contain any water or ice within 5 m of the surface, nor high concentrations of water bearing minerals . Evidence for the water content within the Martian regolith, though, is provided by neutron measurements by Mars Odyssey, which give a lower limit of 3-6 wt% water abundance in the upper metre of Martian regolith near the InSight landing site (Feldman et al. 2004), and analysis of Mars Express infrared reflectance spectra, which finds similar values for the upper surface layer of the regolith in this region (Milliken et al. 2007). The water could be present either in the form of hydrous minerals, which would be stable under Martian P-T conditions (Bish et al. 2003), or adsorped water (Möhlmann 2008). Laboratory measurements on crushed volcanic ash, although with a smaller particle size than expected for the regolith at the InSight landing site, indicate a liquid-like water content of at least two monolayers down to −70 • C (Lorek and Wagner 2013). This "sorption water" is supposed to reside mainly below depths of a few decimetres, outside the range of Martian diurnal and seasonal thermal cycles (Möhlmann 2004). Laboratory measurements have shown that already a few monolayers of adsorbed water can drastically reduce the high Q values observed in outgassed lunar or terrestrial samples (Tittmann et al. 1979). Thus, we assume that Q values of the Martian regolith and shallow subsurface are within the range of one to two times the terrestrial values. Terrestrial values are estimated by using the rule of thumb Q S = v s /10 (Dal Moro 2014Moro , 2015 for the regolith, which is consistent with the low Q S values obtained by borehole measurements in terrestrial sediments (e.g. Parolai et al. 2010;Fukushima et al. 2016), and taking Q S = 400 for the basalt. In the model, we set Q S to 1.5 times these values, resulting in values between 20 and 30 for the regolith. Based on laboratory measurements on dry quartz sands (Prasad et al. 2004), Q P is set to equal Q S for the regolith, and increased to 1.5 and 2 times Q S for the coarse ejecta and the basalt, respectively.

Synthetic Seismograms
Synthetic seismograms simulating the ambient vibration wavefield are calculated by using a modal summation technique (Herrmann 2013) for a multitude of surface sources (e.g.  and the 1-D model developed above (Fig. 2, Table 1). Five thousand sources are randomly distributed at distances of up to 5 km from the station and randomly activated up to 5 times, with randomly varying amplitudes. The source signals are delta-peak force functions. In total, 15,195 such delta forces were applied during a recording time of 30 min for the synthetics. We created two different data sets, using either only fundamental modes or also higher modes of surface waves in the source process for randomly inclined forces, generating both Rayleigh and Love waves. In this way, waves from different directions and with a different amount of Love and Rayleigh waves may reach the recording station at the same time and interfere with each other. In case of the multi-mode wavefield, the relative content of fundamental mode and higher mode energy arriving at the recording station is also variable in time depending on the orientation and distance of the active forces.
The synthetics thus created present a simplification in that no body waves, including those caused by scattering, are considered. Based on examples from other planets, we can assume a significant presence of surface waves in the wavefield caused by surface sources on Mars. Rayleigh waves in the frequency range considered here are routinely analysed in ambient vibrations array recordings of Earth data when studying site effects (e.g. Satoh et al. 2001;Kind et al. 2005;Picozzi and Albarello 2007;Hannemann et al. 2014;Garofalo et al. 2016) and have also been extracted from ambient vibrations in the highly scattering environment of the Moon (Larose et al. 2005;Tanimoto et al. 2008;Sens-Schönfelder and Larose 2010). The two methods presented in the following to extract Rayleigh wave ellipticity from ambient vibrations have both been successfully tested on synthetic wavefields that contain both body and surface waves (Fäh et al. 2009;Hobiger et al. 2009Hobiger et al. , 2012 and applied to Earth data (e.g. Poggi et al. 2012;Gouveia et al. 2016). We thus demonstrate our inversion approach using synthetics that contain surface waves only in a first-order approximation of the actual ambient vibration wavefield.

Extraction of Rayleigh Wave Ellipticity
The standard H/V ratio is calculated by using the squared average of the horizontal signal components. However, if the wavefield contains Love or SH waves, they will be present on the horizontal components only and lead to an overestimation of H/V amplitudes. Accordingly, other methods are needed to directly estimate the ellipticity from the signals. We compare two different methods to extract Rayleigh wave ellipticity from single station recordings. Both make use of the phase shift of π/2 between vertical and horizontal components of particle motion that is characteristic of Rayleigh waves.
The first method, called HVTFA (H/V using time frequency analysis, Fäh et al. 2009) and originally proposed by Kristekova (2006), uses a continuous wavelet transform based on modified Morlet wavelets (Lardies and Gouttebroze 2002) to transform the three signal component into the time-frequency domain. Rayleigh waves are identified by scanning for maxima in the transformed vertical component in each frequency band (Kristekova 2006). Love or SH waves that contain horizontal energy only are thus effectively excluded from further consideration. For each maximum on the vertical component, the corresponding maximum value on the horizontal components with a phase shift of ±π/2 is identified and used to calculate an ellipticity value. All values derived for a given frequency are analysed statistically via filtering of histograms (Fäh et al. 2009). HVTFA is implemented as a module in the GEOPSY software (www.geopsy.org) and requires two input parameters, the Morlet wavelet parameter that controls the wavelet's width in the spectral domain and the number of maxima on the vertical component selected per minute. Based on the study reported by Fäh et al. (2009), we selected a value of 8 for the Morlet wavelet parameter and choose 5 maxima per minute. The above study found that in general, the number of selected maxima per minute should be in the range 1-5 or less, with preference to lower values. We checked that the extracted average curves were comparable for 2 and 5 maxima per minute, and chose the larger number due to the short time window analysed here, compared to two hours in the above study, to get meaningful statistics. Besides, a larger number of maxima allows for a better identification of a higher mode ellipticity curve (see below). The method has previously been demonstrated on synthetic wavefields containing body and multi-mode surface waves (Fäh et al. 2009;Hobiger et al. 2012) and applied to measured data (Poggi et al. 2012) at frequencies up to at least 15 Hz.
The second method is RAYDEC (Hobiger et al. 2009), based on the random decrement technique (Cole 1973). For this technique, the signals are split into short analysis time windows based on the number of zero crossings of the vertical component seismogram within narrow frequency bands. These short time windows are shifted for the horizontal components to accommodate the π/2 phase shift characteristic of Rayleigh waves. Then, an optimum rotation angle for the radial direction of the signals is determined by maximizing the correlation between the rotated horizontal components and the vertical component. Horizontal and vertical components for all time windows are summed, using the correlations as weighting factors, and the ellipticity is obtained by dividing these sums. The weighting factors assure that time windows that do not predominantly contain Rayleigh waves are efficiently down-weighted in the ellipticity calculation. As pointed out by Hobiger et al. (2009), higher mode Rayleigh waves cannot be distinguished from the fundamental mode by this approach, though. The two free parameters in this method are the sharpness of the frequency bands used in filtering the data and the length of the short time windows used for the analysis. We follow the suggestions by Hobiger et al. (2009) in using a time window length of 10/f , but use a somewhat smaller bandwidth of 0.1f , where f is the central frequency of the respective filter band. The method has previously been demonstrated on synthetic wavefields containing body and multi-mode surface waves (Hobiger et al. 2009(Hobiger et al. , 2012 and applied to measured data (Hobiger et al. 2009(Hobiger et al. , 2013Garofalo et al. 2016;Gouveia et al. 2016) at frequencies up to 30 Hz.

Inversion
Inversion of Rayleigh wave ellipticity for shallow subsurface structure is a non-unique problem with a strong trade-off between layer thicknesses and velocities (Scherbaum et al. 2003). Accordingly, this non-uniqueness has to be explored during the inversion to provide a meaningful set of models that can explain the data within their uncertainties while at the same time allowing an estimate of the uncertainty in the model. Here, we use the Conditional Neighbourhood Algorithm implemented in GEOPSY ). The Neighbourhood Algorithm (NA), as introduced by Sambridge (1999), is a direct search algorithm based on Voronoi cells that preferentially samples the regions of parameter space showing a low misfit in a self-adaptive manner. It has the ability to escape local minima and can locate several disparate regions of low misfit simultaneously, while requiring a lower number of tuning parameters than comparable algorithms. The NA has been applied to a diverse range of geophysical inversion problems, including earthquake location (Sambridge and Kennett 2001;Oye and Roth 2003), inversion of receiver functions (Frederiksen et al. 2003;Sherrington et al. 2004), inversion of surface wave dispersion curves Erduran et al. 2008;Yao et al. 2008) and surface wave waveforms (Yoshizawa and Kennett 2002), and inversion of interferometric synthetic aperture radar data (Pritchard and Simons 2004;Fukushima et al. 2005). The Conditional NA adds the possibility to define irregular limits to the searchable parameter space based on physical conditions (e.g. constraints on v P /v S ratio, in addition to independent constraints on v P and v S ), numerical issues, or prior information (Wathelet 2008). Besides, a dynamic scaling is implemented to keep the exploration of the parameter space as constant as possible while the inversion progresses. The Conditional NA has found many applications in site characterization (e.g. Coccia et al. 2010;Kühn et al. 2011;Souriau et al. 2011;Di Giulio et al. 2012Hobiger et al. 2013;Michel et al. 2014;Mundepi et al. 2015;Gouveia et al. 2016).
The choice of model space parameterization (e.g. number of layers, range of velocities, depths and Poisson's ratios, velocity-depth laws) in the inversion process also influences results. In the case of surface wave dispersion curve inversion for shallow subsurface structure, this issue has for example been addressed by  and Di Giulio et al. (2012). As ellipticity curves are rarely inverted on their own, a comparable study focussing only on them is missing. In case of the InSight landing site, some prior knowledge on stratigraphy and layer thicknesses is available from analyses of orbital photography, as discussed above, as well as some constraints on regolith velocities from laboratory measurements. In the following, we investigate further how model parameterization, inclusion of prior information, and the parts of the ellipticity curve used in the inversion influence the results. We follow Di Giulio et al. (2012) in that we evaluate the different models based on the corrected Akaike's information criterion (AICc) and thus combine data misfit and model complexity (i.e. number of degrees of freedom) to rank the models. The AIC is an information-theoretical approach based on the idea to combine the Kullback-Leibler information number, indicating the loss of information when an approximating model is used to explain reality, and the maximum likelihood function (Kullback and Leibler 1952;Akaike 1973). It is expressed by where K is the number of free parameters. The first term in (1) is a measure of the misfit between the approximating model and the true representation of reality, and the second term penalizes model complexity, i.e. a large number of degrees of freedom. In contrast to simply using misfit to rank a model, the AIC also considers the trade-off between bias and variance in model selection, where a larger number of free parameters in the model will reduce the bias (or misfit) at the expense of increasing variance and leading to over-fitting. Models with a lower value of AIC are considered to be better models. A corrected form of the AIC (1) has been proposed for cases of least-square estimation with normally distributed errors and small sample sizes (Sugiura 1978;Hurvich and Tsai 1989): where n f is the number of observations, i.e. in our case the number of samples in the ellipticity curves,ê 2 is the sum of estimated residuals for candidate models divided by n f , which in our case is equivalent to the misfit between observed and modeled data, and K the number of free parameters, i.e. the degree of freedom of the model parameterization. Following Di Giulio et al. (2012), we use the AICc to rank the models resulting from our inversions of the measured ellipticity curves. The model sets we use consist of an increasing number of layers, from one to six, over a half-space. Within the first layer, velocities can be uniform or follow either a linear or a power-law velocity-depth function (Fig. 3). In the latter two cases, the upper-most layer is composed of five sublayers to allow for the velocity increase with depth according to the respective law. This type of parameterization is often used to mirror the increasing compaction of the subsoil with depth in sedimentary environments (Faust 1951). All layers below the first one have uniform, homogeneous velocities. Note that, to portray a realistic situation, the model used to compute the input waveforms is actually more complex than the parameterizations allow for in the inversion: It consists of three layers, but both the first and third , v bot P 0 , and v bot S0 have to be inverted for, which provides two additional degrees of freedom compared to the previous case. (c) Power-law velocity increase in the topmost layer, which also corresponds to four unknown velocity parameters in the topmost layer layer contain a velocity increase with depth, and the transition between layers is gradational. However, both the gradational transitions and the velocity structure of the lower-most layer have a very minor influence on the shape of the ellipticity curve (Figs. 4a and b) and thus cannot be resolved in the inversion, whereas the velocity structure within the first layer has a stronger influence on the shape of the ellipticity curve, especially the first higher mode (Fig. 4c), which might be resolvable. As we do not want to introduce additional, unconstrained complexity in the inversion that might not be warranted, we try to approximate the data by the most simple model, e.g. with first-order discontinuities between different layers, and investigate how well we can constrain its parameters.
For each layer, P-and S-wave velocity as well as layer thickness are allowed to vary. This leads to three degrees of freedom per layer for uniform layers, and five degrees of freedom in the upper-most layer if the velocity follows a power-law or linear dependence with depth, as velocities both at the top and bottom of this layer are free parameters in these cases (Fig. 3). In addition, two degrees of freedom are associated with the basal layer of the model, the bedrock (P-and S-wave velocity). In summary, the models possess a minimum of 5 and a maximum of 20 degrees of freedom in case of uniform layers and between 7 and 22 degrees of freedom if velocities of the first layer follow a power-law or linear dependence. For each parameterization, we consider 5 runs of the Conditional NA that start with different random seeds to assure robust results. For each run, we use n S0 = 250 starting models and N = 5000 iterations, where a new sample is added to each of the n S = 100 cells with the lowest misfit in each iteration step. In total, we thus sample n S0 + N × n S = 500,250 models in each run.

Ellipticity Measurement
In Fig. 5, we compare the results of the standard H/V processing as well as the HVTFA and RAYDEC methods to extract Rayleigh wave ellipticity for our baseline model in the case of a Love and Rayleigh wave field either containing fundamental modes only or including higher modes. All curves show a peak frequency f 0 in agreement with the theoretical prediction of 4.9 Hz. The observed peak frequency is also close to the fundamental peak of the SH resonances for the model, as typically observed for models with a strong impedance contrast (Bonnefoy-Claudet et al. 2008). f 0 also closely agrees with the rule of thumb λ/4 estimate of when using a model thickness z to the top of the gradient between regolith and coarse ejecta of 9.5 m and the travel time-based average S-wave velocity of the regolith v S obtained as . This results in a value of 184 m/s for v S and a resonance frequency equal to 4.84 Hz. A more general formula for multiple layers over a halfspace that also considers the position of individual layers has recently been published by Tuan et al. (2016). Application of this formula to our model provides an estimate of 4.94 Hz for the resonance frequency. The close agreement between these calculated values and the observed peak frequency indicates that the ellipticity peak is caused by the contrast between the regolith layer and the coarse ejecta layer of the model, whereas the lower layers do not seem to have a distinct influence on the curve. The influence of Love waves and the resulting discrepancy between the theoretical Rayleigh wave ellipticity for the model and the measured H/V curves is greatest at frequencies below and up to the peak frequency (Fig. 5a). The influence becomes stronger in the more realistic case that also contains higher modes (Fig. 5d). Both methods for ellipticity extraction lead to a closer fit to the theoretical curve, especially at frequencies above 7 Hz on the right flank of the peak and along the entire left flank of the peak. HVTFA seems to perform slightly better there than RAYDEC (Figs. 5b and c). For the case of the mixed mode wavefield, additional complexity arises from the prominent excitation of the first higher mode that influences the H/V curve through interference near the peak and the addition of a broad secondary peak around 12 Hz (Fig. 5d). Again, both HVTFA and RAY-DEC give a closer representation of the theoretical flanks of the fundamental mode curve, with HVTFA being somewhat better in the estimation of the lower flank (Figs. 5e and f). Some influence of the superposition of peaks near 5 Hz and the additional higher mode peak around 12 Hz remains visible in the RAYDEC result. As noted by Hobiger et al. (2009), this method cannot distinguish between different modes. In contrast, the HVTFA results permit the extraction of ellipticity curves for both fundamental and first higher mode in this case (Fig. 6). Due to the representation of results for individual time windows in histogram shape before averaging, minor but distinct contributions to the measurements at a specific frequency can be identified, i.e. higher modes. If a smaller number of peaks per minute is selected, the higher mode curve is increasingly suppressed in favour of the fundamental one. However, as stated above, we compared the fundamental mode curves for selecting both the 2 or 5 largest peaks per minute and got very consistent results. This makes us confident that selecting 5 peaks per minute to also capture the higher mode does not lead to degradation of the curves. This is also confirmed when comparing to the theoretical predictions (Fig. 5e).
In general, the H/V spectrum does not provide information about different modes. Attempts have been made to model the whole spectrum, assuming a mixture of Love and Rayleigh waves and of different modes (Arai and Tokimatsu 2004;Dal Moro 2015). However, this requires a priori assumptions, e.g. about source types and distribution, and the relative contribution of Love and Rayleigh waves to the noise wavefield.   have successfully extracted higher-mode ellipticity curves from threecomponent array recordings using high-resolution frequency-wavenumber analysis, which, however, requires recordings at a set of at least 10 stations. To our knowledge, the extraction of higher mode ellipticity information from single station recordings has not been demonstrated so far, neither using synthetic nor actual measured data. When comparing ellipticity results from array analysis to single station HVTFA results, Poggi et al. (2012) in fact state that a disadvantage of the latter method is that it is not capable of separating contributions from several different modes. For the configuration considered here, a clear identification of branches belonging to several modes is possible in the HVTFA results, though (Fig. 6).
We proceed by using the HVTFA results as input for ellipticity curve inversion. They are closer to the theoretical curves than the RAYDEC results and, in contrast to them, also allow to study the effect of including higher mode information.

Inversion
Inversion of Rayleigh wave ellipticity is a nonlinear problem, which is also non-unique with a trade-off between layer thickness and velocity. Under these circumstances, model parameterization can significantly influence inversion results. We use the AICc as a way to combine the misfit between measured and modeled data and the model complexity, given by the number of degrees of freedom in the model, in a single number for model ranking.

Unconstrained Parameter Space
In a first step, we allow for a broad range of parameter values and do not impose a priori constraints, e.g. on the regolith thickness or the bedrock velocities, to investigate how the non-uniqueness of the problem is captured by the inversion and how the inclusion of different parts of the ellipticity curve constrains results. The individual parameter ranges allowed in this scenario are given in Table 2. In addition to v P and v S , the Poisson's ratio for each layer is also allowed to cover a wide range, from 0.2 to 0.5, equal to v P /v S larger than 1.633. No low-velocity layers are considered, i.e. velocities are required to increase with depth. Density is fixed to 1600 kg/m 3 in the individual layers and 2000 kg/m 3 in the  bedrock due to the lower sensitivity of the ellipticity curve to this parameter. This parameterization defines a rather large parameter space and, specifically in the inversions with a large number of degrees of freedom, more than 5 inversion runs need to be considered to ensure a good coverage of the solution space. We have two data sets to consider: one is derived from the simulation with only fundamental mode sources (Fig. 5b), the other from the more general simulated wavefield that also contains higher modes (Fig. 5e). Though the fundamental mode curves derived in both cases are highly similar, and both give good approximations of the theoretical ellipticity, the frequency range covered and the estimated uncertainties are slightly different. We thus show inversions of the fundamental mode data for both cases for comparison.
Starting with the fundamental mode simulation, we first invert only the right flank of the measured ellipticity curve. In the case of uniform layers, the minimum misfit found during the inversion drops strongly with increasing number of layers up to three layers over a halfspace (11 degrees of freedom) and stays approximately constant for further increases in layer number (Fig. 7a). The same is true for the parameterizations with a linear or powerlaw velocity function in the topmost layer. However, in these cases the plateau of the misfit function is already reached for two layers over the halfspace (10 degrees of freedom). The AICc reaches its minimum at the start of the plateau area in misfit for each of the different parameterizations, and the minimum values show only minor differences for the various parameterizations (Fig. 7c). The lowest value of AICc is found for the power-law parameterization in the topmost layer, which agrees best with the actual model used to calculate the synthetics.
A more detailed investigation in the case of uniform layers shows that a single layer over a halfspace is not able to provide an acceptable fit to the measured data (Fig. 8a, d). In the case of three layers over a halfspace, though, all models with a misfit of less than 0.4 are within the standard deviation of the data, and adding more layers only results in a smoother transition to higher velocities between 8 and 40 m depth (Fig. 8e, f). Though three is the true number of layers in the model, the inversion shows that a wide range of possible parameter values can explain the data and the layer thicknesses and velocities can deviate significantly from those of the true model. The trade-off between layer velocity and thickness is clearly apparent in all cases, as a wide range of possible regolith S-wave velocities between 50 and Fig. 9 Velocity profiles and fit to the data derived from the fundamental mode wavefield for the best parameterization, corresponding to the minimum AICc, in an unconstrained model space (Table 2). (a) Inverting the upper flank of the fundamental mode ellipticity curve only; (b) inverting both flanks of the fundamental mode ellipticity curve. Mode space and data are drawn as in Fig. 8. The color scale is the same for all subplots. In every case, all models with a misfit of less than 0.44 are judged to satisfy the data 550 m/s corresponds to layer thicknesses between 8 and 20 m. For the inversion with the best model parameterization, S-wave velocities in the regolith span the same range, whereas velocities below 11 m depth are not constrained at all (Fig. 9a). However, this does narrow down the original range of values for the regolith (S-wave velocity between 50 and 2500 m/s and layer thickness of 1 to 50 m) significantly. Thus, a low velocity layer is required by this small part of the fundamental mode curve already, as well as an additional layer of intermediate velocities above the halfspace, but details cannot be constrained reliably. The velocity and also the depth to the lowermost layer, considered as bedrock, is basically unconstrained.
Including the left flank of the fundamental mode ellipticity curve leads to comparable trends for misfit and AICc (Figs. 7b and d). Again, the lowest value of AICc is obtained for two layers over a halfspace with a power-law increase in the topmost layer, with a larger Fig. 10 Trade-off between the S-wave velocity at the bottom of the regolith layer and the regolith thickness for models resulting from the inversion of both flanks of the fundamental mode ellipticity curve using the optimum parameterization (Fig. 9b). White star marks true values for the baseline model (Table 1) offset between the minimum for this parameterization and the ones for a linear increase or a constant velocity. The power-law parameterization is thus more clearly favoured in this case (Fig. 7d). The trade-off between depth and velocity still allows for a large range of models, with regolith thicknesses between 6 and 26.5 m corresponding to S-wave velocities between 135 and 575 m/s (Figs. 9b, 10). Though the misfit curve is not sampled completely by our limited number of inversion runs in the huge parameter space considered, its nearly linear shape is readily apparent, and the baseline model lies on the curve (Fig. 10).
S-velocities are somewhat better constrained than when inverting only the right flank of the curve (Fig. 9). Specifically, the models provide a rough estimate of the S-wave velocity between 10 and 20 m depth, whereas using less data, S-wave velocities beneath 11 m depth are not constrained at all. In both cases, though, P-velocities at depth larger than 6-8 m are unconstrained. This can be understood from analysis of the sensitivity kernels for our baseline model (Fig. 11). The main sensitivity of the data is to velocities in the upper 10 m of the model, and generally sensitivity to S-wave velocity is an order of magnitude larger than to P-wave velocity for the fundamental mode. However, close to the fundamental mode ellipticity peak, where sensitivity to changes is highest, the right flank also shows increased sensitivity to the P-wave velocity of the upper 10 m. The left flank of the curve provides additional sensitivity to the regolith S-wave velocities, specifically between 5 and 10 m depth, and adds some, though limited, sensitivity to structure at larger depth up to 20 m, whereas influence of structure below this depth is very small. Accordingly, in both cases, some models that fit the data well only show S-wave velocities above 2000 m/s indicative of bedrock below 40 m depth. The minimum depth to the bedrock is estimated at 18 m from inversion of both flanks of the dispersion curve, compared to 11 m when inverting the right flank only.
Considering the data derived from the multi-mode wavefield, results are very similar when inverting only the fundamental mode ellipticity curve (Fig. 12a), showing that minor changes in data quality and curve picking do not influence the inversion results. Including higher mode data slightly improves the picture (Fig. 12b). The number of data points de- Fig. 11 Numerical approximation of sensitivity kernels at different frequencies for the fundamental mode (A-E) and first higher mode (F-J) ellipticity curve. Solid lines are for S-wave velocity and dashed lines for P-wave velocity. Note the different amplitude scaling at different frequencies. A-3.5 Hz, B-4.25 Hz, C-4.6 Hz, D-6 Hz, E-8.5 Hz, F-5 Hz, G-6 Hz, H-8.25 Hz, I-12 Hz, J-18.5 Hz scribing parts of the fundamental mode curve and the higher mode curve is almost equal (58 vs. 57), implying an equal weight given to fitting the different modes during the inversion. The lowest AICc is in this case found for a model with four layers of constant velocities over the half-space. Still, the S-wave velocity in the regolith is only constrained between 50 and 450 m/s, and the regolith thickness between 10 and 20 m, whereas the S-velocity in the layer below is barely constrained at 500 to 2500 m/s. Interestingly, the P-wave velocity below the regolith is not completely unconstrained in this case, but estimated to lie between 800 and 4250 m/s.

Fig. 12
Velocity profiles and fit to the data derived from the multi-mode wavefield for the best parameterization, corresponding to the minimum AICc, in an unconstrained model space (Table 2). (a) Inverting both flanks of the fundamental mode ellipticity curve. (b) Inverting both flanks of the fundamental mode and both parts of the higher mode ellipticity curve. Model space and data are drawn as in Fig. 8. The color scale is the same for all subplots. In every case, all models with a misfit of less than 0.44 are judged to satisfy the data The sensitivity kernels indicate that, unlike the case for dispersion curves, where higher modes sample deeper structure than fundamental modes at the same frequency (e.g. Rivet et al. 2015), the higher mode ellipticity curve provides only very limited information on deeper structure and, like the fundamental mode curve, is mainly sensitive to the regolith layer (compare kernels at points D and G in Fig. 11). However, the higher-frequency parts of the higher mode curve (points G-J in Fig. 11) are sensitive to changes in P-wave velocity on a comparable scale as to changes in S-wave velocity, though sensitivity to changes in v S is still largest. In contrast, for the larger part of the fundamental mode curve (points A-C and E in Fig. 11), sensitivity to v P is an order of magnitude smaller than sensitivity to v S .
In summary, even when there are no a priori constraints on the parameter space, the inclusion of higher mode data results in somewhat tighter constraints on the model space, specifically on P-wave velocities at depths between 10 and 30 m. Independent of which part of the data is inverted, the data demand low velocities at shallow depth (v S between 50 and 550 m or 450 m, depending on the inclusion of higher mode data), at least one additional layer of intermediate velocities between the regolith and the bedrock, and an increase to bedrock velocities (v S larger than 2000 m/s) only below 16 m depth, but more detailed conclusions cannot be drawn.

Constrained Parameter Space
In a second step, we introduce a priori constraints to the parameter space and again invert an increasing number of data, from the right flank of the fundamental mode peak only via both of its flanks to the inclusion of higher mode information. The allowed parameter range for the inversions is given in Table 3. Constraints are mainly introduced for the top-most layer, the regolith, and the halfspace at the bottom of the layer, requiring high velocities appropriate for basalt. For the regolith, velocities are constrained based on the results of the laboratory measurements on soil analogues, using the values found at pressures corresponding to 0 and 20 m depth and adding an additional 20 % of uncertainty to the lowest value to include the effect of reduced surface pressure on Mars, and 10 % to the highest value. The regolith thickness is constrained by the information based on rocky crater ejecta analysis and fragmentation theory, and the lower limit set to 5 m, as InSight's heat flow probe HP 3 is expected to penetrate to this depth within 30 days after deployment (Spohn et al. 2012;Grott et al. 2015). If HP 3 penetration encounters no difficulties, it can be assumed that the probe is moving through unconsolidated regolith only, whereas contact with hard rocks would impede the penetration. Thus, the maximum penetration depth of HP 3 can serve to constrain minimum regolith thickness, assuming that other causes that could prevent deeper penetration like non-vertical motion of HP 3 or instrument malfunction can be excluded. Additional analysis of seismic data generated by HP 3 hammering ) could further constrain regolith properties (see below). The parameters of the intermediate layers have larger ranges, in accordance with the limited available prior information. However, their thickness is restricted between 5 and 30 m, based on the reasoning that layers as thin as 1 m cannot be meaningfully constrained by the data, whereas geological information gives a shallow depth for the basalt layer, not supporting overlying regolith greater 30 m. Besides, v P /v S is required to be between 1.63 and 2.08 in the upper layers (meaning Poisson's ratio between 0.2 and 0.35) and between 1.63 and 1.87 in the basalt (translating to Poisson's ratio between 0.2 and 0.3). Inversion results based on the fundamental mode ellipticity curve again favour a model with two layers over a halfspace with a power-law velocity increase in the first layer. Constraining the velocities in the near-surface layer, the parameter to which the ellipticity curves are most sensitive (Fig. 11), helps to put tighter constraints on other parameters, i.e. the thickness of the near-surface layer and the velocity in the layer below, and the depth to the bedrock (Fig. 13). All of these parameters are slightly more tightly constrained if both flanks of the fundamental mode peak are inverted, specifically the depth values (Fig. 13b). The thickness of the regolith layer is estimated at 8 to 12.5 m, versus 6.5 to 12 m when only the right flank is inverted, and the depth to the bedrock is estimated at 17 to 39 m, compared to 16.5 to 42 m. The true P-wave velocity in the layer below the regolith, 1500 m/s, is almost at the center of the possible values obtained from the inversion, which lie between 750 and Fig. 13 Velocity profiles and fit to the data derived from the fundamental-mode wavefield for the best parameterization, corresponding to the minimum AICc, if the model space is constrained according to Table 3. (a) Inverting the upper flank of the fundamental mode ellipticity curve only; (b) inverting both flanks of the fundamental mode ellipticity curve. Model space and data are drawn as in Fig. 8. The color scale is the same for all subplots. In every case, all models with a misfit of less than 0.44 are judged to satisfy the data 2300 m/s. Values for v S between 450 and 1200 m/s show a tendency to overestimation, compared to the true value of 790 m/s. The bedrock velocity is only constrained by the limits on the parameter space, as can be expected from the missing resolution of the data at these depths (Figs. 11, 4c).
The fundamental mode curve extracted from the multi-mode data set has a slightly larger standard deviation at a number of frequencies. As a result, inverting these data using the same parameterization as above results in larger model uncertainties, specifically with regard to the velocities in the sub-regolith layer (Fig. 14a). Again a model with two layers over the halfspace and a power-law velocity increase in the shallowest layer is preferred in all three cases discussed in the following. The regolith thickness is estimated at 6.5 to 13.5 m, and the S-and P-wave velocities in the layer below cover a larger range, 500 to 1500 m/s and 800 to Fig. 14 Velocity profiles and fit to the data derived from the multi-mode wavefield for the best parameterization, corresponding to the minimum AICc, if the model space is constrained according to Table 3. (a) Inverting both flanks of the fundamental mode ellipticity curve; (b) inverting both flanks of the fundamental mode ellipticity curve and the low-frequency part of the higher mode curve; (c) inverting both flanks of the fundamental mode ellipticity curve and both parts of the higher mode curve. Model space and data are drawn as in Fig. 8. The color scale is the same for all subplots. In every case, all models with a misfit of less than 0.44 are judged to satisfy the data 2500 m/s, respectively, with a tendency to overestimate the true velocities. The depth to the bedrock likewise has a higher uncertainty; it lies between 16 and 42.5 m. Including the lowfrequency part of the higher mode curve leads to limited improvements, with the possible parameter range for the depth of the regolith layer and the depth to the bedrock slightly smaller at 9.5 to 13 m and 20 to 42.5 m (Fig. 14b). Inverting all available higher mode information together with the fundamental mode curve has a stronger effect: discontinuity depths and sub-regolith velocities are clearly more tightly constrained (Fig. 14c). Regolith thickness is estimated at 8 to 12 m, which brackets the true thickness of 10 m symmetrically. The depth to the bedrock is constrained at 16 to 31.5 m, which roughly corresponds to the depth between the bottom of the coarse ejecta layer and the top of the basalt layer. The inversion results thus indicate potentially higher velocities than above in this depth range, the fractured basalt layer. Velocities in the sub-regolith layer lie between 400 and 900 m/s and 700 and 1700 m/s, which does include the true velocities of 790 m/s and 1500 m/s, but shows a tendency to underestimate the true value. The P-wave velocity of the bedrock is indeed underestimated, whereas the S-wave velocity is at the upper limit of the inverted range. Similar results are obtained if only the high-frequency part of the higher mode data is inverted together with the fundamental mode. Analysis of seismic signals generated by the hammer strokes of HP 3 can potentially further constrain regolith properties, with an estimated uncertainty in regolith thickness of 20 % ). An average P-wave velocity of the regolith could also be obtained from HP 3 signal stacking at different penetration depths . Here, we also assign a 20 % uncertainty to the resulting P-wave velocity, to capture the uncertainty of the measurement but also to allow for a velocity increase with depth within the regolith due to compaction. The resulting, tighter constraints on the parameter space are outlined in Table 4.
When inverting the fundamental mode data only, in contrast to the previous inversion runs, the uniform models achieve significantly higher minimum misfit values, and consequently higher AICc's, than models with a velocity increase in the first layer for all tested degrees of freedom. A power-law velocity increase in the first layer is again favoured by the comparison of inversion results (minimum misfit and minimum AICc). The further constraints on regolith velocities result in tighter constraints on layer thickness and velocities (Fig. 15a). For the regolith thickness, values between 8 and 11 m are obtained, whereas values for the depth to the bedrock range between 19.5 and 35 m, which brackets the depth to the bottom of the coarse ejecta layer and the top of the basalt layer. Velocities in the subregolith layer are estimated at 450 to 1150 m/s and 850 to 2200 m/s, respectively, which is consistent with the true v S ± 350 m/s and the true v P ± 675 m/s. Including higher mode data puts tighter constraints on discontinuity depths, at 8.5 to 10.5 m for the regolith thickness and 17 to 26.5 m for the lower limit of the second layer, which contains the true depth to the interface between coarse ejecta and fractured basalt at 20 m depth (Fig. 15b). Velocities in the sub-regolith layer are again underestimated, though, at 450 to 850 m/s and 750 to 1600 m/s, respectively, and the same is true for the bedrock velocities. Again, similar results are obtained when only inverting the high-frequency part of the higher mode curve together with the fundamental mode data.
In summary, placing some basic constraints on the regolith and bedrock velocities results in some useful information on the regolith thickness, estimated with 20 % uncertainty, and  Fig. 14c). Model space and data are drawn as in Fig. 8. The color scale is the same for all subplots. In every case, all models with a misfit of less than 0.46 are judged to satisfy the data the velocities in the layer below, estimated with uncertainties of 50 % or less. Including the complete available higher mode information leads to tighter constraints and improved estimates of discontinuity depths. The estimated velocities, however, show a tendency to lie below the true model velocities. A possible reason for this is that the actual velocity structure in the upper 10 m contains more details than the model parameterization can capture in 5 layers and deviates from a power law, and trade-offs with properties of the layer below result from this. Comparing the results when inverting different parts of the dispersion curve, e.g. the fundamental mode only and both fundamental and higher mode, might help to identify this kind of bias: the high-frequency part of the first higher mode is very sensitive to details of the velocity in the regolith layer (Fig. 11 I-J) and might respond strongly to unmodelled complexity. Including potential information from seismic analysis of HP 3 data will lead to tighter constraints on the sub-regolith structure, including both velocities and layer thickness, especially if only fundamental mode data is available. The improvement is smaller if both fundamental and higher mode data are inverted as this already results in significantly better retrieval of model parameters in the less-constrained case.

Regolith Properties
We investigate how variations of the parameters of our velocity model within a reasonable range influence the measured Rayleigh wave ellipticity. One obvious parameter to vary is the regolith thickness, as it is only constrained to lie between 2.4 and 17 m from orbital imagery and crater ejecta analysis. Figure 16 shows the influence of variations in regolith thickness on the ellipticity curves. The thickness of the coarse ejecta layer and the fractured basalt layer are scaled in accordance with the regolith thickness in each model, based on the assumption that the thickness of each of these layers is related to the cratering history. Thus, areas with more cratering should have a thicker regolith as well as a thicker coarse breccia and fractured basalt layer. Peak frequencies of the ellipticity curve vary from 16.3 Hz for 3 m of regolith to 3.3 Hz for 15 m of regolith. The frequency decrease with thickness is less pronounced for thicker layers. The shape of the individual ellipticity curves for different modes does not vary. It is thus reasonable to assume that the processing scheme outlined above for the baseline model will also work for all regolith thicknesses within the range expected at the InSight landing site and that the peak frequency can give a direct estimate of regolith thickness when regolith velocity is reasonably well known (Fig. 10).
Other parameters with some uncertainty are the velocities in the regolith and in the coarse ejecta layer. The regolith velocities are based on laboratory measurements at reference pressures. Ambient pressure on the surface of Mars will be lower, though, which could also lower the regolith velocities at the surface, before the effect of soil compaction becomes dominant (Murdoch et al. 2016). P-wave velocities of lunar regolith obtained from active seismic experiments at Apollo 14 and 16 are even lower than the ones used here, with values close to 100 m/s to more than 10 m depth (Kovach et al. 1971(Kovach et al. , 1972Cooper and Kovach 1974). However, the low velocities of lunar regolith, like the high Q, are at least partly due to its extreme dryness after outgassing under vacuum conditions (Tittmann 1977;Tittmann et al. 1979), removing even adsorbed water between mineral grains, which did not occur on Mars. Accordingly, we expect that  (Table 4). (a) and (c) inversions of complete data set for baseline model (compare Fig. 14b); (b) and (d) inversions of complete data set for model with constant velocity in the regolith layer (compare Fig. 18). In each figure, the blue line and dots correspond to a parameterization with uniform layers, the red line and dots to a parameterization with a linear velocity increase in the topmost layer, and the green line and dots to a parameterization with a power-law velocity increase in the topmost layer the lab data provide the best available approximation of the actual regolith velocities on Mars.
The inversion results for the baseline model seem to indicate that it is possible to distinguish between a constant velocity and a velocity increase within the regolith layer based on the ellipticity curve. In all above inversions in which the parameter space is constrained, and even in almost all unconstrained inversions, a model with a power-law velocity increase in the regolith layer is preferred, in keeping with the actual velocity structure of the baseline model. To corroborate this finding, we performed a multi-mode wavefield simulation for a model with a constant regolith velocity, averaged over the velocity increase in the baseline model, and analyze the resulting data. Again, the different parts of the higher mode curve are clearly discernible in the data. Inverting all extracted data for a model space constrained by HP 3 results (Table 4) results in a preferred model of two constant-velocity layers over a halfspace (Figs. 17b and d). The misfits and AICc values for all model parameterizations follow a common trend if two or more layers over the halfspace are considered, but the AICc minimum is clearly found for a constant regolith velocity. This is in contrast to the results for the baseline model (Figs. 17a and c), where models with a constant velocity in the topmost layer show an inferior fit to the data for low numbers of degrees of freedom below ∼ 14. The thickness of the regolith is estimated at 8.5 to 11.5 m (Fig. 18). Similar to the results for the baseline model, P-wave velocities directly below the regolith tend to be underestimated at 900 to 1850 m/s, but the estimated S-wave velocities nicely constrain the true value of 790 m/s at 550 to 1050 m/s. This also corroborates our interpretation that the underestimation of velocities in the sub-regolith layer for the baseline model when inverting a similar data set (Fig. 15c) is mainly due to unmodelled structure in the regolith layer. In case of the model with a constant regolith velocity, the structure of the regolith layer can be perfectly matched during the inversion, and here, no underestimation occurs. The thickness of the sub-regolith layer is constrained to 22 to 39 m, i.e. slightly overestimated.

Fig. 18
Velocity profiles and fit to the data for the best parameterization, corresponding to the minimum AICc, for a model with a constant velocity in the regolith layer if the model space is constrained according to Table 4. Both fundamental mode and higher mode ellipticity curves are inverted. Model space and data are drawn as in Fig. 8. All models with a misfit of less than 0.44 are judged to satisfy the data

Sub-regolith Velocity
In contrast to the regolith velocity, the seismic velocities of the coarse breccia layer in the baseline model is not constrained by dedicated lab tests, but only by a terrestrial analog. On the Moon, sub-regolith P-wave velocities between 250 and 300 m/s have been found in the active seismic experiments at Apollo 14 and 16 (Kovach et al. 1971(Kovach et al. , 1972 as well as in the seismic profile at Apollo 17 , and interpreted as indicating brecciated, shattered basalt or impact debris (Kovach et al. 1972Cooper and Kovach 1974). The seismic profiling experiment constrained the thickness of this layer to more than 200 m . Again, the very low velocities are related to the dryness of the material, which is specific to the vacuum conditions on the Moon. Besides, recent high-resolution lunar gravity data find a high porosity of the lunar crust, on average 12 %, to a depth of a few kilometers . This high porosity also results in significantly lower seismic velocities than expected for Mars. Still, somewhat lower velocities in the coarse ejecta layer than assumed in the baseline model cannot be excluded based on the information presently available.
Reducing the assumed P-wave velocity of the coarse ejecta layer, and scaling all other parameters accordingly, will reduce the peak frequency of the ellipticity curve, but can also change the shape of the curve itself (Fig. 19). If the velocity decreases past a certain threshold, the impedance contrast between the coarse ejecta layer and the fractured basalt layer becomes also relevant for the shape of the ellipticity curve: it shows two peaks that are associated with the contrasts at the top and bottom of the coarse ejecta layer, and these narrow peaks have rather high amplitudes. The resolution of actual data will likely be insufficient to separate the two narrow peaks and outline the amplitude minimum inbetween, though. We investigate this further by computing a synthetic multi-mode data set for a model with a P-wave velocity of 650 m/s in the coarse ejecta layer, and all other parameters in this layer down-scaled accordingly. Figures 20a and b show the standard H/V and HVTFA results for this model, respectively. Both fundamental mode curves show a broad peak rather than imaging the two adjacent peaks with a minimum inbetween. Compared to the baseline Black lines are theoretical ellipticity curves for the fundamental and first higher mode. Dashed curves give H/V (blue), ellipticity (green), and theoretical curves (gray) for baseline model model, the peak is shifted to lower frequencies. HVTFA again provides a more reliable estimate of the ellipticity, especially of the left flank of the fundamental mode peak. The higher mode is less well resolved in this case-the ascending branch at frequencies below 5 Hz can only be measured rather poorly, whereas the broad plateau around 10 Hz can be recovered reasonably well.

Fig. 21
Velocity profiles and fit to the data for the best parameterization, corresponding to the minimum AICc, for a model with reduced velocities in the coarse ejecta layer if the model space is constrained according to Table 4. (a) Inverting both flanks of the fundamental mode ellipticity curve; (b) inverting both fundamental mode and higher mode ellipticity curves. Model space and data are drawn as in Fig. 8. The color scale is the same for all subplots, and all models with a misfit of less than 0.3 are judged to satisfy the data Inversion of the HVTFA results indicates that a model with two layers over a half-space and a power-law velocity increase in the uppermost layer explains the data best. When inverting the fundamental mode data only, the results show a high uncertainty, likely because important data points around the peak frequency could not be extracted reliably from the HVTFA results and are missing (Fig. 21a). Compared to the results for the baseline model (Fig. 15a), there is a trend to lower velocities in the sub-regolith layer, though. Indeed, the model space that is compatible with the data extends to the lower boundary of the allowed parameter space. Including the higher mode information (Fig. 21b) results in significantly tighter constraints, with velocities in the sub-regolith layer estimated at 340 to 635 m/s and 600 to 1250 m/s compared to 300 to 1050 m/s and 500 to 1950 m/s. Besides, the thickness of the sub-regolith layer is better constrained at 22 to 30 m, compared to 21 to 42 m. Again, the velocities in the sub-regolith layer extend very close to the parameter space boundary. An underestimation of velocities in this layer, as observed for the baseline model (Fig. 15b), is prevented by these constraints. When encountering cases like this with actual data, it would be advisable to extend the parameter space to lower velocities, if no independent, prior constraints are available, to fully capture model uncertainty. Here, we wanted to investigate if we can differentiate between different models using a consistent parameterization, though, and indeed, results point to lower velocities in the sub-regolith layer in this case.

Attenuation
Finally, it has been noted that Q can have a strong influence on H/V amplitudes (Lunedei and Albarello 2009;Dal Moro 2015). However, the fundamental mode ellipticity curve for our model stays constant, regardless of increased Q values in the regolith and sub-regolith layers (Fig. 22). Some changes occur in higher mode curves if Q is increased, which is also mirrored in higher amplitudes of the higher order SH resonances (Fig. 22). Changes between the models with a Q of 100 or 500 in the regolith are negligible, though, as also observed by Dal Moro (2015) in H/V curves for the Moon. We also computed and analysed synthetic seismograms for the model with a regolith Q of 100 (Figs. 20c and d). The right flank of the standard H/V curve is actually closer to the theoretical ellipticity in this case, and the peak amplitude is larger by about a factor of two compared to the baseline model, consistent with the observations of Albarello (2009) andDal Moro (2015). In contrast, the HVTFA curve is indistinguishable from the one of the baseline model. This indicates that HVTFA indeed provides an estimate of the Rayleigh wave ellipticity, whereas the standard H/V does not, and cannot be used to constrain damping. The higher mode branch is barely visible in the HVTFA results, compared to the good visibility in case of low Q (Fig. 6b), and can only be traced at high frequencies between 9 and 18 Hz. If, in a realistic situation, only the fundamental mode curve could be reliably extracted due to high Q values, or little Rayleigh wave excitation at high frequencies, this would probably lead to higher model uncertainty in the inversion results (compare Figs. 15a and b).

Discussion
SEIS data will be recorded at sampling rates of 20 Hz for the three-component VBB sensor and 100 Hz for the three-component SP sensor, respectively. Due to limits on downlink capacity for a mission that has to use and share existing orbiters around Mars, continuous seismic data will be relayed to Earth at a reduced rate, though, and full-range data will only be available upon requesting specific event time-windows, with a maximum volume of 8 Mbit per sol. Selection of these time windows has to be based on the available continuous data streams. The background model and its variations studied here show peak amplifications at frequencies between 3 and 17 Hz, which are above the sampling frequency of SEIS's continuously transmitted three-component VBB recordings at 2 Hz. A combined VBB/SP channel will be continuously transmitted at 10 Hz, but as this is a vertical channel only, no site amplification effects have to be expected. In an extreme scenario, using the maximum expected regolith thickness of 17 m and the extremely low regolith velocities measured on the Moon, the fundamental mode ellipticity peak would be around 1 Hz, within the range of continuously transmitted VBB data. This is however considered highly unlikely for In-Sight as the very low lunar regolith velocities could only be reached by outgassing in a hard vacuum, which did not occur under the atmospheric conditions of Mars. Besides, a regolith thickness as large as 17 m is probably the exception rather than the rule within the selected landing ellipse (Warner et al. 2016b;. Thus, we do not expect any visible influence of regolith resonances on teleseismic recordings by InSight. If the regolith is thicker than 10 m, frequencies below 5 Hz, which could be important for regional events, are likely to be affected, though, whereas recordings of local events would probably also show site amplification for a smaller regolith thickness. In each considered scenario, site amplification would be measurable with the full rate SP data sampled at 100 Hz, whereas the full rate VBB data sampled at 20 Hz would likely show site effects for regolith thicknesses larger than 5 m. Amplitudes of the horizontal components can then be expected to be considerably higher than the vertical component amplitude around these frequencies, as observed in the Apollo lunar seismic data (e.g. Lammlein et al. 1974).
In the absence of oceans, a main source of ambient noise on Mars is expected to be the direct interaction between the atmosphere and the solid surface of the planet. Indeed, wind noise is the primary signal registered by the Viking 2 seismometer (Anderson et al. 1977;Nakamura and Anderson 1979). However, in that case the seismometer was placed on top of the lander, not on the Martian surface, and the wind was transmitted to the sensor via interaction with the lander and not through the ground. Murdoch et al. (2016) estimated the wind environment at frequencies below 1 Hz at the InSight landing site and studied how this may influence mechanical noise transmitted to SEIS through the ground by wind interacting with the InSight lander. On Earth, with its much denser atmosphere, wind has also been identified as a direct source of ambient noise at higher frequencies, including the range used here to study ellipticity (Withers et al. 1996;Mucciarelli et al. 2005;Naderyan et al. 2016). Both Withers et al. (1996) and Naderyan et al. (2016) specifically found wind effects in seismic recordings at locations with little topography and vegetation, similar to the InSight landing site . Quiros et al. (2016) observe Rayleigh waves generated by wind gusts along a geophone line at frequencies between 1 and 10 Hz. A high frequency component has also been observed in the seismic recording of a dust devil on a terrestrial desert playa (Lorenz et al. 2015) and identified as surface waves (Kenda et al. 2016). Dust devils are a frequently observed phenomenon on Mars, and dust devil tracks have been mapped from orbital imagery in the InSight landing region (Reiss and Lorenz 2016). Modeling based on large eddy simulations also indicates that a number of seismically observable vortices might occur in the landing area during daytime (Kenda et al. 2016). Thus, at least during high-wind regimes, winds on Mars can be expected to generate a background noise wavefield that could be used for ellipticity measurements. Though the InSight landing site was selected to lie outside of storm tracks, reducing the overall wind noise, the landing will take place during the later part of the global dust storm season ). This could result in favourable conditions for ambient vibration based measurements directly after landing, meaning ellipticity measurements could be available simultaneously with HP 3 hammering results, ideal for a joint interpretation. Murdoch et al. (2016) analyzed the influence of dynamic pressure and winds on lander mechanical noise transmitted to SEIS at frequencies below 1 Hz. In this study, we neglect the potential influence of high-frequency (> 1 Hz) lander generated noise, e.g. due to eigenmodes of solar panel vibrations in response to wind load, in favour of first understanding general first-order effects in a homogeneous background wavefield. Potentially, the lander as a very close source could be an important contributor to the noise wavefield at 1-20 Hz. How it may affect the observed spectra and spectral ratios remains subject of future study.
On the Moon, constant bombardment by small meteorites in combination with long coda duration has been predicted to create a continuous background "seismic hum", though with amplitudes too low to be observable by Apollo seismometers (Lognonné et al. 2009). This source of a background ambient vibrations wavefield can likely be excluded for Mars, though. In contrast to the situation on the Moon, the Martian atmosphere affects meteoroids by deceleration, ablation, and fragmentation, resulting in a minimum meteoroid size to reach the Martian surface (Popova et al. 2003). Taking into account current crater production functions and wave propagation characteristics on Mars, Teanby (2015) finds only 1-3 regional impacts per year with signal amplitudes in the 1-16 Hz frequency range above the SP noise floor for InSight. In addition, coda length due to scattering is expected to be significantly reduced on Mars compared to the Moon.
Another possible source of ambient vibration surface waves at the InSight landing site are thermal events. On the Moon, numerous high-frequency weak repeating events have been observed locally at all Apollo sites between sunrise and sunset and attributed to variations in diurnal thermal stresses (Duennebier and Sutton 1974). The events have been estimated to occur within 4 km distance of each station, with slumping of small amounts of soil triggered by thermally induced stresses as suggested source mechanism (Duennebier and Sutton 1974). Duennebier (1976) was able to determine source locations with the help of the Lunar Surface Profiling Event of Apollo 17 and found that sources do not correlate with steeply dipping surfaces, but seem to be associated with craters, implying that gravitational energy is not necessary to trigger thermal moonquakes. Criswell and Lindsay (1974) suggest a different source mechanism akin to booming dunes in terrestrial deserts, which emit a characteristic sound during slumping. Rayleigh waves extracted from ambient vibration cross-correlations at the Apollo 17 array between 3.5 and 11.5 Hz, the only dispersive surface waves ever recorded in the highly scattering environment of the Moon, are based on thermal quakes as ambient vibration sources (Larose et al. 2005;Tanimoto et al. 2008;Sens-Schönfelder and Larose 2010).
Though the InSight landing site was selected to have low slopes, there are numerous small craters (Warner et al. 2016b) that could be locations of thermally triggered soil slumping. Temperatures at the lunar surface may increase by almost 300 K during solar heating, with the largest difference of about 200 K within 24 h occurring during sunrise (Langseth et al. 1973). In contrast, diurnal temperature variations measured by the Mars Science Laboratory rover in Gale crater approximately 550 km south of the InSight landing site are only about 90 K (Hamilton et al. 2014), however at a significantly higher thermal inertia (300-350 Jm −2 K −1 s −1/2 ) than measured for the InSight landing site (200 Jm −2 K −1 s −1/2 , . Diurnal temperature variations could thus be several 10s of K larger at the landing site, with the complete temperature increase occurring during less than 10 hours. As the details of the temperature variations on the Moon are not resolved by taking just one data point per day, it is not clear if the expected temperature changes on Mars will approach lunar values when considering comparable time scales and the amount of resulting thermal stress will be sufficient to trigger thermal quakes. Thus, while thermal quakes are likely to generate high frequency surface waves, their abundance on Mars near the InSight landing site remains to be determined. In this synthetic study, 30 minutes of data were sufficient to estimate ellipticity curves. However, this of course depends on the source distribution and activity. Thus, for real data, a longer time span of measurements, up to several hours, might be desirable, especially if surface sources generating Rayleigh waves are infrequent on Mars, to obtain better statistics in the HVTFA evaluation and be able to resolve the contribution of different modes. The data do not necessarily have to be acquired as one continuous recording, though, but could in principle also consist of several shorter time spans combined for analysis. This allows more flexible operations, e.g. first requesting a shorter time period and, if this turns out to be insufficient, backing results with additional data recorded later. In the synthetic test, HVTFA provided the best estimates of the actual ellipticity. For a measured data set, it will however be useful to apply both HVTFA and RAYDEC, as suggested by Hobiger et al. (2012), to get an idea about the reliability of results and to identify potential contamination by higher modes.
In our inversion tests for fundamental mode data, the tightest constrains on sub-regolith properties can be obtained when using information on the regolith that might be obtained from analysis of SEIS recordings of HP 3 hammering. In this case, the information from HP 3 reduces model uncertainty significantly and it would be highly desirable to use it (compare Fig. 14a and Fig. 15a). An alternative could be the inclusion of higher mode information in the inversion, especially at high frequencies, where only little Rayleigh wave energy might be available, though. In that case, the improvement between inversions for a parameter space constrained based on general a priori information or constrained based on HP 3 analysis is less significant (compare Fig. 14c and Fig. 15b) and results are clearly better constrained than when using fundamental modes only. However, it might be useful to compare models derived from different parts of the data to identify potential biases in velocity estimations. Other constraints that could be used to a priori reduce the size of the parameter space in the inversion are a map of the regolith thickness derived from rocky ejecta craters that should be available for the whole landing ellipse before landing . Once the actual landing position has been determined, this map will provide an initial estimate of regolith thickness. In addition, fragmentation theory can provide an estimate of the maturity index of the regolith from crater counts and surface rock abundance within the landing region, which will allow an independent estimate of regolith thickness. Results from these two methods could be used to narrow down the possible range of regolith thickness at the landing site and constrain this parameter in the ellipticity inversion, reducing the tradeoff between velocity and thickness for the regolith layer (Fig. 10) and leading to tighter constraints on sub-regolith properties.
If only fundamental mode data are available, using both flanks of the ellipticity peak provides tighter constraints than using the right flank only. Specifically, data samples close to the peak itself are useful in constraining the solution, if they can be estimated reliably. On the other hand, data gaps should also be taken seriously as they can result from superposition of different modes or from peak splitting, and interpolating through them can result in gross errors. If no well-founded prior constraints are available, it is also advisable to extend the parameter space to get a reliable estimate of model uncertainty if initial inversions extend to the parameter space boundary.
In the constrained inversions, ellipticity cannot only provide estimates of the sub-regolith structure (existence of additional layer over the halfspace, sub-regolith velocity, layer thickness), but also distinguish between a constant velocity in the regolith or a velocity increase with depth (Fig. 17). From the example shown above, it might be assumed that the ability to distinguish between the two cases rests on the availability of higher mode data, as they are influenced most distinctly by differences in the velocity structure of the top-most layer (Fig. 4c). There is a visible influence on the fundamental mode curve, too, though, around 8 Hz (Fig. 4c). To understand how much the resolution of regolith structure depends on the availability of higher mode data, we ran inversions for different model parameterizations using fundamental mode data only. Comparing misfits and AICc values for two layers over a halfspace, again the distinctly lowest values are obtained for constant regolith velocities. Accordingly, using AICc to rank models, the fundamental mode ellipticity curve is sufficient to distinguish between a constant-velocity regolith layer and one where velocity increases with depth. This is an information that will not be available from HP 3 analysis, and could for example neither be obtained from the active seismic experiments on the Moon, which only resulted in a sparsely sampled travel-time curves that were interpreted in terms of a stack of constant-velocity layers (Cooper and Kovach 1974).
Trans-dimensional inversions could be an alternative to the model ranking approach using AICc as shown here. It offers a way to directly include the dimension of the parameter space as a variable to be solved for in the inverse problem in a Bayesian framework (e.g. Malinverno 2002;Sambridge et al. 2006;Dettmer et al. 2012). However, it would be advantageous to keep some parameterization options that are specific to site characterisation, e.g. a power-law velocity increase in the sediment layer, and the possibility to constrain Poisson's ratio in addition to individual constraints on v P and v S , when inverting InSight data.
While H/V peak amplitudes can potentially help to distinguish between underground structure with high and low Q-values, the fundamental mode ellipticity, and correspondingly its approximation by HVTFA, are independent of Q. On the one hand, this removes additional complexity from the inversions, but on the other hand it means that the processing outlined here cannot help to determine the regolith Q. A complementary approach would be necessary, e.g. by using modeling based on the method of Arai and Tokimatsu (2004) and further developed by Lunedei and Albarello (2009). Results of an ellipticity inversion as described above could provide a priori constraints on velocities and layer thicknesses in that case. This method would also require assumptions on the relative contribution of Love and Rayleigh waves to the ambient wavefield, though. Another option could be based on diffuse field theory, which has been applied to H/V spectral ratios calculated from both ambient vibrations Kawase et al. 2015) and earthquake data (Kawase et al. 2011;Salinas et al. 2014). For ambient vibrations as well as for earthquakes located up to hundreds of kilometers away from the station, the results are compatible with a 3D diffuse field model that is sensitive to Q Salinas et al. 2014). While the application to earthquake data requires extensive stacking to approximate a diffuse field, which might not be possible on Mars due to limited seismic activity (e.g. Knapmeyer et al. 2006), ambient vibration data could potentially be modeled that way, again using prior constraints from ellipticity inversions. Besides, our modeling results indicate that, at least for the model range considered here, the clear visibility of a higher mode ellipticity curve over an extended frequency range could serve as a first-order indication against high Q values. The absence of a clear higher mode curve does not necessarily indicate a high Q, though, as variations in the model and the available sources might also influence the visibility of higher modes (e.g. Fig. 20b).

Conclusions
We constructed a plausible model of the shallow subsurface at the InSight landing site, based on laboratory measurements and analysis of orbital data, and investigate how parameters of this model and reasonable variations of it can be retrieved from ambient vibration Rayleigh wave ellipticity. We consider two different methods to calculate the ellipticity from the wavefield and find that, while both provide better estimates than the standard H/V curve, the method based on time-frequency analysis gives the best results in our case and also allows for the extraction of higher-mode information. This information proved subsequently very useful to constrain model parameters in an inversion. We use model ranking based on the AICc to select a preferred model parameterization and show that ellipticity data can distinguish between different velocity-depth functions in the shallowest layer, e.g. a constant regolith velocity or a velocity increase with depth due to compaction. This information might not be obtainable from other InSight data. Either the combination of fundamental and higher mode data and some reasonable a priori constraints on the parameter space or the use of fundamental mode data only and tighter constraints on regolith properties, i.e. from HP 3 hammering analysis or a priori regolith thickness maps, result in useful information on sub-regolith properties. While unconstrained inversions already give an idea about the existence of a low-velocity layer and additional layers above the bedrock, constrained inversions can determine velocities in the sub-regolith layer at uncertainties of less than ±400 m/s for v S and ±700 m/s for v P and the thickness of the sub-regolith layer to within a few meters. Ellipticity measurements cannot constrain Q, but could provide useful information about the subsurface model for wavefield or spectral modeling techniques that can. Either alone or in combination with results from seismic analysis of HP 3 signals, ellipticity measurements can provide important constraints on properties of the regolith and sub-regolith layers and an estimate of the minimum depth to the bedrock.