Crustal S-Wave Velocity from Apparent Incidence Angles: A Case Study in Preparation for InSight

Retrieval of crustal structure and thickness of Mars is among the main goals of InSight. Here we investigate which constraints on the crust at the landing site can be provided by apparent P-wave incidence angles derived from P-receiver functions. We consider receiver functions for six different Mars models, calculated from synthetic seismograms generated via Instaseis from the Green’s function databases of the Marsquake Service, in detail. To allow for a larger range of crustal thicknesses and structures, we additionally analyze data from five broad-band stations across Central Europe. We find that the likely usable epicentral distance range for P-wave receiver functions on Mars lies between 35∘$35^{\circ}$ and the core shadow, and can be extended to more than 150∘$150^{\circ}$ by also using the PP-phase. Comparison to models for the spatial distribution of Martian seismicity indicates that sufficient seismicity should occur within the P-wave distance range around InSight within the nominal mission duration to allow for the application of our method. Apparent P-wave incidence angles are derived from the amplitudes of vertical and radial receiver functions at the P-wave onset within a range of period bands, up to 120 s. The apparent incidence angles are directly related to apparent S-wave velocities, which are inverted for the subsurface S-wave velocity structure via a grid search. The veracity of the forward calculated receiver functions and apparent S-wave velocities is ensured by benchmarking various algorithms against the Instaseis synthetics. Results indicate that apparent S-wave velocity curves provide valuable constraints on crustal thickness and structure, even without any additional constraints, and considering the location uncertainty and limited data quantity of InSight. S-wave velocities in the upper half of the crust are constrained best, but if reliable measurements at long periods are available, the curves also provide constraints down to the uppermost mantle. Besides, it is demonstrated that the apparent velocity curves can differentiate between crustal velocity models that are indistinguishable by other methods.


Introduction
Since their introduction in the 1970s (Vinnik 1977;Langston 1979), receiver functions have become a standard tool to study the crustal and upper mantle structure of the Earth using the coda of distant earthquakes. As receiver function analysis is basically a single station method, it is also envisioned to be very useful during NASA's InSight mission, which will place a single seismometer, SEIS (Seismic Experiment for Interior Structure), on the surface of planet Mars in November 2018 (Banerdt et al. 2013;. Crustal thickness and structure of Mars are among the main scientific objectives of InSight. Orbital topography and gravity data can constrain relative variations in Martian crustal thickness, but need to assume the crustal thickness at a specific tie-point, usually a global thickness minimum in the Isidis impact basin (e.g. Neumann et al. 2004;Plesa et al. 2016). An independent estimate of crustal thickness at the InSight landing site could thus help to better constrain Moho depths all across Mars. The receiver function method has already been applied to extraterrestrial seismograms: S-receiver functions have been calculated for one station of the Apollo lunar seismic network, though their interpretation in terms of crustal thickness remains non-unique (Vinnik et al. 2001;Lognonné and Chenet 2003).
As a travel-time method, receiver functions generally suffer from a trade-off between the depth of a seismic discontinuity and the velocity above (Ammon et al. 1990). Besides, the seismic activity of Mars is expected to be significantly smaller than Earth's (Knapmeyer et al. 2006;Plesa et al. 2018), so the number of usable teleseismic events might as well be limited. A common strategy to resolve the depth-velocity trade-off in receiver function inversion is the joint inversion with Rayleigh wave phase velocities (e.g. Özalaybey et al. 1997;Du and Foulger 1999;Juliá et al. 2000). For a single-station situation like InSight, with no local dispersion measurements readily available, a joint inversion with Rayleigh wave ellipticity measurements is a more realistic option (Chong et al. 2016;. Alternatively, the measurement of frequency-dependent P-wave polarization offers a possibility to derive crustal S-wave velocities directly from receiver functions. As already demonstrated by Wiechert (1907), the apparent P wave incidence angle φ P at the free surface of a half space, i.e. the measured polarization of an incident P wave, is twice the angle φ S of the reflected SV wave. The measured P wave incidence angle is an apparent rather than the true one due to the superposition of motion of the incoming P wave and the reflected P and converted SV waves. φ S is related to the half-space S-wave velocity via Snell's law, so the measured P-wave polarization can be used to constrain this velocity. For the more realistic case of crustal layering instead of a half space, a frequency dependence of φ P , which could be used to constrain this layering, has been proposed by Phinney (1964). The S-wave velocities obtained at various periods are then apparent rather than true sub-surface velocities, with more and more information from larger depths contributing at increasingly longer periods. Those apparent S-wave velocities can be inverted for the true S-wave velocity structure of the subsurface (Svenningsen and Jacobsen 2007).
Measuring the apparent incidence angle directly in the raw recorded seismograms, e.g. by using hodographs, might be hampered by complex P-wave trains that include the influence of the source-time function and near-source wave propagation effects. Signals from nuclear explosions, which have comparatively simple source-time functions and well-known source properties, have been used in the study of P-wave polarization to avoid these complexities (Krüger 1994). More recently, Svenningsen and Jacobsen (2007) have suggested to use P-wave receiver functions instead of the raw seismograms to measure φ P . The deconvolution included in receiver function processing effectively removes source and distant wave propagation effects, so this approach allows using a larger amount of data and a more automated processing.
There are several potential benefits in applying this method to InSight data. Svenningsen and Jacobsen (2007) showed that the linearized inversion of apparent S-wave velocity curves calculated from P-wave polarization for crustal and upper-mantle S-wave velocities is independent of the starting model, in contrast to the direct inversion of the receiver function waveforms (Ammon et al. 1990). After adjusting the method to the boundary conditions at the ocean bottom, Hannemann et al. (2016) successfully applied it to an OBS data set with one to five receiver functions per station, thus demonstrating its usefulness when only a small number of event recordings is available. It has also been shown that a priori S-velocity information deduced from P-wave polarization measurements can be useful when attempting to model or invert the actual receiver function waveforms (Peng et al. 2012;Hannemann et al. 2017), even when polarization information at long periods is missing (Kieling et al. 2011). Schiffer et al. (2015Schiffer et al. ( , 2016 used the frequency-dependent apparent S-wave velocities as stabilizing information in a joint inversion with receiver function waveforms, whereas Chong et al. (2018) directly inverted the measured receiver function amplitude ratios instead. This, however, requires a sufficiently large amount of data with comparable ray parameters, which might not be available for InSight.
Here we calculate receiver functions and measure frequency-dependent apparent S-wave velocities from synthetic seismograms for a number of Mars models as well as for a complimentary data set of terrestrial data from Central European stations. We investigate the distance range likely usable for the calculation of P-wave receiver functions on Mars and demonstrate that the derived apparent S-velocity curves are distinct for the different models, even taking into account InSight's maximum event localization uncertainty. After comparing different forward computation schemes for Mars receiver functions to identify an accurate and fast method, we apply a grid-search based inversion approach to both synthetics and measured data and show that we are able to recover the main characteristics of the different models and stations well.

Mars Synthetics
The synthetic seismograms we employ here are generated using Green's functions (GF) databases that were produced for 16 different one-dimensional Mars models (Fig. 1). These GF databases are computed using the axi-symmetrical spectral elements code, AxiSEM (Nissen-Meyer et al. 2014), and are publicly available within Marsquake Service at ETH Zurich to allow any interested parties to generate on-demand synthetic waveforms , http://instaseis.ethz.ch/marssynthetics/). The GF databases can be accessed using the Instaseis package (van Driel et al. 2015) for computing synthetics within a duration of seconds, for arbitrary source-receiver pairs and moment tensors. The resulting waveforms are accurate down to a period of 1 s, and total simulation duration is ∼ 30 minutes; hence very suitable for body-wave studies. The study of  includes detailed information on generating on demand seismograms and available web services. The family of structural models used for the GF databases results from combinations of (1) thin (30 km, C30) and thick (80 km, C80) crust, (2) low (VL) and high (VH) crustal velocities, (3) two different mantle structures derived from the bulk composition of Mars by Sanloup et al. (1999) (SNL) and Taylor (2013) (T13), and (4) two different areotherms taken from Bertka and Fei (1997) (BF) and Khan et al. (2016) (AK). The areotherm profiles employed in structural inversions mainly influence lower mantle structure, i.e. velocity gradients below ∼ 1000 km where orthopyroxene and olivin phase transitions occur (Fig. 1). Additionally, for the Bertka and Fei (1997) areotherm, seismic Q is distinctly lower within the crust, and also slightly lower in the mantle (Fig. 2). The two different compositional models, on the other hand, mainly result in different velocities in the upper-most mantle, influencing the velocity contrast at the Moho (Fig. 2). The strength of this contrast, together Table 1 List of terrestrial stations from which receiver functions were used, including network code (EE-Estonian Seismic Network, GR-German Regional Seismic Network, PL-Polish Seismic Network) and time frame analyzed. Crustal thickness estimates are taken from Knapmeyer-Endrun and Krüger (2014). A map with station locations in their geological context can for example be found in Becker and Knapmeyer-Endrun (2018)  with the traveltime within the crust, is also the main difference between the "fast" and "slow" models. Furthermore, all models share the same core radius of 1589.5 km and core structure ( Fig. 1). More information on the inversion procedure for the models can be found in Khan et al. (2016) and references therein. In contrast to the mantle, these models do not consider mineralogy or temperature for the crustal structure. The thin and thick crusts with different velocity contrasts at the Moho represent 1-D global end-member models, rather than what is expected beneath the InSight landing site. Additionally, all the models consider a two-layer crust that varies in density and seismic velocities only, with an upper-crustal thickness of 10 km. Orbital topography and gravity data cannot resolve absolute crustal thickness uniquely and their interpretation depends on assumptions about the average crustal and upper mantle density and the minimum crustal thickness of Mars (Plesa et al. 2016). However, a commonly adopted model predicts an approximately 30 km thick crust at the landing site (Neumann et al. 2004), based on an crustal density of 2900 kg/m 3 . This thickness could increase to 40-50 km if the average crustal density is larger either globally or locally, i.e. if there exists a density dichotomy between the southern highlands and the northern lowlands of Mars (Plesa et al. 2016).
For the purpose of this study, we calculated synthetic seismograms using the aforementioned GF databases and Instaseis for events at epicentral distances between 15 • and 180 • with an increment of 1 • . We use a dip-slip source with an angle of 45 • , assuming normal faulting as a typical source type for Mars (Knapmeyer et al. 2006). Each event has a depth of 5 km and is located due north of the InSight station (4.5 • N, 136.0 • E), resulting in a natural separation of radial and transverse motion on the N and E component, respectively.

Terrestrial Data Set
To allow for a wider range of possible crustal structures and to study the effect of measured vs. noise-free synthetic data, we also analyzed data from five stations in Central Europe with known differences in crustal structure, including thickness, velocities, and presence of sediments. The considered stations are listed in Table 1, with information on crustal thickness derived by applying the stacking method of Zhu and Kanamori (2000) to P-receiver functions (Knapmeyer-Endrun and Krüger 2014), and their tectonic setting. Station BFO is located on the thinned crust of the Upper Rhine Graben, a part of the European Cenozoic Rift System, whereas station GEC2 sits on the thicker crystalline Variscan crust of the Bohemian Massif. Station BSEG is situated in a cave within an outcropping Permian salt diapir in the sediments of the North German Basin (Bormann et al. 1997). In contrast to the three previous stations that are located in Phanerozoic western Europe, both stations SUW and VSU are situated on the East European Craton. The craton is characterized by a three-layer instead of a two-layer crust, with an additional fast, mafic lower crustal layer and a more gradational transition to the mantle (e.g. Grad et al. 2003). It is also covered by a thin layer of slow sediments, which leads to strong reverberations in P-receiver functions for station SUW (Knapmeyer-Endrun and Krüger 2014;Wilde-Piórko et al. 2017, compare Fig. 7(b)). The receiver function data set used here is the same as in Knapmeyer-Endrun and Krüger (2014), covering 4.5 to 6 years of data. Due to the variable station quality, this resulted in similar numbers of individual usable receiver functions per station (Table 1). Data from station BFO have already been used as an example by Svenningsen and Jacobsen (2007) in their study of lithospheric S-wave velocities derived from P-wave polarization measured by receiver functions.

Receiver Function Calculation
For the terrestrial data set, P-receiver functions are calculated by first restituting the instrument response and filtering between 5 Hz and 50 s. Then, the horizontal seismogram components are rotated into radial and transverse direction, using backazimuths determined by polarization analysis (Jurkevics 1988). The synthetic seismograms do not require the removal of any instrument response, but they are filtered between 1 Hz and 50 s, 1 Hz being the upper frequency limit of the synthetics. Additionally, due to the alignment of source and receiver, these data are already in the ZRT system. The influence of uncertainties in the measured azimuth on Mars on the results is discussed below. For both data sets, the vertical component is deconvolved from the vertical and radial components in the time domain using a Wiener filter to obtain the receiver functions, as discussed in detail by Hannemann et al. (2017).

Apparent S-Wave Velocities
Apparent P-wave incidence angles φ P are estimated from the amplitudes of vertical and radial receiver functions, designated ZRF and RRF, respectively, at time t = 0, using the relation Applying Snell's law for the reflected SV-wave and the dependence between the apparent P-wave incidence angle and the true incidence angle of the SV-wave φ P = 2φ S results in the following equation linking the measured quantity φ P and the apparent S-wave velocities v S,app : with p being the slowness. Whereas the slowness for a given source and receiver location is well known for terrestrial data, the influence of the uncertainty in slowness on Mars, due to both location uncertainty and uncertainty in the velocity model, on the derived v S,app (T ) values is discussed in more detail below. A set of second-order zero-phase Butterworth low-pass filters is applied to the receiver functions to obtain the variation of apparent incidence angles with period (Hannemann et al. 2016). Following Svenningsen and Jacobsen (2007), the corner periods of the filters T f are selected to be logarithmically distributed. When directly comparing results with Svenningsen and Jacobsen (2007), it has to be taken into account that, for an equivalent filter band, the corner period of the Butterworth filters used here is twice the corner period of the cosine filters used by Svenningsen and Jacobsen (2007). To avoid measurements at periods shorter than the corner period of each event's source spectrum and to correct for any effect of the employed filter on the P spike of the receiver function, we measure the dominant period T rf of this spike on the vertical component receiver function for each trace. We then discard filter periods smaller than T rf and correct for the filter effect by adjusting the actual corner period of the filter T f to T 2 f − T 2 rf if this adjustment is larger than 1% of T f , as proposed by Hannemann et al. (2016).
To prevent contamination by noisy data, especially at long periods, we consider only measurements with a signal-to-noise ratio larger than 5 on both the vertical and the radial component at each separate period. The signal-to-noise ratio is estimated as ratio between the mean squared amplitudes of the filtered traces over time windows from 40 s to 25 s before the P onset for the noise and within ±10 s round the P onset for the signal. We calculate a median curve from these measurements for all periods for which at least 10 individual data points are available.

Inversion
Two main inversion approaches have been used in the literature: Either an iterative linearized least-squares method (Svenningsen and Jacobsen 2007;Schiffer et al. 2015Schiffer et al. , 2016, or a grid search over the parameter space (Kieling et al. 2011;Hannemann et al. 2016). Applications of linearized least-squares methods use a fast forward computation approach for the impulse response of a layer stack, which is applied to a number of starting models to account for a possible dependence of the inversion results on them. In the case of Hannemann et al. (2016), the grid search was implemented in a step-wise manner, simultaneously varying at most two parameters, due to the large computational demands when calculating the full teleseismic wavefield, even in a 1D structure. Furthermore, Chong et al. (2018) used a different approach by applying a simulated annealing algorithm to invert receiver function amplitude ratios.
As the crustal models in the Mars synthetics span a large range, and prior information is limited, we prefer the exhaustive grid search approach, in order to also characterize uncertainties and possible trade-offs. After identifying a suitably fast and accurate forward calculation algorithm (Sect. 4.3.1), we extend the inversion approach by Hannemann et al. (2016) to a complete grid search as performed by Kieling et al. (2011). Specifically, for the Mars models as well as for stations BFO and GEC2, we conduct a search over two layers over a half space, and for stations BSEG, SUW and VSU over three layers over a halfspace, including a sediment layer. The parameters of the grids are given in Table 2. For the Mars models, the grid is wide enough to include all variations within the data base down to 63 km depth, whereas for the terrestrial data, the selection of the grid parameters was guided by comparison with the v S,app (T ) curve for a standard Earth model (see Sect. 4.2.2). Resolution of the v S,app (T ) curves decreases with depth; for example, a velocity change between 1 and  Note the different range of the X-axis in both plots 2 km depth in the reference model IASP91 (Kennett et al. 1995) mainly influences periods above 4 s, whereas a velocity change between 15 and 16 km depth influences a significantly larger period range between 4 and 30 s, and a change between 21 and 22 km depth influences periods from 20 s to more than the 108 s, the maximum period considered in the calculation of sensitivity kernels (Fig. 3). In accordance with this observation, the spacing of the inversion grids is coarser for the lower crust and mantle. The grid spacing could be increased in a second step, though, focusing on the best solution obtained in this first step to study resolution in more detail, if desired. The maximum crustal thickness investigated for the Mars models is less than the 80 km Moho depth of the models with a thick crust. However, sensitivity kernels indicate that sensitivity to velocity changes at depth below 60 km is small and restricted to periods below 70 s ( Fig. 3(b)), which are covered barely, if at all, in the measured v S,app curves for these models (Fig. 8(c) and (f)). The period range of highest sensitivity in the kernels corresponds to the time range containing the multiples from a velocity change at that specific depth (Fig. 3). The polarity of the kernels is different for the mantle range due to a complex interplay between phase amplitudes and timings when changing the velocity in a thin layer in the mantle. In the grid search, we employ the additional constraint that velocity has to increase with depth. The question if the number of crustal layers can be determined from the v S,app (T ) curve alone is addressed in the discussion. For the Instaseis Mars synthetics, we invert for the median curve derived from receiver functions calculated at 12 epicentral distances, from 40 • to 95 • , in 5 • intervals, as no single curve proved to be a perfect match for the median curves shown in Fig. 8, and averaging over a limited set of curves from several distances is a more realistic scenario for InSight. For the terrestrial data, we likewise selected the median over 6-7 curves that provides a close representation of the over-all median curve in each case to limit the amount of synthetic receiver functions that needs to be calculated for each model. Again, no single curve covers the whole period range while simultaneously being a close match to the average curve at all periods. We evaluate the results of the forward calculation on the parameter grid in terms of a misfit function. For each calculated curve v mod with N period samples, the misfit to the observed curve v obs is calculated as

Usable Distance Range for Receiver Functions on Mars
For P-receiver function studies on Earth, the commonly used epicentral distance range is between 30 • and 90-100 • , to meet the requirement of steeply incident P-waves and avoid secondary arrivals, e.g. from mantle triplications or PP, within the time window of interest (Burdick and Langston 1977). Very few studies use a wider range, including a few events at smaller distances as well as the weak Pdiff phase (e.g. Chevrot and Girardin 2000;Levin and Park 2000). At larger epicentral distances, PKP has been used, especially to image dipping structures (e.g. Levin and Park 2000;Endrun et al. 2005;Knapmeyer-Endrun et al. 2013). The use of PP to increase the coverage in ray parameter and backazimuth as done by Gurrola et al. (1994) is rather rare, likely due to the higher noise level and lower frequency content of this phase, and possible overlap with other phases like PKP. On Mars, phase transitions in the mantle will occur at significantly larger depths than on Earth, if at all, due to the reduced pressure and temperature, meaning triplicated phases will influence seismograms at greater distances than on Earth. Besides, the smaller size of the planet results in a smaller travel time difference between P and PP, and the likely absence of an inner core has a distinct influence on core-related phases. Thus, it is not self-evident that the same phases and distance ranges as on the Earth can be used for receiver function calculation on Mars.
Exemplary synthetic seismogram sections (Fig. 4) show that P and PP arrive within less than 25 s from another for a source at 5 km depth and epicentral distances smaller than 45 • , Fig. 4 Radial component synthetic seismograms for a strike-slip source in 5 km depth and models (a) C30VH_AKSNL, with high velocities in the crust and a small velocity contrast at the Moho, and (b) C30VL_AKSNL, with a low-velocity crust and a large velocity contrast at the Moho. Travel times for P, PP and PKP phases are indicated as calculated by TTBOX (Knapmeyer 2004) whereas the travel time difference is still only about 45 s at 60 • distance. No strong triplicated phases are expected from the orthopyroxene phase transitions around 800 km depth, as the velocity increase, based on thermodynamic modeling of the phase equilibria (Khan and Connolly 2008;Rivoldini et al. 2011), is rather gradational and small, extending over about 25 km. The velocity increase due to the olivine-wadsleyite phase transition around 1000 km depth extends over a larger range, but is accompanied by a larger velocity increase, and an additional phase is predicted by the ray theoretical travel time curves between 64 • and 73 • for P and 129 • and 146 • for PP. It does not show up distinctly in the radial component seismograms, though (Fig. 4). Similar to Earth, the core shadow for P starts at about 105 • , with some Pdiff energy propagating to larger distances. The PKP phase shows a characteristic pattern caused by the absence of a solid inner core. The single PKP branch arriving between 146 • and 168-170 • is bent sharply when exiting the liquid core and emerges on the opposite hemisphere. Strong focusing of energy at the antipode is also apparent. Note that, contrary to Earth, PKP is a prominent first arrival only beyond about 156 • distance. The difference between models with high and low crustal velocities is clearly visible in the body wave trains of the synthetic seismograms. As noted by , a large velocity contrast at the Moho, present in the models with lower crustal velocities, results in long codas due to energy trapped in the crust. The coda after the direct onsets of P, PP and PKP is reduced significantly for the models with high crustal velocities, for which the velocity contrast across the Moho itself is rather minor (compare Figs. 2 and 4).
Receiver function sections for two representative end-member models, C30VH_BFT13, with the least distinct Moho phases, and C80VL_AKSNL, with the largest number of crustal phases, are shown in Figs. 5 and 6. Note that the time axes are different as converted phases arrive significantly later for the model with a thick, slow crust. We calculated receiver func- mark the direct converted Ps phase, and, with increasing travel time, the PpPs multiple and the PsPs+PpSS multiple from the intracrustal discontinuity, respectively. Vertical orange lines separate the distance ranges in which P, PP, and PKP phases are used for receiver function calculation tions over the complete distance range from 15 • to 180 • , based on different P-wave phases: direct P up to a distance of 120 • in the case of C30VH_BFT13 and 107 • in the case of C80VL_AKSNL, followed by PP, clearly recognizable by the lower frequency content of the corresponding receiver functions, to a distance of 156 • in both cases, and PKP afterwards.
In both cases, the direct converted phases (Ps) and the two multiples (PpPs and PsPs+PpSs) from the intracrustal discontinuity at 10 km depth arrive within the first 10 s and are clearly visible. The amplitude of the direct phases is strongly reduced between 76 • and the start of the PP receiver functions, whereas the multiples are only clearly visible for distances larger than about 30 • in the case of C30VH_BFT13 (Fig. 5). The Moho phase and its multiples are not recognizable at all in the data for this model, which is not surprising as the velocity contrast at the Moho is rather small, −171 m/s in v P and −80 m/s in v S . In contrast, the receiver functions for model C80VL_AKSNL show many additional phases, some of which also show up on the vertical component (Fig. 6). The direct Moho conversion appears at 12-15 s, with decreasing move-out with epicentral distance, whereas the multiples, with the opposite move-out direction, can be identified near 45 and 60 s, respectively. The timing of additional phases that occur between the direct Moho conversion and its multiples and show a move-out comparable to the multiples aligns with a P-wave reflection at the Moho (between 27 and 33 s on both components), the same phase converted to S at Both examples show some common characteristics in the evolution of the receiver functions with epicentral distance that are also observed for other models. At distances less than 35 • , there is a strong interference with later phases and not all multiples are readily apparent. Around 70 • and 140 • , the signal complexity at the onset at zero time increases due to triplications. This leads to a somewhat increased noise level, but the receiver functions still seem usable. Starting from somewhere between 105 • and 120 • , PP becomes the phase to use for receiver function calculation. This results in a significantly lower frequency content, with an average dominant period of the zero-time impulse on the ZRF function of 4.0 s compared to 2.2 s when using the direct P-phase. This means that the intracrustal conversion is no longer separable from the direct wave at zero time for high-velocity models (Fig. 5). At distances beyond 156 • , there is an overlap between PP and different PKP branches. The branch traveling through the opposite hemisphere is clearly recognizable by the inverted sign of the converted phases around 165 • . However, these are the only clear PKP Fig. 7 Stacked RRF sum traces for (a) the Mars synthetics for six different models, using epicentral distances from 35 • to 95 • at 1 • intervals, and (b) for the terrestrial test data set, using the P-wave range. RRFs are normalized by the direct P-wave peak at relative time zero on the corresponding ZRFs and not move-out corrected before summing. Traces are labeled by model in (a) and by station code in (b), and clearly identifiable Moho Ps-conversions are marked by gray triangles in the corresponding traces receiver functions and the distance range in which they might be observable is rather small, dependent on the mantle model, and will also vary for different core sizes. In the cases examined here, PKP only generates the first, dominant arrival for distances beyond 156 • . We conclude from this that PKP is likely not a useful phase for receiver function analysis on Mars, whereas P and PP might be used between 35 • and the core shadow, and within the core shadow, respectively.

Mars Synthetics
We concentrate on the direct P-phase for the measurement of apparent incidence angles. We use the epicentral distance range from 35 • to 95 • , where the v S,app (T ) curves derived for each distance are very similar for each model, as outlined by the data densities in Fig. 8. We limit the investigation to six models that cover the range of crustal and upper mantle structures contained within the ETH model pool. Receiver function sum traces for these models are shown in Fig. 7(a), and v S,app (T ) curves are compared in Fig. 8. Moho phases are clearly identifiable only in the sum traces of the slow models, with a significant velocity contrast across the Moho, whereas the RRF amplitudes at zero time are distinctly larger for the fast models. The difference in upper crustal velocities between the fast and slow model families is readily apparent in the apparent S-wave velocities at periods shorter than 5 s. Apparent S-wave velocities at the shortest periods covered, around 2 s, agree very well with the actual S-wave velocities in the shallowest layer of 2.0 km/s and 2.75 km/s, respectively. The increase to higher apparent S-wave velocities occurs at shorter periods, around 5 s, for the models with high crustal velocities ( Fig. 8(a)-(c)) compared to those with low crustal velocities, where apparent S-wave velocities only increase from about 12 s onwards ( Fig. 8(d)-(f)). For the models with a thin, high-velocity crust, the apparent S-wave velocities at long periods, below about 40 s, converge to the upper mantle S-wave velocity of 4.116 km/s (model C30VH_AKSNL, Fig. 8(a)) and 3.636 km/s (model C30VH_BFT13, Fig. 8(b)), respectively. For the models with a thin, slow crust, the mantle velocities are not quite reached within the period range covered (Fig. 8(d), (e)), whereas for the models with an 80 km thick crust, no information on the mantle is contained in the v S,app (T ) curves. They converge to the velocity of the lower crustal layer at long periods ( Fig. 8(c), (f)). All considered, the curve for each of the six models is clearly distinct and provides information on the crustal, and, in some cases, upper mantle velocity structure.

Terrestrial Data
As for the Mars synthetics, we focus on data from the teleseismic P-wave distance range. Receiver function sum traces for the five stations are depicted in Fig. 7(b). The uncertainty in the v S,app (T ) curves for actual measured data is somewhat larger than for the synthetics assuming perfectly known locations, especially for the more noisy stations located on thin sediments (Fig. 9). For real data, apparent P-wave incidence angles can be measured up to somewhat higher frequencies than for the synthetics that are band-limited at 1 Hz. The increase in frequency content is largest for the stations located on sediments, with minimum dominant periods of the zero-time peak on the ZRF of 0.6 s for stations VSU and SUW, and 0.9 s for stations BFO, BSEG, GEC2. Average dominant periods increase from 1.15 s at SUW via 1.4 s at VSU, 1.6 at BSEG and 1.8 at GEC2 to 2.0 at BFO. All of these values are lower than the 2.2 s average dominant period in the synthetics.
V S,app (T ) curves for all stations show clear deviations from the one obtained from synthetics for the global Earth model IASP91 (Kennett et al. 1995, see below for calculation of the synthetics), and are distinguishable from one another. For both stations BFO and GEC2, the curves agree with the IASP91 model for the shortest periods, but show higher apparent velocities for longer periods ( Fig. 9(a), (c)). The period at which apparent velocities larger then 4 km/s are first obtained is 12 s at station BFO, but 22 s at station GEC2, indicating a thicker crust at the later station. At stations BSEG, SUW and VSU, apparent velocities at short periods are significantly lower than derived for IASP91, caused by the sedimentary cover at these stations. The lowest apparent velocities are measured for station SUW ( Fig. 9(d)), whereas velocities lower than those of the reference extend to the longest periods, 7.5 s, for station BSEG, indicating a thicker stack of sedimentary rocks in the North German Basin ( Fig. 9(b)). At stations VSU and SUW, apparent velocities are consistently higher than expected for IASP91 at periods below 1.2-2.5 s, indicating a different crustal structure beneath the craton ( Fig. 9(d), (e)). Apparent velocities above 4 km/s are only obtained at periods of 30 s and longer, corresponding to a thicker crust.

Forward Calculation of v S,app (T ) Curves
To set up an inversion scheme for the data, an accurate and fast way to perform forward calculations of receiver functions is needed. The selected method not only has to provide the arrival times and amplitudes of the direct Ps conversions and their multiples, but also account for additional complexity, e.g. the Pp-phases visible in the Instaseis Mars receiver functions, or deconvolution effects, that will contribute to the receiver function amplitudes averaged over long time windows, i.e. filtered at long periods, and thus to the measured apparent incidence angles. To make sure we choose a forward calculation algorithm suitable for the task, we compare results for different algorithms used in the literature to the results of the Instaseis full waveform synthetics, which have been defined as a benchmark for In-Sight . As the Instaseis synthetics rely on AxiSEM databases, which take several thousands of CPU hours to calculate on a supercomputer with a frequency content of up to 1 Hz for Mars models (van Driel et al. 2015, Fig. 15), this method is not a feasible option for an inversion scheme that requires calculations for tenthousands of models.
To avoid additional issues that might be caused by adjusting algorithms designed for terrestrial seismogram or receiver function calculation to Mars, we start the comparison with a terrestrial model, IASP91 (Kennett et al. 1995), for which Instaseis synthetics are available via the Syngine web service at IRIS (Incorporate Research Institutions for Seismology, . In contrast to the Mars synthetics, these seismograms are only complete to 2 s. We calculated receiver functions from synthetic Instaseis seismograms at 30 • to 90 • epicentral distance at 1 • intervals for a 10 km deep dip-slip source, using the same processing as described above for the Mars synthetics. The resulting waveforms as well as the v S,app (T ) curves are compared with results from alternative methods, specifically the full wave field propagator matrix method QSEIS (Wang 1999), a 2.5D pseudo-spectral method (Endrun et al. 2005), and the reflectivity P-SV impulse response of a layer stack as implemented by Shibutani et al. (1996).
QSEIS can compute complete synthetic seismograms for a stratified 1D-Earth model and allows to separately compute the wave propagation from the source up to a specified depth on the receiver side and the near-receiver propagation, which increases efficiency when doing calculations for many different crustal models. Calculations are performed for the same source and receiver parameter as in the case of Instaseis, and processing to generate the receiver functions is the same. QSEIS has been used in a step-wise grid-search modeling of v S,app (T ) curves by Hannemann et al. (2016). Due to the relatively high computational demands, i.e. several days of CPU time per model even when only calculating near-receiver propagation, a more complete search in 5 dimensions was not possible in that study.
The 2.5D pseudo-spectral method (PSM) investigated here solves the wave equation in spherical coordinates using a Chebychev approximation. The teleseismic P wave coda is simulated by successively triggering impulsive explosion sources at the lower boundary of the model to generate a plane P wave front, with the timing of the triggering dependent on a predefined ray parameter. Synthetics are calculated here for the 1D model IASP91 for 13 different ray parameters of the incident wave field corresponding to epicentral distances of 30 • to 90 • in 5 • increments. Data from 60 receiver points are averaged before calculating the receiver functions using the same processing in terms of filtering and deconvolution as above. This method has previously been used to investigate 2D effects on receiver functions at both crustal and upper mantle depth (Endrun et al. 2005;Knapmeyer-Endrun et al. 2013), but not in terms of a formal inversion.
The simple forward calculation implemented by Shibutani et al. (1996) computes the impulse response of a layer stack in the P-SV system and directly outputs receiver functions, i.e. P-multiples on the vertical component are not taken into account. It has been used to invert receiver functions in the original implementation of the Neighborhood Algorithm (Sambridge 1999), and, related to apparent incidence angles, in the inversion of receiver Gray shading in the background indicates data density taking into account the complete set of Instaseis curves based on receiver functions from 1 • distance increments between 30 • and 90 • function amplitude ratios (Chong et al. 2018). It is the fastest method considered here, with a single calculation performed in a matter of seconds. However, Schiffer et al. (2015) caution that synthetic receiver functions based on the assumption of an impulsive event do not properly account for the observed waveform widths and complexity, and convolve their synthetic Z-and RRFs with observed ZRFs before calculating v S,app (T ) curves (Svenningsen and Jacobsen 2007).
In Fig. 10, we compare Z-and RRFs at an epicentral distance of 75 • as well as the v S,app (T ) curves averaged over 13 receiver functions calculated for epicentral distances between 30 • and 90 • at 5 • increments. The resulting curve for Instaseis is indistinguishable from the one averaging over 61 curves for a distance spacing in 1 • increments, and the reduced amount of data is both more realistic for InSight and easier to handle in calculations. Basic features of all receiver functions are the same. However, the ones assuming impulsive incident P-waves (PSM and simple modeling) miss the added complexity on the vertical component due to source complexity and deconvolution, and any deconvolution artefacts before zero time on the radial component. The corresponding v S,app (T ) curves deviate from the median Instaseis curve early on and lie outside of the spread of the Instaseis curves already around 10 s. The simple layer response calculations tend to overestimate apparent velocities, while the PSM calculation severely underestimates them at long periods. The QSEIS curve is close to the Instaseis one and lies within the spread of curves across the whole frequency range, though deviations increase below 40 s. Remarkably, the curve Gray shading in the background indicates data density taking into account the complete set of Instaseis curves based on receiver functions from 1 • distance increments between 35 • and 95 • ; color range is the same as in Fig. 10 derived from the simple forward calculations convolved with the ZRFs from Instaseis is virtually identical to the Instaseis curve.
Next, we extend the comparison to Mars models. We did not adapt the PSM code to Martian parameters (e.g. different planetary radius), as the above comparison indicates that it does not provide a good approximation of the Instaseis results, even for Earth. Comparisons for both fast (Fig. 11) and slow models (Fig. 12) indicate that, while using results from full wave propagation modeling is possible, the simple layer-stack forward calculation is also a viable and fast option, as long as the results are convolved with the measured ZRFs to approximate source complexity, additional phases, and deconvolution effects. Note that receiver function waveforms are more complex for the Mars examples than for IASP91, including later phases like PP in the case of C30VH_AKSNL, and a more complex set of reverberations and conversions in case of C30VL_AKSNL. For another of the slow models considered, C30VL_BFT13, the QSEIS results are closer to the Instaseis curve than any other curve, but the results from layer-stack modeling convolved with the Instaseis ZRFs still lie within the spread of the data, which is rather large at long periods (compare Fig. 8(e)).
In conclusion, we select the layer-stack impulse response, calculated with the code of Shibutani et al. (1996), and convolved with the measured ZRFs, for the forward modeling of the v S,app (T ) curves. This method has the benefits of being both good at reproducing the benchmark results and fast, with each forward calculation taking only a few seconds.

Grid Search
The grid search results for the Mars synthetics are displayed in Fig. 13. The minimum misfits achieved lie between 0.034 and 0.070. All models that have a misfit up to 0.1 larger than the minimum are used to calculate the data density and the median model and curve in Fig. 13. Depending on the model, these are between 865 and 1822 models. The same threshold was used by Hannemann et al. (2016), and the spread of the resulting curves gives a good approximation of the spread of the measured data (Fig. 8). Note that the maximum period to which the measured curves are available is different for the different models, with longer periods missing especially for models C30VH_BFT13, C30VL_BFT13, and C80VL_AKSNL (Fig. 8).
The best inverted curves give a very good approximation of the measured data, though some of the sharper kinks in the curves are not reproduced. The median curves are also close to the inversion input, though their misfits are up to 50% larger than that of the best-fitting curve.
The velocity of the upper crustal layer is well characterized by the best-fitting model in each case, and the depth of this layer is recovered within ±2 km of the true value. Resolution decreases with depth and the range of acceptable models increases. Moho depths of the best models are within 10 km of the true value for those shallow models that show a strong velocity contrast at the Moho (Fig. 13(d) and (h)). For models with a weak contrast, velocities in the second layer seem to be tighter constrained than for those with a strong contrast, but average over the velocity values in the lower crust and the upper-most mantle ( Fig. 13(b) and (f)). The median models do a better job in matching lower-crustal velocity and Moho depth in these cases, and are also slightly closer to the actual Moho depth for the models with a strong contrast. The 80 km Moho depth is not contained in the searched parameter space as we do not expect to be able to resolve velocity changes at this large depth (compare Fig. 3), and the v S,app (T ) curves for these models do not reach upper-mantle velocities ( Fig. 13(i) and (k)). However, the best fitting models for both deep cases have Moho depths at the lower boundary of the parameter space.
The upper mantle velocity is not well resolved in any of the considered cases, with models within the given misfit range covering the whole parameter space and the best fitting model often approaching or touching the parameter space boundary. The median models, however, give a good indication of the velocity directly below the Moho, also in cases with a gradient (Fig. 13(f) and (h)).
Results for the terrestrial stations are displayed in Fig. 14. Again, the calculated curves show a close agreement to the measured data, with minimum misfits between 0.0526 and 0.0761, comparable to the values achieved for the synthetic Mars data. The one exception is station SUW, where the values at the shortest periods, above 2 s, are not well reproduced even by the best model, and the minimum misfit is about twice as large as for the other stations at 0.1303. The median curves cannot capture all of the details of the data in some cases, but still provide a good fit to the measured curves. Due to the larger spread in the measured data, we consider all models within 0.15 of the minimum misfit when calculating the data density and median curves and models. This range contains between 4384 and, in the case of station BSEG, 33066 models. This case is exceptional, though, with less than 7000 models considered at all other stations.
Layer thicknesses and velocities are similar between best-fit and median model in most cases. Exceptions are the intercrustal layering at stations BSEG and SUW. For station BSEG, the best model has almost identical velocities in the second and third layer, which might indicate that the data could also be explained by a two-layer model. The data density plot indicates that velocities between 7 and 20 km depth are not well constrained, though the distribution is tighter in the lower-crustal layer. This can also explain the rather large number of models with an acceptable misfit found for this station. Velocity distributions in the lower crust are tightly constrained for stations GEC2, SUW and VSU. For station SUW, the best model deviates from the median significantly with regard to the thickness of the sediment layer and the sub-sediment velocity. However, the larger discrepancy between the measured and modeled curves at short periods indicates that the full complexity of the sedimentary layer cannot be captured by the current model parameterization. Either an even finer grid, or, more likely, a velocity gradient in the sediments would be needed to achieve a larger misfit reduction. The problems with fitting this part of the curve also result in a data density maximum that is located away from the best fitting model at shallow depth. We did not pursue the case of SUW further by using more complex models, as a detailed inversion of the sedimentary structure beneath this station is not the main aim of this study.
For all stations, the difference between the best and the median model is 3 km or less for the depth of the inner-crustal discontinuity and 6 km or less for the depth of the Moho. The best-fitting model for station BFO also shows close agreement with the model derived by Svenningsen and Jacobsen (2007), using linearized inversion and the timing of the Moho peak in the receiver function waveform as additional constraint. The main difference is that they used data up to 0.2 s at station BFO, which we were not able to derive from our receiver functions, and thus were able to constrain a sediment layer of less than 1 km thickness that influences the curve above 0.8 s and is therefore not resolvable with our data. Those stations that provided data to the longest periods also show tight constraints on the upper mantle velocity (Fig. 14(b), (d) and (f)), whereas for station SUW, with the most limited extent of the v S,app (T ) curve towards long periods, the best model lies at the boundary of the parameter space and good models cover the whole range of possible mantle velocities (Fig. 14(h)). Moho depths for the best and the median model are 26 km and 28 km, respectively, for BFO, 30 km and 32 km for GEC2, 28 km and 28 km for BSEG, 43 km and 37 km for SUW, and 46 km and 43 km for VSU. With the exception of station GEC2, where our inversion suggests a shallower Moho depth, the values for the best models are within 3 km of those provided for reference in Table 1. The upper mantle velocity for the best and median models is 4.6 km/s for stations in Phanerozoic western Europe, identical to the value obtained for station BFO by Svenningsen and Jacobsen (2007), whereas the best models show higher velocities of 4.8 km/s for station VSU and 5.0 km/s for station SUW on the East European Craton. Due to missing long period measurements, these velocities are however not well constrained compared to those obtained for the western stations. Crustal velocities are close to those of IASP91 for stations BFO and GEC2 in both upper and lower crust for both best and median models. The S-wave velocity of the upper-most, 6 to 8 km thick layer at station BSEG is significantly lower at 2.2 to 2.4 km/s. Both upper and lower crust show faster velocities than those in IASP91 for the cratonic stations, at 3.45 to 3.75 km/s and 3.9 km/s, respectively.

Distance Range for P-Receiver Functions and Spatial Distribution of Martian Seismicity
The level of Martian seismicity is assumed to be considerably lower than on Earth due to the absence of plate tectonics. The spatial distribution of seismic events is as yet unknown and can only be estimated. In particular, Knapmeyer et al. (2006) mapped 8500 globallydistributed surface faults on Mars and determined their maximum age to simulate spatial distributions of seismicity. Plesa et al. (2018), on the other hand, predict the seismic activity and the spatial distribution of seismicity of Mars based on 3D numerical thermal evolution models considering stresses due to mantle cooling and convection. They use different assumptions in these modeling, in particular two different isotherms to define the depth of the seismogenic lithosphere, different seismic efficiencies, and different crustal thickness models. We compare estimates for the spatial distribution of Martian seismicity from both studies with the usable distance ranges for various P-phases in calculating receiver functions in Fig. 15. We only show results for two of the three crustal thickness models used by Plesa et al. (2018), as the results for the third model (Neumann et al. 2004) are broadly similar to the ones for the DC model (Fig. 15(c) and (e)). The crustal thickness models are based on current gravity and topography data and different assumptions on the bulk crustal density.
The DC model has a density dichotomy in the crust, with a crustal density of 2900 kg/m 3 in the southern highlands compared to 3100 kg/m 3 in the northern lowlands. Model HC (Fig. 15(d) and (f)) has a uniform, but relatively high crustal density of 3200 kg/m 3 . Considering all surface faults mapped by Knapmeyer et al. (2006), seismic activity can be expected at all distances from the InSight landing site, specifically also in the range of 35 • to 105 • -120 • which, based on the considered models, can be used for P-receiver functions ( Fig. 15(a)). Considering only faults on surfaces dated as younger than 600 Myr drastically reduces the spatial extent of faulting, with extensional faulting centered on the Tharsis region east of the landing site, and compressional faulting limited to several regions in the northern hemisphere. As Tharsis is close to the edge of the core shadow for P-waves in the considered models, this might limit the useful P-wave energy from events in this area. Of all scenarios considered, this one is the least favorable in terms of using P-wave receiver functions to measure apparent incidence angles.
For the seismicity distributions based on thermal evolution models, the annual seismic moment budget is highest for the HC model and the 1073 K isotherm. Seismicity is nearly uniformly distributed globally in this case, which is the most favorable one. The distribution is also roughly homogeneous for the DC model with this thermal constraint, but with a lower annual moment release. For models based on the 573 K isotherm, planetary cooling is the dominant factor in the thermal evolution models, which tends to concentrate seismicity in regions of thin crust (Plesa et al. 2018), i.e. in the northern hemisphere and the Hellas impact basin. These models also produce seismicity at all distance ranges from InSight, but favor certain azimuths, and have an altogether lower annual moment release than the models based on a thicker seismogenic lithosphere. The thermal evolution models do not consider stresses caused by lithospheric loading, which might add an enhanced level of seismicity in the Tharsis region.
With the exception of the model only considering the youngest faults as causing any seismic activity, and thus restricting seismicity to Tharsis and a few isolated regions in the northern hemisphere, all of the models predict significant activity within the P-wave distance range that could be used for measuring apparent incidence angles within InSight's nominal mission duration of one Martian year. We thus expect the method presented here to be applicable to the InSight data set.

Influence of Location Uncertainty on Apparent Velocity Estimates
The v S,app (T ) curves for the Mars synthetics were calculated assuming that the event backazimuth as well as the slowness of the P-wave phase at each epicentral distance is perfectly known. This will likely not be the case for InSight on Mars, though. The L1 mission requirements specify a minimum location accuracy of ±25% for the epicentral distance and ±20 • for the backazimuth ).
The slowness not only depends on the accuracy with which the epicentral distance is known, but also on how accurate the velocity model used to calculate the slowness at a given epicentral distance is. To investigate what the combined uncertainty of both location and velocity model means in terms of slowness, we used TTBOX (Knapmeyer 2004) to calculate the slowness of P, PP and PKP phases for epicentral distances between 0 • and 180 • , hypocentral depths between 0 km and 100 km in 5 km steps, and a variety of Mars velocity models (Fig. 16). To better capture the uncertainty in seismic velocities of Mars, we did not only consider the models of the ETH Instaseis database here, but also the four models by Mocquet et al. (1996) with molar fraction of iron in the Martian mantle between 10% and 40%, varied in 10% increments, the two models by Sohl and Spohn (1997), model M6 (Gudkova and Zharkov 2004) and M13 (Zharkov and Gudkova 2005). A maximum hypocentral depth of 100 km has also been assumed in generating seismicity catalogues for testing procedures for InSight , though recent thermal evolution modeling indicates that under certain conditions, the seismogenic depth of Mars might extent down to 450 km (Plesa et al. 2018). The range of slownesses covered by the various models and event depths for the same distance approaches 1 s/deg for both P and PP (Fig. 16). Additionally considering a ±25% uncertainty in distance, we estimate a combined slowness uncertainty of ±1 s/deg for the P phase. However, the location uncertainty might well be significantly smaller, as  obtained distance errors of less than 5% and backazimuth errors of less than 1.5 • for both terrestrial data and Mars synthetics at regional to teleseismic distances with the single-station location algorithm proposed for InSight.
To investigate the effect of the various uncertainties associated with the InSight Mars data on our calculated v S,app (T ) curves for the Mars synthetics, we first randomly assigned backazimuth deviations normally distributed with a maximum absolute value of 20 • to the seismograms. We rotated the data into the ZRT system using these new backazimuths and re-calculated the receiver functions. The same deviations were used for all Mars models to ensure comparability between the results (Fig. 17(b), (e)). Next, we added the uncertainty in Fig. 16 Slowness versus epicentral distance for phases P, PP, and PKP, source depths between 0 and 100 km in 5 km increments, and various Mars velocity models calculated using TTBOX (Knapmeyer 2004). In addition to the ETH Instaseis models (Fig. 1), we considered the models by Mocquet et al. (1996), Sohl and Spohn (1997), M6 by Gudkova and Zharkov (2004), and M13 by Zharkov and Gudkova (2005) slowness by randomly changing the slowness values for each event by a value drawn from a normal distribution scaled to a maximum absolute value of 1 s/deg. The perturbed values were than used to calculate the v S,app (T ) curves using Eq. (2) (Fig. 17(c), (f)).
Including backazimuth errors does not noticeably change the spread of the distribution of v S,app (T ) curves ( Fig. 17(b), (e)). The median curve is also barely affected, the only clear change is a decrease in values at short periods of less than 3 s. Using an erroneous backazimuth in the receiver function calculation will only affect the horizontal components, i.e. in this synthetic noise-free, isotropic, 1D case it will decrease the amplitude of phases on the radial component. Lower amplitudes on the radial component at zero time will decrease the estimated apparent incidence angle of the P-wave and thus also decrease the estimated apparent S-wave velocity at this period (Eq. (2)). The largest phase on the radial component is the P onset at zero time, so changes will be most pronounced here. At short periods, only this phase is sampled when taking the amplitude ratio of the low-passed seismograms at zero time, whereas at periods longer than about 3 s, part of the energy of the Ps conversion from the intracrustal discontinuity, which will decrease less in absolute terms for an erroneous backazimuth, is already included. For real data, amplitudes on the transverse component are unlikely to be zero across the whole time window considered here. While the effect of noise might average out when considering longer time windows and several events, anisotropy or dipping layers might lead to characteristic amplitude patterns on the transverse component, which could result in amplitude errors on the radial component if rotated with an incorrect backazimuth. Strong Moho topography can be expected along the dichotomy boundary (Neumann et al. 2004), but this is located about 10 • south of the InSight landing site and will not be sampled by the receiver functions. We thus consider the effects of dipping layers and anisotropy to likely be of second order for InSight receiver functions.
The introduction of erroneous slowness values visibly increases the spread of the distribution of v S,app (T ) curves, though the effect on the median curve is negligible (Fig. 17(c), (f)).

Fig. 17
Apparent S-wave velocity curves for two different Mars models, successively taking into account the effects of location uncertainty. Plots show data density and median curve using receiver functions from 1 • distance increments between 35 • and 95 • at periods with at least 10 data points (black lines), and median curve using a reduced number of data, i.e. receiver functions from 5 • distance increments between 40 • and 95 • (green lines). Additionally, the gray lines in (b), (c), (e) and (f) indicate the original curve obtained for perfectly known event locations, equivalent to the black lines in (a) and (d). Top row is for model C30VH_BFT13, bottom row is for model C30VL_AKSNL, and the color scale is the same in each plot. Curves are calculated assuming perfectly known event locations ((a) and (d)), taking into account a ±20 • uncertainty in backazimuth ((b) and (e)), and additionally taking into account a ±1 sec/deg uncertainty in slowness ((c) and (f)) They are still very close to the median curves obtained for the data set with perfectly known event location and velocity model; even the offset at short periods caused by the erroneous backazimuth values is somewhat reduced. It has to be considered that these averages are obtained from 61 individual receiver functions at a large number of distances, though, while the data set expected for InSight might be significantly smaller. Therefore, Fig. 17 also shows the median curves derived from only 12 events at distances of 40 • to 90 • at 5 • intervals. A homogeneous event distribution with distance is in agreement with the available predictions for the Martian seismicity distribution (Fig. 15). Here, deviations at long periods are somewhat larger when including uncertainty in both backazimuth and ray parameter, but the median curves are still well within the uncertainty range of the data.
In our consideration of location errors, we assume a random distribution of these errors, i.e. no systematic location bias. This seems to be reasonable, though, as  find that realistic 3D models of Mars have little influence on the timing of P-and S-wave Gray shading indicates data density for the best fitting models as shown in Fig. 13(b) and (d). Red circle indicates the true solution, blue cross the best model, and green x the median of the best fitting models. Black boxes outline the investigated parameter space. Color range is the same for all plots arrivals compared to 1D models, and that locations based on fitting these arrivals in spherical symmetric models will effectively be equivalent to locations that rely on more complex models.

Trade-Offs
To further investigate the capability to resolve different model parameters with the v S,app (T ) curves, we considered trade-off plots, using the parameters of the best fitting models. Examples for models C30VH_AKSNL (Fig. 13(b)) and C30VL_AKSNL (Fig. 13(d)) are shown in Fig. 18. It is apparent that some parameters are less well resolved than others; specifically, models that explain the data well can be found for all values of the upper mantle velocity in the case of C30VL_AKSNL, with little concentration in the parameter space, especially when also considering variations in the properties (velocity and thickness) of the lower crustal layer (Fig. 18(b)). Some parameters show a clear trade-off, e.g. thickness vs. depth of the second layer, where a higher velocity is compatible with a larger thickness and vice versa. A similar trade-off is visible between the thickness of the first layer and the velocity in the second. This might explain the strategy of Hannemann et al. (2016), who search for the velocity of the second layer last, after determining thicknesses and velocities of all other layers, in their step-wise grid-search of a model with two layers over a half-space.
Parameters of the shallow-most layer show the largest concentration of models in specific regions of the parameter space, i.e. appear to be best resolved. Other parameters show a larger scatter, but also a concentration in specific regions, e.g. the Moho depth z 2 . In the case of C30VH_AKSNL, there are two regions of higher data density, one around 27 km, close to the true crustal thickness of 30 km, and one around 53 km, close to the best-fitting model ( Fig. 18(a)).
Interestingly, the true model is not always located in the region of highest data density in the grid search results. In the case of C30VH_AKSNL (Fig. 18(a)), the median model provides a much better approximation of the true model than the best-fitting model, especially considering parameters of the deeper layers. However, this is not the case for the properties of the shallow-most layer in case of the slow model, C30VL_AKSNL. Observations are similar for the second model family with a thin crust, but different mantle properties (BFT13 models), whereas differences are less clear for the model family with a thick crust. For the terrestrial data set, median and best model are very close together for a number of cases (stations BFO, GEC2 and VSU), whereas the differences seem to reflect trade-offs between the middle and lower crustal layers for station BSEG and problems in achieving an acceptable model fit at short periods for station SUW. In summary, considering not only the best model, but also the median as a representation of the whole ensemble of models that can explain the data helps to obtain a better idea of the uncertainty in the results and might often give a closer approximation of the true model.

Tighter Constraints on the Moho Depth
As some trade-offs between model parameters became apparent, we considered the question of whether having tighter constraints on the Moho depth would help to better resolve other model parameters. For the terrestrial data, some prior information on the Moho depth is available, while for InSight, other studies, e.g. on near-by impacts (Daubar et al. 2018) or ambient noise auto-correlations (Becker and Knapmeyer-Endrun 2018), might also provide a crustal thickness range tighter than the 40 km used here. We filtered the best fitting models by assuming that the Moho depth is known within ±10 km of the true value for InSight and is within ±3 km of the value given in Table 1 for the terrestrial data. Exemplary results for four models are shown in Fig. 19. For the Mars model C30VH_AKSNL with a small velocity contrast at the Moho, constraining the Moho depth significantly improves the inversion results in providing a better fit not only to the true Moho depth, but also to the true lower crustal and upper mantle velocities for both best and median model ( Fig. 19(a) and (b)). However, for the other Mars models with a weak velocity contrast at the Moho improvements are less striking, except in terms of Moho depth. For the models with a strong velocity contrast at the Moho, the best fitting models remain the same as the Moho depth was already resolved to within ±10 km without any additional constraints. Changes in the median models are likewise minor.
For the terrestrial data set, we consider stations BFO and GEC2, as these stations showed the largest difference between the inversion results and the prior information on Moho depth, while avoiding the additional complexities of station SUW. Constraining the Moho depth for station BFO makes the median model follow the results of Svenningsen and Jacobsen (2007) more closely, who also employed additional information on the Moho depth in terms of the timing of the Ps conversion in the receiver function waveforms. The best fitting model stays the same, as it already recovered a Moho depth within ±3 km of the a priori value. For  Table 1 with an uncertainty of ±3 km. (g) and (h) are results for station GEC2, assuming the Moho depth listed in Table 1 with an uncertainty of ±3 km station GEC2, a deeper Moho within the prior constraints leads to an increased, unusually high velocity of 4.05 km/s in the lower crust for the best model, mirroring the observed trade-off between these two parameters (Fig. 18). Velocities in the median model are not affected. It provides a worse fit to the exact curvature of the measured curve around 10 s and 40 s, as was already the case without constraints on Moho depth, but still is well within the measurement uncertainty.
In summary, tighter constraints on the Moho depth do not seem to improve the inversion results for other parameters much in most cases. The exception could be models with a very weak velocity contrast across the Moho that is not well imaged in receiver functions. However, it is unclear which other seismic methods would be able to resolve this weak contrast and could thus provide better constraints to use in the inversion in that case.

Model Parameterization
When selecting the model parameterization for the grid search, we assume the correct number of layers for the Mars models, and a two-layer crust for the terrestrial data, with an Fig. 20 (a) Stacked RRF sum traces for two different terrestrial crustal velocity models, using synthetic events between 30 • and 90 • epicentral distance in 1 • intervals. (b) Corresponding velocity models. (c) v S,app (T ) curves calculated from receiver functions for the two models. Depicted is the combined data density as well as the median curves in blue and green, corresponding to the colors used in (a) and (b) additional sedimentary layer if appropriate. However, the number of layers is in principle also a variable in the inversion process , and it would be interesting to know if this information can be obtained directly from the data. Peng et al. (2012) suggest to look for maxima in the derivative dv S,app (T )/dT to identify the depth of layer boundaries based on some scaling relation, but for more complex models, these boundaries do not correctly reflect the number of layers in the models, and the derived depths show large deviations as the v S,app (T ) curves only constrain the average velocity over a certain depth range (Fig. 3).
In Fig. 20(c), we compare v S,app (T ) curves for synthetic receiver functions calculated with QSEIS for two Earth models, one with two crustal layers and one with just one crustal layer ( Fig. 20(b)). Curves are averaged over 61 events at epicentral distances between 30 • and 90 • sampled at 1 • intervals. While the curves are not exactly identical at short as well as at long periods, the differences are so small that they would fall within the uncertainty of measured data (Fig. 9). However, crustal thickness, which is one of the main parameters we would like to constrain with the v S,app (T ) curves, is rather different between the two models, 24 km vs. 35 km. A more promising way to constrain the number of layers from the data is to also consider the actual receiver function waveforms (Fig. 20(a)), as for example done in the joint inversion approach by Schiffer et al. (2015) and Schiffer et al. (2016). The stacked traces clearly show an additional phase around 2 s, and the corresponding multiples at 7-9 s, for the first model, which are missing in the receiver functions for the second model.
In actual data, the number of crustal layers might not be that easy to discern, but RRFs for all stations considered here show more than one positive phase, aside from the direct wave peak at zero time, within the first 5-6 s ( Fig. 7(b)). For station BFO, the phases are not clearly separate, but the peak around 2-4 s is very broad, and for station SUW, multiples of the sediment layer completely overprint the first 8 s of the RRF waveform. Small sediment phases are also discernible for stations BSEG and VSU within the first ∼ 1 s. On the other hand, the RRFs for the fast Mars models do not contain any well-defined Moho phases ( Fig. 7(a)). To test the influence of the parameterization on the inversion results and see if, in fact, the data could also be fit by models with less free parameters, we conducted additional grid-searches for one layer over a half-space for seven data sets: C30VH_AKSNL, C30VL_AKSNL, C80VH_AKSNL, C80VL_AKSNL, BFO, GEC2, BSEG (Figs. 21 and 22). For station BSEG, we additional considered the inversion results for a model with two layers over a halfspace, since the intracrustal layer does not seem to be well resolved by the data (Fig. 14). The parameters of the single layer were selected broad enough so as to include the complete velocity and depth range covered by both crustal layers in the previous inversion. Results are displayed slightly differently for the Mars models here-as only a very small amount of models, between 22 and 31, fit the data within the assumed uncertainty (minimum misfit + 0.1), no data densities are calculated, but the individual models are plotted (Fig. 21). The median model is also of less significance here as in most cases, Fig. 22 Inversion results for a smaller number of layers for terrestrial stations. (a) Best fitting curves from inversion of v S,app (T ) curve for station BFO (minimum misfit + 0.15), assuming one layer over a halfspace. Gray shading indicates data density, red line is the measured curve, blue line the curve for the best model, and green line the median of all best fitting models. (b) Models belonging to curves in (a). Gray shading again indicates data density, blue line the model with the lowest misfit, green line the median of all models with a low misfit. Thin dashed lines in dark blue and green indicate best and median model of the original inversions with a parameterization of two or three layers over a halfspace (Fig. 14(b), (d), and (f)). Black lines give the outline of the sampled parameter space. (c) and (d)-same as (a) and (b) for station GEC2. (e) and (f)-same as (a) and (b) for station BSEG. Data density and median are calculated for models within 0.1 of the minimum misfit. (g) and (h)-same as (e) and (f), assuming two layers over a halfspace. Range of gray scale is chosen individually for each plot all models that fit the data are permutations of just a few distinct parameter values. For the Earth data, the number of models close to the minimum misfit (within 0.15 for BFO and GEC2, and within 0.1 for BSEG) is larger, between 78 and 374 models, so data density is used here (Fig. 22). In all cases, the misfit of the best-fitting model increases when only one crustal layer is inverted for. However, the increase is quite moderate, from 0.0618 to 0.0719, for C80VH_AKSNL, whereas the minimum misfit more than doubles for C30VH_AKSNL, C30VL_AKSNL and C80VL_AKSNL. A misfit increase of 21% and 37% is observed for stations GEC2 and BFO, respectively, whereas the misfit increase for station BSEG is 77% when using two crustal layers and more than a factor of 5 when using just one layer over a halfspace.
For all Mars models with a clear misfit increase, it is obvious that the shallow part of the model between 2.5 s and 5 s, in case of a 30 km crust ( Fig. 21(a) and (c)), and 3 s and 8 s, in case of a 80 km crust ( Fig. 21(g)), is not reproduced well by the inversion results. Misfits at longer periods (10 s to 30 s) are also large for the best fitting model in cases of a strong velocity contrast across the Moho (Fig. 21(c) and (g)). The depth of the single discontinuity for the best models is located roughly in the middle between the true depth of the intracrustal discontinuity and the Moho in these cases, with a shallower depth for the model with a thick crust, and the fit to the upper-crustal velocity is still quite close. This is similar to the synthetic experiment (Fig. 20), where both models shared the same uppercrustal velocity and the discontinuity in the one-layer model is located at 24 km, in the middle between the intra-crustal discontinuity (12 km) and Moho (35 km) of the two-layer model. For C80VH_AKSNL, the best-fit model is very close to the true model in terms of thickness of the first layer and velocities in the first and second layers. The small velocity increase at the Moho in 80 km depth is not resolved with the data, so that a comparable fit to the measured curve can be achieved with a one-or two-layer model. For the other three models, a significantly better fit is obtained when implementing the additional layer, though.
For the terrestrial data from stations BFO and GEC2, the difference in misfit is not that striking and there is no obvious large discrepancy between the best modeled curve and the measured data ( Fig. 22(a) and (c)). The single discontinuity is again located roughly in the middle between the two discontinuities of the two-layer models, whereas the estimated upper-crustal velocities remain similar ( Fig. 22(b) and (d)). For station BSEG, it is obvious that the single-layer model cannot explain the data both at periods shorter than 3 s and longer than 30 s (Fig. 22(e)), and the minimum misfit is very large in this case at 0.2968. Using a model with two crustal layers, the minimum misfit is still distinctly larger than for three layers, and the fit of the data is worse for periods below 20 s, although it still seems acceptable ( Fig. 22(g)). The velocity in the topmost layer is close to the one obtained for a three-layer model, and the thickness of the layer is identical (Fig. 22(h)). However, the velocity in the lower crust is considerably larger at 4.05 km/s, and the Moho depth of the best-fitting model is also more than 10 km larger than in the three-layer case, which is not deemed realistic.
For all considered terrestrial cases, a complimentary detailed analysis and interpretation of the receiver function waveforms could provide valuable constraints on the number of distinct sub-surface velocity contrasts. Additionally, the comparison of the terrestrial data to the theoretical curve for IASP91 provided some indication of deviations from this model and in turns implied that two to three layers are involved (Fig. 9), but for Mars such prior knowledge on crustal structure is missing. The v S,app (T ) curves calculated for the different Mars models in this study could offer a first baseline for comparison, though. Besides, the application of the corrected Akaike's information criterion (Sugiura 1978;Hurvich and Tsai 1989), which combines the misfit function with an expression penalizing model complexity and overfitting, could be used to identify the best model parameterization. The criterion has been previously applied in geophysical inversion of surface wave dispersion (Di Giulio et al. 2012) and Rayleigh wave ellipticity  for shallow subsurface structure. In the cases considered here, the criterion always provides lower values for, i.e. favors, the models with more parameters.

Long-Period Noise
Long-period noise might limit the ability to correctly measure apparent incidence angles at long periods, i.e. larger than several tens of seconds. A prominent cause of this kind of noise is tilt, which might also affect InSight's seismometer SEIS once it is placed on Mars. The noise model created for SEIS considers tilt noise caused by both thermoelastic tilt, i.e. due to differential thermal expansion or contraction of various parts of SEIS, and ground tilt due to atmospheric pressure fluctuations . While the pressure sensor carried by InSight allows for the decorrelation and removal of the latter, as demonstrated by , this is not the case for the former. Estimates for the vertical component indicate that the complete instrument noise, including the thermoelastic tilt, accounts for less than 20% of the total system noise at 0.1 Hz, whereas the pressure noise contributes more than 50% . At 0.01 Hz, the relative contribution of both noise types decreases, and the difference between them is less striking at approximately 10% and 20%. The relative contributions may be different for the horizontal components, though, where the noise is larger in general, especially during day time. As we select the used v S,app curves based on their signal-to-noise ratio on both vertical and radial component receiver functions at each period, we would probably exclude data contaminated by tilt noise. Though thermoelastic tilt does not seem to be a dominant source of noise, it will increase in relevance when other noise sources (e.g. pressure and magnetic noise) are corrected for. This might mean that less data will be available at periods long enough to constrain velocities at several tens of km depth for the actual SEIS recordings compared to the synthetic case considered here.
Another possible influence on the measured apparent incidence angles is exerted by longperiod phases arriving shortly after the direct P-wave, e.g. the W phase (Kanamori 1993). Hannemann et al. (2016) cite the W phase as a possible reason for a decrease in v S,app at long periods in their synthetic data. We would expect this phase to also show up in the full waveform synthetics we calculated for the Mars models if they might cause a problem there (i.e. those waveforms calculated by AxiSEM/Instaseis and QSEIS), but we do not observe any specific waveform features or influence on the v S,app curves (Figs. 11 and 12). Chong et al. (2016) show, for synthetics as well as real data, that the joint inversion of Rayleigh wave ellipticity and receiver functions might be insufficient to correctly and uniquely resolve average crustal velocities. Specifically, they present the pathological case of three synthetic models that can be derived from one another by scaling P-and S-wave velocities as well as layer thicknesses by an equal factor. The ellipticity curves for these models are indistinguishable, and Chong et al. (2016) show that the receiver functions are also quite similar and an inversion might converge to a wrong mean crustal velocity. To investigate whether v S,app (T ) curves can resolve the differences between these models, we compare v S,app (T ) curves from synthetic receiver functions for a single event at 75 • epicentral distance, calculated with QSEIS. The three curves clearly differ across the whole frequency range considered (Fig. 23), and even taking into account the spread in measured data, at least the slowest and fastest model should be readily distinguishable from actual measurements. Thus, including v S,app (T ) curves as additional data has the potential to improve an inversion based on receiver function waveforms and Rayleigh wave ellipticity curves.

Conclusions
Based on the Green's function data base at the Marsquake Service, we calculated receiver functions and v S,app (T ) curves for six Mars models. We find that these curves, without any Fig. 23 (a) Three velocity models scaled to each other as used by Chong et al. (2016) to illustrate the limitations of the joint inversion of Rayleigh wave ellipticity and receiver functions in constraining absolute crustal velocities. (b) v S,app (T ) curves for the three models, derived from QSEIS synthetics at 75 • epicentral distance additional constraints, already provide useful information on crustal thickness and structure, with especially good resolution of the S-wave velocities in the upper half of the crust. This remains the case when taking into account the maximum location uncertainties of InSight, and the likely limited amount of available seismograms. These results are further corroborated by the complementary analysis of data from five stations in Central Europe with a variety of crustal thicknesses and structures. In combination with results obtained from other methods, e.g. using P-wave travel-times of independently located impacts, that mainly constrain P-wave velocities, the retrieved S-wave information could be used to constrain crustal rheology.
We define the usable distance range of P and PP phases for receiver function calculation for the Mars models and find that, based on predictions of the seismicity distribution on Mars, it is very likely that a sufficient number of P-receiver functions to use for apparent incidence angle calculations will be available during InSight's nominal mission duration of one Martian year.
The large number of receiver functions pre-calculated in this study could be used as a starting point for the interpretation of actual measured Mars data. Interpretation of the v S,app (T ) curves is somewhat non-unique regarding the number of layers in a model, as they average velocity information over certain depth ranges. The selection of an appropriate number of layers for an inversion would therefore best use also additional information when available, e.g. from the detailed interpretation of receiver function waveforms.
The method is tested in a stand-alone fashion here, to see if it could be applicable to In-Sight data from Mars and provide useful information. However, it could well be combined in an inversion with other data that aim at resolving crustal structure, e.g. receiver function waveforms or Rayleigh wave ellipticity curves, especially as it can potentially resolve differences between models that are indistinguishable by these methods. and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.