MSS/1: Single‐Station and Single‐Event Marsquake Inversion

SEIS, the seismometer of the InSight mission, which landed on Mars on 26 November 2018, is monitoring the seismic activity of the planet. The goal of the Mars Structure Service (MSS) is to provide, as a mission product, the first average 1‐D velocity model of Mars from the recorded InSight data. Prior to the mission, methodologies have been developed and tested to allow the location of the seismic events and estimation of the radial structure, using surface waves and body waves arrival times, and receiver functions. The paper describes these validation tests and compares the performance of the different algorithms to constrain the velocity model below the InSight station and estimate the 1‐D average model over the great circle path between source and receiver. These tests were performed in the frame of a blind test, during which synthetic data were inverted. In order to propagate the data uncertainties on the output model distribution, Bayesian inversion techniques are mainly used. The limitations and strengths of the methods are assessed. The results show the potential of the MSS approach to retrieve the structure of the crust and underlying mantle. However, at this time, large quakes with clear surface waves have not yet been recorded by SEIS, which makes the estimation of the 1‐D average seismic velocity model challenging. Additional locatable events, especially at large epicentral distances, and development of new techniques to fully investigate the data, will ultimately provide more constraints on the crust and mantle of Mars.


Introduction
Because of its higher-resolving power relative to other geophysical methods for sounding the interior of a planetary body, seismology has played a prominent role in the study of the interiors of Earth and Moon. This is one of the primary reasons for landing a seismometer on Mars with the InSight mission (Banerdt et al., 2013). The InSight (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport) lander successfully delivered a geophysical instrument package to the Martian surface on 26 November 2018, including broadband and a short-period seismometer instrument package called SEIS (Seismic Experiment for Interior Structure) . Although not buried but deployed at the Martian surface under a Wind and Thermal Shield, SEIS is specifically designed to record marsquakes and meteoritic impacts under Martian conditions, with a very low noise level (see Lognonné et al., 2020, for noise levels, as compared to the prelaunch estimation provided by Mimoun et al., 2017).
Most of our knowledge about the internal structure of Mars has been inferred from geodesy data that has been supplemented with assumptions about the bulk composition of Mars. From the precise tracking of ©2020 The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. spacecraft orbiting Mars, the static gravity field and dynamic gravity field have been determined, and the tracking of surface landers allows for the determination of the planet's precession rate. By combining the static gravity field and the precession rate, the moment of inertia is determined from which the mass distribution within the planet can be constrained. Supplementing the gravity field data with topographic data allows constraining the structure of the crust, whereas the dynamic gravity field or tides allow constraining the core radius and the rigidity of the mantle. Constraints about the bulk composition of the planet are deduced from a large set of Martian meteorites, in situ rock analysis obtained from Martian rovers, surface spectroscopy, and assumptions about how the planer formed. All these constraints have allowed for several estimates of the internal mechanical properties and compositional structure of Mars (e.g., Baratoux et al., 2014;Khan et al., 2018;Mocquet et al., 2011;Neumann et al., 2004;Plesa et al., 2015;Rivoldini et al., 2011;Smrekar et al., 2019;Sohl & Spohn, 1997;Verhoeven et al., 2005). However, these models suffer from the intrinsic nonuniqueness of any gravity data inversion, in particular with respect to the structure and thickness of the crust, structure, and composition of the mantle, including the existence and depth of discontinuities, and the core size.
The InSight mission extends planetary seismology to Mars. Seismic data will be integrated with geophysical measurements to determine details of the internal structure and evolution of another terrestrial body for the first time. However, the determination of the internal velocity structure and the location of the seismic sources using only one station is a challenge. Routine operations in the InSight team are split into two services: the Mars Structure Service (MSS) and the Marsquake Service (MQS), which are responsible for defining structure models (Panning et al., 2017) and seismicity catalogs (Clinton et al., 2018), respectively. While the two tasks are intimately related and require constant feedback and interaction, these two services provide a structure to ensure the mission will meet its science goals.
Following early works (Khan et al., 2016;Panning et al., 2015Panning et al., , 2017, the MSS team has developed several other different inversion algorithms in order to retrieve the first 1-D-averaged models of Mars from single station seismic data. To deal with the large inescapable uncertainties associated with the quake parameters and 1-D structure model from a single station, we developed probabilistic inversion strategies with several types of seismic observables: body wave phase arrivals, surface waves, and receiver functions (RFs). In all cases, significant efforts have been devoted to the modeling approaches: when data are limited, they are indeed keys for understanding the significance of the resulting models.
The goal of this paper is to present the results of a blind test in terms of structure, in order to show how MSS can investigate Mars's interior using a variety of well-suited methods performed on a single high-quality seismogram. We took as "blind data" a synthetic seismic event and its associated broadband seismogram, computed within a 3-D crust overlaying a 1-D model from the Moho discontinuity to the core. All parameters used in the forward modeling were unknown to those in charge of inversions. This single synthetic data set was then used to perform inversions for inferring both Marsquake parameters (location, depth, origin time, and moment tensor) and interior velocity models.
This study provides a clear framework to describe the methods employed by the MSS demonstrated on a data set for which results can be compared to the ground truth. Clearly, this remains a "best case" scenario for the resolving power of a single event, but it is an important demonstration of the planned approaches. The different algorithms handle surface waves and/or body waves, RFs, and consider different parameterizations. We first describe the synthetic input data, the traveltime measurements, and then detail the different inversion techniques for investigating the deep structure of Mars, and quake location. Constraints from RFs, which provide information on crustal structure directly below the receiver, and from surface waves and body waves, which average the structure over the great circle path between source and receiver, are considered independently here. We first present, discuss, and compare the results of six different methods estimating the 1-D average model over the great circle path between source and receiver. Second, we described the results from a local study below the station, by inverting RFs. The results clearly highlight the nonuniqueness of the solution and enhance the approach considered by the MSS of using several complementary methods to fully investigate the data. We also demonstrate the feasibility to estimate the attenuation of Mars, using a single event. Lognonné et al., 2020). Only two marsquakes show clear P and S arrivals with picking errors smaller than 2 s for P and S, which means that several methods developed in the framework of the MSS cannot yet be used with SEIS data. The core of the proposed methods will, therefore, be applicable if Martian seismicity provides soon larger quakes with body wave phases and first orbit surface wave dispersion, or one event large enough to record multiple orbit surface waves. A section is therefore dedicated to the inversion results using only body wave arrival times. The results show a strong trade-off between the seismic velocity profile and the quake location, which prevents a clear estimation of the crust and mantle velocity structure, as well as the crustal thickness. In the absence of surface waves, the work of the MSS becomes even more challenging, and new methods need to be developed in order to fully take advantage of all the information contained in the SEIS data.

Synthetic Waveforms
Two blind tests were performed prior to landing to validate the SEIS data processing. The first one was for the detection and characterization of events and was linked to the MQS activities (Clinton et al., 2018). We refer to Clinton et al. (2017) and van Driel et al. (2019) for further details on the MQS tests results. The MSS blind test was on its side designed for practicing and improving the procedures and methods aiming to invert Mars's internal structure from SEIS data.
To mimic expected Mars conditions before InSight touchdown, the synthetic test data set included four Earth day long seismic recordings (three-component VBB data at 2 and 20 sps, single-component combined VBB vertical channel at 10 sps) as well as auxiliary channels such as pressure, magnetometer, wind speed, wind direction, and atmospheric temperature.
A single marsquake is hidden in the continuous signal, which is contaminated by Martian noise based on prelaunch hypotheses (Mimoun et al., 2017). The event's parameters are detailed in Table 1. The background structural model, or 1-D base model, for quake simulations was chosen from a suite of 14 models . In order to make this blind test challenging, an anomalous 1-D model is considered, with a temperature profile in the crust and mantle close to the liquidus (Khan et al., 2018). Such a temperature profile gives a seismic velocity profile located in the extreme lower bound of the expected velocity models of Mars (Smrekar et al., 2019). The source parameters were randomly chosen.
In the MQS exercise focused on event detection and location van Driel et al., 2019), the seismograms were computed in a 1-D model, which greatly simplified the analysis. The goal here is to include unexpected behavior as well as 3-D effects on surface waves due to crustal heterogeneities (Bozdaǧ et al., 2017). To this end we performed a full 3-D computation (Afanasiev et al., 2019) for S waves and surface waves to the shortest period of 5 s and merged the trace in the time domain just before the S wave arrival with 1-D synthetics (Bozdaǧ et al., 2017;Nissen-Meyer et al., 2014;van Driel et al., 2015) accurate to 1.5 Hz. This is motivated by the observation that the attenuation of the model would not allow shorter period teleseismic S waves above the noise level in any case. For consistency, the 1-D synthetics are aligned to the 3-D P wave arrival time by cross correlation of the 1-D and 3-D synthetics in the frequency range where both are well resolved. The seismic noise model comes from the SEIS noise model (Mimoun et al., 2017). This comprehensive model includes prelanding estimates of the key contributors to the seismic noise: the pressure noise (Kenda et al., 2017;Murdoch et al., 2017a), the wind-generated mechanical noise (Murdoch et al., 2017b), the thermal and magnetic noise (Mimoun et al., 2017), and the instrument self-noise . Further details about how this noise was integrated into the synthetic waveforms are provided in Clinton et al. (2017).

MQS Estimates
Since raw synthetic waveforms are supposed to mimic the worst-case scenario planned for data transmission of the continuous signal acquired on Mars, they are provided at a sampling rate of 2 Hz. ote. MsM is the magnitude derived for Mars using the surface wave amplitudes (Böse et al., 2017).
The probabilistic location algorithms that MQS utilizes are explained in Panning et al. (2015) and Böse et al. (2017) in detail. MQS uses a collection of approximately 2,500 Mars models, which is a combined set of inputs from all members of the InSight science team. The MQS algorithms operate on a traveltime database for the most common body wave phases at all distance ranges, as well as surface wave arrivals in predetermined frequency bands. Traveltimes are computed using the TauP package of Crotwell et al. (1999). Table 1 summarizes the computed location and the input parameters. Although direct and major arc surface wave arrivals were visible, it was not possible to identify a clear arrival of multiorbit surface waves (R3). Further, 3-D traveltime corrections for surface waves were not implemented at the time of the test. Therefore, the MSS location was computed solely using body wave phases. The event depth (∼36 km) was constrained using a clear pP arrival 14.4 s after the P wave. The event magnitude (MsM ¼ 3.7) is computed using the direct Rayleigh wave arrival. Following Böse et al. (2017), this value scales to approximately Mw ¼ 4.2-4.3.

Body Waves
Onset times are picked by various contributors, they are gathered by the institution names (IPGP, ISAE, ETH, MQS, and LPG), and the major seismic phases are represented in Figure 1. Concerning the compressional seismic phases (i.e., P, pP, and PP), there is a very good agreement between all manually performed picks, whereas we observe a slightly larger discrepancy for shear wave onset times.
The mean value for the arrival time of the direct P wave is 2019-01-03T15:09:54.5 (UTC) with a standard deviation lower than 0.1 s; the corresponding value for the arrival time of the direct S wave is 2019-01-03T15:18:31.1 (UTC) with a standard deviation of 3.6 s. The larger uncertainty for the direct S phase is mostly due to a low-amplitude, emergent direct shear wave arrival, which is moreover almost concomitant with a high-frequency signal, which has not been identified, and no clear S wave coda (Figure 2).
To determine back azimuth (in the horizontal plane) and arrival (in the vertical plane) angles of the direct P wave ray, we use an approach based on the analysis of the polarization of the incoming P wave energy. The seismogram is band-pass filtered between 0.2 and 0.5 Hz in order to enhance long-period particle motion around the hand-picked P wave onset time. The window of interest starts 1 s before the picked P wave arrival and ends 5 s after.
The analysis is based on the approach from Jurkevics (1988); that is, we solve an eigenvalue problem on the 3-D particle motion to find the eigenvector that is a good approximation of the P wave vector. Its orientation and polarity are then used to recover the back azimuth of the P wave ( Figure 3a). To strengthen the analysis, which is very sensitive to noise, we also perform a Monte Carlo exploration of the 3-D particle motion obtained from the 6 s analysis window. Two hundred sets of samples from between 60% and 90% of all the particle motion samples are randomly selected, and they are used to compute the average back azimuth. This process is repeated 200 times so that the preferred back azimuth value (Figure 3b) is given by the mean and the standard deviation of the distribution (given by a kernel density estimation). The inferred value for the back azimuth is in agreement with the source location.
The arrival angle is computed following the same strategy ( Figure 4). The back azimuth angle is used to rotate the seismogram in order to use a vertical/radial/transverse reference frame. We also limit the analysis to the radial/vertical particle motion as almost all of the P wave energy should be confined to this plane.

Surface Waves
Both Love waves on the transverse component and Rayleigh waves on the two other components (Figure 2) are easily recognizable in the seismograms. Most of the surface wave energy is between 10-and 60-s period. The rotation in the great circle plane, defined by the source and the receiver, efficiently isolates Love waves that travel faster than Rayleigh waves. Two types of observables can be derived from raw data: discrete group velocity values at different frequencies and dispersion analyses using probability density functions (pdfs) of the energy.
As for body waves, arrival times for Love and Rayleigh waves are manually picked and are converted into group velocities for the different wave trains. Figures 5a and 5b exhibit the group velocity curves, with the corresponding uncertainties, between 10-and 60-s period, obtained from the G1 and L1 wave trains (travel path along the minor arc between the source and the receiver). They are computed using the origin time and the quake location as given by the MQS (Table 1). The major arc wave trains (i.e., G2 and R2) are picked with good consistency between the various contributions, but since they are not used for inversions, they are not shown here. Only two contributors picked G3 and R3 arrival times in a very narrow frequency band (between 25-and 30-s period).
For both Love and Rayleigh waves, all contributions are very consistent over the whole frequency range. Almost all group velocity values are within the uncertainties defined by the ETH and the MQS. This feature is consistent with a mean model for which body wave velocities are increasing as a function of depth.

Dispersion Analysis
In addition to the arrival time picks for the different wave trains of surface waves, we present a dispersion analysis of the first fundamental mode waveforms (Figures 6a and 6b). Both Love and Rayleigh waves are

10.1029/2020EA001118
Earth and Space Science filtered using 20 narrow Gaussian band-pass filters, and each trace is converted into energy using the envelope of the signal. The signal is then interpolated on a regularly discretized group velocity axis, assuming an origin time value and an epicentral distance (the location from MQS is used, Table 1), and the whole window is finally converted into probabilities, which is a slightly different version of the one described in Panning et al. (2015). The collection of pdfs of group velocities displayed in Figure 6, therefore, can be seen as a dispersion diagram. Two histograms, corresponding to two different frequencies, are displayed above the dispersion diagrams to illustrate this approach.
When looking at the dispersion for periods greater than 12 s, it is clear that group velocities are unimodal, and the maximum of each pdf matches almost perfectly the picked values shown in Figure 5. At shorter periods, we observe a sharp discontinuity in both pdfs. They are highlighted by bimodal distributions (see red histograms in Figures 6a and 6b). This feature affects both Love and Rayleigh waves and prevents inversions for periods shorter than 10-12 s. We do not observe any robust correlation between Love and Rayleigh signals that could indicate the presence of anisotropy (Babuska & Cara, 1991), and we interpret this short-period energy in terms of 3-D propagation effects. Since we do observe a 2π phase shift between   longitudinal and vertical components (see Figure 7 for the fundamental mode), we are facing retrograde particle motion. This spurious signal can, therefore, be associated with multipathing since there is a large crustal thickness anomaly (up to 20-25 km) in the source neighborhood. Hence, this energy can propagate along the same travel path and therefore reach the receiver with no significant arrival angle anomaly (Woodhouse & Wong, 1986).

Particle Motion
As explained in Panning et al. (2015), the particle motion of the Rayleigh wave fundamental mode can be used for back azimuth determination but it can be used as well to investigate the physical properties of the crust. In a velocity model showing a substantive velocity increase as a function of depth, Rayleigh waves are characterized by a retrograde elliptical particle motion. In many cases on Earth, the particle motion of the Rayleigh wave fundamental mode is retrograde while higher modes propagate in prograde motion. Very few studies report prograde observations of the fundamental mode (Tanimoto & Rivera, 2005), and this effect due to propagation through a thick sedimentary basin. In the case of a retrograde motion, once the horizontal components are rotated into the source-receiver great circle reference frame, the longitudinal (in the direction of positive minor arc propagation) is π/2 phase-shifted with respect to the vertical. It means that horizontal displacements are in phase with the Hilbert transform of the vertical component multiplied by −1 (Gribler & Mikesell, 2019). Thus, such a comparison between longitudinal and vertical waveforms can bring some information on the near-surface properties. The result of the particle motion analysis is displayed in Figure 7. For 170 time samples (blue and red points), which corresponds approximately to the duration of the Rayleigh fundamental mode wave train, the normalized cross correlations between the longitudinal component and the Hilbert transform of the vertical multiplied by −1 (associated with blue color) and +1 (in red) are computed. Each color point represents the zero-lag value of these cross correlations. In order not to be affected by the length of the time window used for the computations, all possible lengths between 150 and 350 s are tested and the best zero-lag value is retained.
It is obvious that (i) all blue points are positives which is consistent with a retrograde particle motion and (ii) the largest value is almost equal to one and it corresponds to the beginning of the Rayleigh wave train (thick blue line). The top graph on the left (Figure 7) shows for this time window that the normalized cross For each frequency, the envelope of the filtered signal is converted into a probability density function of group velocities and are plotted in gray scales. Two histograms (red for short periods and green for long periods, as indicated by the colored lines in the dispersion diagrams) are shown above the dispersion diagrams. For Love waves and vertical Rayleigh waves, the red histograms are shown for period values of 10.7 and 12.9 s, respectively. The green histograms corresponds to a period of 30.8 s. correlation between the longitudinal and −HðZÞ (blue signal) is almost identical to the autocorrelation of the longitudinal component (black signal). On the other hand, the best value for the cross correlation between longitudinal and positive Hilbert transform is for the last point (top graph on the right), which means that the corresponding time window is after the end of the Rayleigh fundamental mode wave train.

Methods
Six independent inversion algorithms were developed in order to retrieve the 1-D average model along the minor arc. The goal here is not to produce a single model, but rather a family of models that fit the data and are consistent with the most recent set of prior constraints. The main characteristics of each method are described in Table 2. This study spans a relatively wide range in terms of model parameterization from the standard seismic parameterizations over fully self-consistent thermodynamic methods. Indeed, two main approaches are considered. One set of models (called M1) are parameterized in seismic velocity as a function of depth. A second set (called M2) is obtained with models parameterized by geodynamical parameters like temperature and composition. Assuming thermodynamic equilibrium, the seismic velocity profiles can then be calculated using first thermodynamics principles with experimentally derived parameters for candidate minerals. The strength of the M1 models is that they are able to mimic "exotic" models if the composition and/or temperature are variable along the wave path, or if the equilibrium assemblage is not reached. Their weakness is that they do not take into account constraints from mineral physics, nor geophysical data (e.g., moment of inertia, tidal response, and thermal evolution of the planet). The M2 method permits the application of tight constraints on velocity structure with a relatively limited data set, by producing stable velocity models through the whole planet, in contrast to the M1 models which give some constraints only at the depths where the data are sensitive. The M2 modeling approach is extremely powerful if we have very good prior constraints. However, the M2 models will not be representative of Mars if the prior assumptions on mineral physics turn out to be false or if the equilibrium assemblage is not reached. The advantage of using these two complementary approaches is that comparisons between these two families could indicate regions of the models that are not well constrained by the data, or inconsistent with the physical assumption of the M2 modeling approach.
The inverse problem consists in retrieving the seismic velocity profiles as a function of depth from body waves and surface waves traveltimes measurements. The problem is, however, underdetermined, in the sense that several combinations of velocity profiles and locations can give similar arrival times. Due to the underdetermined nature of the problem, we mainly use Bayesian inversion techniques to obtain robust pdfs of seismic velocity profiles. This technique allows the investigation of a large range of possible models and provides a quantitative measure of the models' uncertainty and nonuniqueness. As such, it is well suited to our problem given the still poorly known nature of the Martian interior, as well as the low amount of identified phase arrivals recorded by SEIS at this time Lognonné et al., 2020). The forward problems are computed with 1-D structure, which allows the computations of several million forward problems. The pdf of the a priori distribution, for each of the six methods, is displayed in Figure 8. In the following, the six methods are described in details. The M1 methods were initially developed for Earth applications and recently modified for Mars in preparation for the InSight mission. The M1a and M1b methods were modified from the previous work of Drilleau et al. (2013) and Panning et al. (2015Panning et al. ( , 2017. The M1c method is derived from Xu and Beghein (2019). The M2a method was already published in Khan et al. (2016Khan et al. ( , 2018, whereas the M2b and M2c methods were recently developed in the framework of the MSS.

M1 Models
M1a. The inversion algorithm is mainly based on the work of Drilleau et al. (2013) and Panning et al. (2015Panning et al. ( , 2017. The 1-D V S models are parameterized with three layers in the crust, and with six Bézier points in the mantle, which are interpolated using polynomial C 1 Bézier curves. The inverted parameters are the depth of the three layers in the crust and their V S values, the depth and V S values of the Bézier points in the mantle, and the V P /V S ratios in the three layers of the crust and the mantle. In total, 21 parameters are inverted to describe the velocity model. The a priori conditions on the depths of the crustal layers are that the depth of the first and third layers cannot exceed 10 and 100 km, respectively. The Bézier points are randomly located between the Moho depth and the top of the core. The V S profiles are randomly sampled within relatively wide prior bounds, as shown in Figure 8a. In order to ensure a velocity jump at the Moho, we impose that V S of the first Bézier point in the mantle must be higher than V S in the third layer of the crust. The V P /V S ratios are allowed to vary between 1.5 and 2.2. Simultaneously, a relocation of the quake is performed by moving the epicentral distance and the depth of the quake. For the blind test exercise, we choose to set the prior bounds relying on the uncertainties found by the MQS, between 80°and 95°and 25-45 km for the epicentral distance and the depth of the quake, respectively.
radial anisotropy mantle: composition mantle and core of a priori models parameters and temperature convective temperature based on geodynamical viscosity, parameters activation energy, activation volume 10.1029/2020EA001118

Earth and Space Science
To solve the inverse problem, we employ a Markov chain Monte Carlo (McMC) approach (e.g., Mosegaard & Tarantola, 1995). For each sampled model, we rely on the ray tracing algorithm of Shearer (2019) to compute body wave traveltimes. The surface wave velocity dispersion curves are calculated using the MINEOS package (Masters et al., 2011) and are then converted to traveltimes using the randomly sampled epicentral distance value. Since the origin time t 0 is typically unknown, we use differential times relative to the P wave phase arrival. The cost function is then defined as follows: where C computes the misfit between the observed and computed differential arrival times t S − t P , t pP − t P and the sum of the misfit between the observed and computed differential arrival times t S − t R at each frequency, taking into account the error bars σ P , σ S , σ pP , and σ R on P, S, pP, and Rayleigh waves arrival times, respectively. Superscripts throughout refer to observations (obs) and computed data (calc). Inversion output consists of an ensemble of internal structure models that fit the cost function. M1b. In contrast with the previously described method for M1a models, M1b models are inferred using information carried by surface waves only, considering a fixed location of the marsquake (see Table 1, computed origin). The anisotropy is taken into account in a joint Love/Rayleigh inversion scheme. The data space is composed here by two distinct ensembles of a priori pdfs of both Love and Rayleigh group velocities (Figures 6a and 6b; see section 3 for details). The parameter space in our case is composed by four subspaces (compressional and shear velocities, density, and ξ) as a function of depth. The inverse procedure relies on McMC explorations of the whole parameter space in order to compute a posteriori probability densities for each subspace. It is an improved version of the concepts presented in Panning et al. (2015Panning et al. ( , 2017 in the sense that the forward problem can now take into account for anisotropy. Under the hypothesis of a transversely isotropic medium (Love, 1892), the joint inversion of the two types of surface waves allows constraining where V SH and V SV are the horizontally and vertically polarized shear wave velocities, respectively.
For each 1-D trial model defined by ρ(z), V P (z), and V SV (z) and V SH (z), a modal summation theory (Masters et al., 2011) is used to compute the corresponding Love and Rayleigh group velocities. Each group velocity curve can be compared directly to the dispersion diagrams ( Figure 6), which quantifies the relevance of the trial model in the data space. This means that, in contrast to widely used misfit computations within Bayesian explorations and relying on gaussian assumptions, the goodness of fit is measured by the likelihood, The likelihood function is then the product of the individual probabilities sampled by the group velocity curves (which represent the state of a given parameter configuration a 1 , … , a m evaluated at each frequency ν i ).
To compute a large ensemble of model shapes, 70 Markov chains are running in parallel, and each chain uses a given geometrical constraint for the generated models between 0-and 140-km depth. As introduced by Drilleau et al. (2013), each model shape is controlled by piecewise C 1 Bézier curves (Bézier, 1977;Farin, 1993), based on several anchor points. Two anchor points at 0-and 140-km depth bound the exploration range and all intermediate points are shared by two following Bézier curves. To ensure a broad exploration of the model space, each Markov chain is associated with a unique random seed and an amount of Bézier points, which vary between 5 and 9. This amount varies between all chains but does not vary within a single chain during iterations. The minimum authorized distance between each anchor point is set to 8 km, which allows to generate models varying smoothly over the whole depth space as well as sharp discontinuities.
After a first stage of large wavelength exploration (cold runs), the posterior probabilities are constructed with the 6,000 models showing the lowest misfit values inferred during the 70 × 5,000 iterations (once they have been downsampled in order to prevent covariances).
M1c. The technique used here is a waveform modeling method based on a hierarchical transdimensional Bayesian approach to measure the dispersion of fundamental and higher mode surface waves. It was initially developed by Xu and Beghein (2019) for Earth applications and recently modified for Mars in preparation for the InSight mission. A McMC technique is used to seek a distribution of 1-D shear wave velocity models that represent the phase velocities of the fundamental and/or higher mode surface waves recorded between a seismic source and a receiver on a single seismogram. The distribution of 1-D models is then used to calculate dispersion curves and uncertainties, and tests are performed to assess the reliability of the measurements. Fundamental mode surface waves are generally much easier to isolate on the seismogram than the overtones, and thus, their dispersion is often more reliably measured.
An advantage of using a hierarchical transdimensional Bayesian approach lies in the fact the depth parameterization does not have to be fixed a priori since the algorithm lets the data control the complexity of the 10.1029/2020EA001118 Earth and Space Science solution while being parsimonious. Another great advantage is that it can also fit the data noise, which reduces the risk of mapping unknown noise into the velocity model and associated phase velocities. In addition, the source parameters can be included among the unknowns, thereby allowing source estimate uncertainties to be propagated into the model uncertainties. Since we are working with only one seismometer on Mars, uncertainties in the source may be larger than we are used to on Earth and it is well known that waveform inversions can be affected by trade-offs between source and structure. Here, we present results based on the source parameters obtained by the MQS (Table 1).
A notable difference between performing phase velocity measurements with waveform modeling on Mars compared to Earth is that we do not yet have a reliable reference model for Mars. For Earth applications, a reference dispersion curve calculated for a reference interior model is often used to prevent cycle skipping, which can affect the measured phase velocities. In the absence of a reliable reference model for Mars, we included the envelope of the waveform in the cost function instead, as was done previously by Yoshizawa and Ekström (2019).
The vertical component of the blind test data were filtered in different frequency bands, but no clear higher mode Rayleigh waves were visible. Fundamental mode Rayleigh waves, however, were clearly seen at periods between 25 and 50 s. The method employed here is not fully nonlinear due to the high computational cost of the forward problem. First, a synthetic reference seismogram and corresponding eigenfunctions are calculated using the fully nonlinear normal mode summation code MINEOS (Masters et al., 2011) and a starting model. The synthetic seismogram is then updated at each iteration using a linear approximation.
To find a reference model, we first tested several published 1-D Mars interior models (Sohl & Spohn, 1997;Zheng et al., 2015) and found that one of the Zheng et al. (2015) models predicted a synthetic seismogram that resembles the filtered blind data in the same period range the best, apart from a time shift of about 2 s. We then performed a rough grid search to modify the V S profile in order to bring the synthetic waveform and the blind test waveform closer together in time. The misfit was calculated with a L 2 norm. We used the resulting V S model as our reference model.
The V S profile is described by a variable number k of interpolation points, the vertical and horizontal positions of which define the depths at which V S is perturbed and the amount by which V S is perturbed relative to a reference model, respectively. The prior for V S is a uniform distribution of ±10% around the reference model (Figure 8c), and the Moho is allowed to vary by ±5 km around the reference value and two different reference values are tested to account for the prior uncertainty on crustal thickness. Perturbations in P wave velocity are assumed to be proportional to those in V S , as often done for Earth applications (e.g., Yoshizawa & Ekström, 2019). We compared results for which density anomalies were neglected and scaled to dV S (dρ/ρ ¼ 0.3dV S /V S ) and found no significant change in the resulting V S models.
We compared V S profiles obtained with reference models with different Moho depths, namely, 50-and 75-km depth. Note that by changing the crustal thickness in the reference model, we also had to adjust the reference V S in order to time shift the synthetic waveform and bring it closer to the observed blind surface waveform. We used a uniform prior for source parameters with a range of allowed values based on the uncertainties estimated by the MQS. Only the source latitude and longitude were kept constant. On average, the misfit for models with a 50-km-thick crust is smaller than that of models with a 75-km-thick crust. However, after performing F tests (Menke, 2012) on several of the best models in each case, we found that the misfit difference is not significant due to the difference in the number of model parameters. The waveform and envelope inversion is therefore unable to distinguish between models with a 50-or 75-km-thick crust. In the discussion in section 4.2, only the results using the 50-km-thick crust reference model are shown.
Because the envelope fit can be strongly affected by the noise level, we also tested whether the inclusion of group velocity data instead of the envelope affects the results. To do this, we modified the cost function to include group velocity measurements. We find that the velocity profiles do not strongly differ from those using the envelope measurements only, but the range of allowable V S models is slightly larger with group velocities than with the envelope.

M2 Models
M2a. In the following, we describe a method for inverting P and S wave body wave traveltimes and surface wave dispersion data jointly as outlined in Khan et al. (2016). To compute radial profiles of density, P and S

10.1029/2020EA001118
Earth and Space Science wave velocity, and shear attenuation, we rely on an average bulk Martian mantle composition and model areotherm using thermodynamic principles, mineral physics data, and viscoelastic modeling as described in Khan et al. (2018). Our Martian model is spherically symmetric and consists of three layers: a silicate crust and mantle and a metallic core. Seismic properties in the three layers are determined in the following manner.
1. Mantle. We use the free-energy minimization strategy of Connolly (2009) to determine stable mineralogy, elastic moduli, and density along self-consistently computed mantle adiabats for a given bulk composition. The thermodynamic formulation of Stixrude and Lithgow-Bertelloni (2005) including parameters as in Stixrude and Lithgow-Bertelloni (2011) are employed for this purpose. Bulk moduli are estimated by Voigt-Reuss-Hill averaging, while the pressure profile is obtained by integrating the load from the surface boundary condition p ¼ 10 5 Pa. Mantle compositions are explored within the Na 2 O-CaO-FeO-MgO-Al 2 O 3 − SiO 2 (NCFMAS) chemical model system, which accounts for more than 98% of the mass of the mantle of the experimental Martian model of Bertka and Fei (1997).Estimates for the Martian mantle composition derive from geochemical studies of a set of basaltic achondrite meteorites, which are believed to originate from Mars (e.g., Dreibus & Wanke, 1985;Lodders & Fegley, 1997;McSween Jr, 1994;Mohapatra & Murty, 2003;Sanloup et al., 1999;Taylor, 2013;Treiman, 1986). Based on these analyses, the Martian mantle is found to contain ∼17 wt% FeO, implying a Martian mantle Mg# of 75 (100 times molar Mg/Mg+Fe). 2. Crust. The crust is subdivided into three layers that are parameterized in terms of P and S wave velocity, and density, in addition to Moho thickness. To emulate the effect of porosity, we computed the aforementioned physical properties by multiplying the thermodynamically computed seismic wave speeds and density in crustal layer i by a variable parameter ϕ. The ϕ i is determined from ϕ i ¼ ϕ 0 + (1 − ϕ 0 )(i/N) with ϕ 0 being variable surface porosity and N the total number of crustal layers. This parameterization ensures that crustal properties increase from the surface down to the Moho where porosity is expected to vanish due to pressure. 3. Core. Sulfur is believed to be the dominant light element in the core of Mars because other elements (e.g., silicon, oxygen, and carbon) do not have sufficient solubility in iron-rich liquid at the relatively low pressures that are expected to have been maintained during core formation (e.g., Stevenson, 2001). Evidence in support of this comes from the observed depletion of chalcophile elements in the Martian meteorites, notably sulfur and the large value of the degree-2 gravitational potential Love number. To compute depth-dependent thermoelastic properties for the core, we use equations of state for liquid iron and liquid iron-sulfur alloys, relying on the parameterization of Rivoldini et al. (2011) for a well-mixed and convecting core. 4. Attenuation. The dissipation model adopted here is based on laboratory-derived torsional forced oscillation data on melt-free polycrystalline olivine and is described in detail in Jackson and Faul (2010). The extended Burgers model of Jackson and Faul (2010) is preferred over other rheological models because of its ability to describe the transition from (anharmonic) elasticity to grain size-sensitive viscoelastic behavior as a means of explaining the observed dissipation in the experiments on olivine.For present purposes, computations were conducted employing a single shear wave attenuation (Q) model at seismic frequencies (1 s) and a grain size of 1 mm. For the Martian crust and lithosphere, we fixed shear wave Q to 600 after PREM (Dziewonski & Anderson, 1981). As the core is assumed liquid, no shear attenuation is needed in the core. Dissipation in bulk is neglected and we assume Q κ ¼ 10 4 in line with terrestrial applications (e.g., Durek & Ekström, 1996). Anelastic P and S wave speeds as a function of pressure, temperature, composition, and frequency are estimated from the expressions for the viscoelastically computed temperature-, pressure-, and frequency-dependent moduli (further details may be found in Khan et al., 2018).
M2b. In this approach, we use a McMC joint inversion of body waves and surface wave seismic data, where the modeling of Mars's thermochemical history is part of the forward problem. In order to make reliable predictions about the seismic data from interior structure models the present-day thermal state of Mars is of fundamental importance. To obtain a plausible present-day thermal state of Mars, we simulate for each interior model its long-term thermal evolution. The thermal evolution is calculated with a parameterized thermochemical model that depends on a small set of parameters and an initial temperature profile. Such approach, allows the long-term planetary evolution to be accurately modeled at a reasonable computational cost.
Our forward problem consists in sampling the model space by computing different thermochemical evolutions of a Mars-like planet divided into several concentric and spherically symmetric envelopes ( Figure 9a): a convecting iron-sulfur liquid core (Rivoldini et al., 2011), surrounded by a silicate envelope. The latter consists of an adiabatic mantle (with top and bottom thermal boundary layers) convecting under a stagnant (i.e., essentially diffusive) lithospheric lid of thickness, D l , which includes a crust of thickness D cr , enriched in heat-producing elements. The thermochemical evolution was computed using a parameterized approach (e.g., Breuer & Spohn, 2006;Hauck & Phillips, 2002;Schubert & Spohn, 1990;Stevenson et al., 1983, and references therein) detailed in Samuel et al. (2019). In brief, the evolution of this layered planet is computed by numerically integrating a set of coupled differential equations expressing the conservation of internal energy within each convecting envelope, which includes internal heat production by radioactive elements and latent heating/cooling effects.
Heat transfer between the planetary envelopes is strongly controlled by the value of the mantle viscous rheology. Therefore, the viscosity of the Martian mantle plays an key role and is assumed to depend on both temperature and pressure, following an Arrhenius relationship (Karato & Wu, 1993). The temperature and the pressure dependence of viscosity is expressed by E * and V * , the effective activation energy and activation volume, respectively. These quantities can account for viscous deformation in the diffusion creep regime, or in the dislocation creep regime (Christensen, 1983;Kiefer & Li, 2016;Plesa et al., 2015;Samuel et al., 2019).
In the purely conducting lithospheric layer (in which the crust is embedded, see Figures 9a and 9b), the temperature profile is obtained by integrating the time-dependent heat equation with radially dependent heat sources, density, and thermal conductivity, to account for differences between the enriched and buoyant crust and the depleted residual lithosphere. The lithospheric and crustal thicknesses are not constant but evolve (Figure 9d) as a function of the time-dependent thermochemical state of the planet (Figure 9c). The lithospheric thickness is determined by considering an energy balance between the convective heat flux at the top of the mantle, the conductive heat flux out of the lithosphere, and the energy consumed to transform a portion of convective mantle into additional viscous lithosphere material, and vice versa (Schubert et al., 1979;Spohn, 1991, and references therein). The crustal thickness evolves by accounting for a time-dependent crustal production rate, the latter being a function of shallow mantle temperature and convective velocity (Breuer & Spohn, 2003;Spohn, 1991). Finally, present-day areotherms are obtained by conducting the parameterized thermochemical calculations for 4.5 Gyr (Figure 9b).
Along the McMC inversion procedure, we varied the values of the following governing parameters for the thermochemical evolution: • The mantle rheology: effective activation energy (E * ), volume (V * ), and reference viscosity (η 0 ); • the initial thermal state: temperature below the lithosphere (T m0 ) and core temperature at the CMB (T c0 ); • the core radius (R c ).
The values of all other physical parameters governing the thermochemical evolution of Mars were those listed in Samuel et al. (2019) (Tables S1-S3 in the supporting information), with the exception of the mantle and core average densities, thermal expansion, and specific heat. These quantities were obtained via the successive application of the thermal evolution and the mineral physics models in a fixed-point iteration fashion, until a simultaneous convergence of the values of these physical parameters and the thermal history was reached (typically less than ten iterations were necessary). This ensured consistency between the mineralogical and thermodynamic models, and the physical parameters used to compute the thermal evolution.
The parameterized approach described above reproduces the thermochemical evolution of a Mars-sized stagnant-lid planet in spherical geometry well at both transient and steady-state stages, including the

10.1029/2020EA001118
Earth and Space Science effects of complexities such as temperature and pressure-dependent mantle viscosities and the presence of an enriched evolving crust (Samuel et al., 2019;Thiriet et al., 2018, and references therein).
The present-day thermal profiles resulting from 4.5 Gyr of evolution are then used to compute seismic velocity profiles based on the mineralogical model described below. Underneath the crust, the areotherms were used to compute the mantle densities and seismic velocities using the Perple_X Gibbs energy minimization software (Connolly, 2005), which uses the thermodynamic formulation and thermodynamic properties at reference conditions of mineral phases of Stixrude and Lithgow-Bertelloni (2011), and assuming the mantle bulk composition of Taylor et al. (2006). To account for the influence of Mars's composition, we performed additional sets of inversions for which we considered different bulk Martian mantle compositions (Lodders & Fegley, 1997;Sanloup et al., 1999). Above the Moho depth, we parameterized the crust by considering several crustal layers of variable thickness, in which the body wave velocities and density are uniform. The uppermost part of the crust consisted of a 2-km-thick bedrock layer with a reduced density of 1,900 kg/m 3 and reduced seismic velocities (with V P and V S being set to 0.6 and 0.5 times the value of the corresponding velocities in the layer directly below it; Smrekar et al., 2019). In the liquid core, seismic velocities were computed following the approach underlined in Rivoldini et al. (2011).
Different thermal histories result in a variety of thermochemical structures, notably in terms of crustal and lithospheric thicknesses, or in mantle and core thermal states. These differences yield distinct stable mineralogical assemblages and therefore result in distinct seismic velocity profiles. Consequently, instead of varying independently the values of the density, or body wave velocities along a given profile depth during the inversion process, we sample the model space by varying the values of the governing parameters mentioned above (mantle rheology, initial thermal state, and core size) that control Mars's thermochemical history. This approach yields mantle and core density, seismic velocity profiles, and attenuation, as illustrated in Figure 10. Following the approach in Samuel et al. (2019) bulk attenuation is based on the PREM value, while shear attenuation is determined assuming a power law frequency dependence and a pressure and temperature Arrhenius dependence of the corresponding local quality factor Q μ , and by requiring that the present-day ratio of the planet's quality factor, Q, to its degree-2 Love number, k 2 , to be equal to 559 (Zharkov & Gudkova, 2005) at the tidal period of Phobos, together with an upper bound for Q μ in the mantle of 600 (Khan et al., 2018). In the liquid core, Q μ ¼ 0. The k 2 and Q values were computed following the approach for a viscoelastic medium described in Khan et al. (2004).
Unlike the mantle, the structure of the crust, its layering, and its seismic velocities are not computed from its chemical composition by Gibbs energy minimization because its composition is not well known, likely not in thermodynamic equilibrium, and heavily altered by surface processes. For this reason, we directly vary and invert for the crustal seismic velocity structure (both in terms of layering, and in terms of the values of seismic velocities within each crustal layer) instead of deriving it purely from thermal and mineral physics considerations because our mineralogical model does not apply for crustal conditions, and our approach only constrains crustal thickness and its density. However, the crustal density (within the range of [1,900-3,000] kg/m 3 ) together with the core sulfur content were adjusted together via a bisection method to satisfy the mass and moment of inertia constraints within uncertainties. Hence, these two parameters are not directly sampled by the inversion algorithm but are constrained through the aforementioned geodetic observations. In the cases where the bisection algorithm converged toward a crustal density outside of the above range, the corresponding evolution was excluded from the set of models. In addition, crater counting (Hartmann et al., 1999) indicates the presence of recent volcanism on Mars, which suggests that the interior of Mars is convectively active. Therefore, only evolutions that led to a convective mantle were retained (i.e., cases for which the mantle went subcritical were excluded; these correspond to cases for which the combination of temperature and rheological parameters led to a thin and relatively viscous convecting mantle). We also allowed for a constant shift in the obtained mantle seismic velocity profiles within ±5% in order to account for uncertainties in the thermochemical and mineralogical models. This was achieved by inverting for mantle correction factors, whose values ranged between 0.95 and 1.05.
In addition to allow for a better self-consistency than varying independently seismological parameters along the inversion process (Moho depth, seismic velocities, etc.), the built-in geodynamic frame significantly reduces the parameter space, by accounting for the interdependencies between various quantities (e.g., temperature, composition, and rheology) and their influences on crustal thickness. This approach can also 10.1029/2020EA001118 Earth and Space Science constrain the value of physical quantities inaccessible to direct measurements, such as the rheology of the Martian mantle, and provides the entire history of the planet associated with each model. These advantages are also tied to the main underlying assumptions in the forward geodynamic model. In particular, we assume is that stagnant lid mantle convection has operated for most of Mars history and until present day. In addition, we assume that Mars's mantle is homogeneous in composition. Compositional heterogeneity would affect both the thermal evolution and the seismic structure of the mantle (Smrekar et al., 2019). While the above assumptions are reasonable, it is important to keep in mind such framework when interpreting the results from these inversions.
M2c. This approach uses a set M of geophysically consistent a priori models provided by different authors (Khan et al., 2018;Rivoldini et al., 2011;Samuel et al., 2019). Generally, all models were constructed so as to satisfy current observations of mass, moment of inertia, and tidal response. The reader is referred to Giardini et al. (2020, SI3) for further details. This set originally consists of more than 20,000 models. To reduce this number for location operations, the following approaches were used: • Calculate a set of nine seismic variables for each model in the set. These observables are minimum and maximum thickness of the crust over the whole planet, crustal thickness at the landing site, P traveltime at 5°and 80°, S traveltime at 80°, extend of S shadow and surface wave traveltimes at two different frequencies. This creates a nine-dimensional space S of "seismic observable" with an injective projection M ← S. • Run clustering algorithm over all models, creating K clusters, such that each cluster has equal volume in S. This means different number of models N i in each cluster i, since they were importance-sampled to geophysical parameters by the authors. • Randomly select the same number of models n from each cluster to create a subset of S, called S sel . This subset spans a wide range of potential seismic observables, but the importance sampling property is lost. • Reapproach the importance sampling character by giving models in each cluster the weight w i ¼ N i /n. For operations, MQS selected 2,500 models of the full set by cutting S into K ¼ 10 clusters. Once seismic observations are available, a grid search over possible depths d, distances Δ, and m j ∈ S sel is done to calculate likelihoods for each combination L(d, Δ, m j ). By multiplication with the prior weights w i and integration over depths and distances, a marginal probability for each model is computed.

Inversion Results and Discussion
The inversion results of the six methods are displayed in Figure 11. They show the additional gain in information obtained through inversion, compared to the prior distributions ( Figure 8). The thick black line represents the 1-D base model used to compute the synthetic seismograms, whereas the gray area shows the minimum and maximum V S values encountered over the great circle path between source and receiver. The inverted data sets are the body waves and surfaces waves arrival times estimated by the MQS (Table 1 and Figures 1 and 5), except for the M1b and M1c methods, which use group velocity dispersion diagrams ( Figure 6) and both waveforms and group velocities, respectively. The data fit is shown in Figure 12. Figure 11 clearly highlights the nonuniqueness of the solution. Depending on the depth, some distributions are unimodal, bimodal, or multimodal. The results show that none of the distributions perfectly fit the input model, but depending on the method, they share some common feature with the model to retrieve.
Let us first review the results of the three M1 methods. The pdf of the M1a method (Figure 11a) is spread within the parameter space, because both the seismic velocity profiles and the quake location are simultaneously inverted, with relatively large prior bounds. The maximum of the pdf between 15-and 70-km depth is located in the vicinity of the input model, which means that V S in the crust is well constrained. Below 80-km depth the distribution enlarges and has a similar shape to the a priori pdf (Figure 8a), which means that the sensitivity to the data decreases at these depths (the highest period of surface wave considered with this method is 40 s; Figure 12). The discontinuity located at 100-km depth is due to the prior assumptions and the loss of sensitivity of surface waves at these depths. The pdf of M1b (Figure 11b), obtained using surface waves measurements only, is narrower compared to the previous method, mostly because the source location and the origin time are fixed to the MQS values. Although surface wave data are known to be less sensitive to the sharpness of the discontinuities, two inflexion points are observed near 15-and 60-km depth, close to the two discontinuities of the input model. The V S distribution in the crust lies within the 10.1029/2020EA001118 Earth and Space Science minimum and maximum V S values encountered along the first orbit Rayleigh wave path. Similarly to the M1a method, the pdf is spread below 70-km depth. This enlargement is also visible on the results obtained with the M1c method (Figure 11c), which invert both surface wave group velocities and surface wave waveforms. The retrieved V S value in the crust is constant with depth, overestimated in the upper crust and slightly underestimated in the lower crust, which means that the M1c inversion is unable to recover the detailed layering of the crust, which is not surprising considering the sensitivity of Rayleigh waves at the periods measured (between 25 and 55 s; Figure 12). The M1c pdf shows a discontinuity clearly located at 50-km depth, in the range of the model to retrieve, because the crustal thickness is allow to slightly vary between ±4 km around the 50-km-thick crust 1-D model chosen as the starting model of the inversion. The M1c method is also able to provide the strike, dip, slip, and the focal depth of the

10.1029/2020EA001118
Earth and Space Science marsquake. Figure 13 shows that the strike's distribution is centered near 110°, in agreement with the MQS values. In return, the focal depth is not constrained by the data, and changes in the marginal distributions for the dip, slip, and Moho depth indicates trade-offs among the parameters.
By analyzing the output distributions of the M1 models, we conclude that (1) V S in the crust is relatively well constrained and consistent to the mean path averaged values; (2) the structure is not reliably constrained below 70 km; (3) the crustal thickness is difficult to retrieved by inverting both the V S profile and the quake location (M1a), but it could be approximated if the quake location is fixed (M1b) or if small perturbations around a 1-D starting model with a Moho depth already close to the model to retrieve are performed (M1c).
Concerning the three methods handling geophysical and geodynamical parameters, V S in the crust is slightly underestimated for M2a models (Figure 11d) and overestimated for M2c (Figure 11f). Note that

10.1029/2020EA001118
Earth and Space Science posteriori distribution of the Moho depth is almost equal to the a priori distribution. This result is in good agreement with the M1 outputs, showing that the Rayleigh waves have a poor sensitivity at these depths. The M2c models (Figure 11f) display a distribution of Moho depth compatible with the range of Moho depths encountered along the first orbit (Figure 14f ). It should be noticed that the maximum value of the a priori distribution is located in the range of the true Moho depth values. For the three M2 methods, the seismic velocities in the mantle are narrower than for M1 methods, because they rely on geodetical, geophysical, and thermochemical modeling. These assumptions considerably limit the number of possible solutions in the mantle, as shown for the M2a models, where the seismic velocities in the mantle are nearly fixed for all the sampled models (Figures 11d and 8d). An interesting feature in the M2b models is that V S in the mantle are systematically higher than the model to retrieve (Figure 11e). Indeed, the a priori pdf of V S in the mantle (Figure 8e) is significantly higher than the input model, which means that the algorithm will never reach the input model. It seems that the modeling of Mars thermochemical history is not compatible with such an extreme 1-D base model, build using a hot temperature profile. Concerning M2c models, the V S distribution in the mantle is bimodal, similarly to the a priori distribution (Figure 8f) but the maximum of the pdf value is located in the vicinity of the profile to retrieve (Figure 11f).
By investigating the a posteriori distributions of the M2 models, we highlight the following points: (1) V S in the crust is constrained if the prior bounds are large enough; (2) the geophysical and/or geodynamical a priori assumptions could have some difficulties to retrieve the model if Mars thermal properties are greatly anomalous, compared to the knowledge we have of the Earth (M2b); (3) considering the weak sensitivity of surface waves near Moho depth, the crustal thickness is approximated if the a priori distribution is relatively close to the true value (M2c).
Although it is reductive to represent the complexity of multimodal pdfs using a simple mean and a 1σ standard deviation, in order to compare the results, the mean models of each method and their associated standard deviation are displayed in Figure 15. Figure 12a shows the mean differential arrival time between the arrival times of Rayleigh waves at the different periods and the arrival times of the P wave for all the sampled models. We observe that almost all the models fit the data within uncertainty bounds, which means that Figure 15. Mean V S models of the pdfs shown in Figure 11, for each of the six methods. Solid and dashed lines are the mean profiles ± mean absolute deviation obtained from all sampled models. The thick black line represents the 1-D base model used to compute the synthetic seismograms, whereas the gray area shows the minimum and maximum V S values encountered along the first orbit Rayleigh wave path.

10.1029/2020EA001118
Earth and Space Science almost all the sampled models are able to explain the data. The crust is mainly constrained by surface wave group velocity, and the M1 models ( Figure 15) show a relatively good agreement of V S in the crust. The enlargement of the three M1 distributions near 70-km depth indicates that the structure is poorly constrained at these depths. The uncertainties of the M1a model are large compared to the M1b and M1c models, because the quake location is inverted at the same time. In the crust, the M2a and M2c models are located on the extrema, mainly due to the a priori conditions considered (Figures 8d and 8f). The M2b model is closer to the input model due to the large a priori values considered in the crust (Figure 8e). However, the mean V S model obtained by the M2b algorithm has higher seismic velocities than expected in the mantle, which explain why the maximum of the t S − t P distribution is lower than the measurement (Figure 12b), as well as the t R1 − t P mean value at 33.6 s (Figure 12a). It should be noticed that reaching the input V S value in the mantle with this algorithm should be possible by modifying the prior condition, but this would enhance mantle melting and result in thicker crust. The small uncertainties of M2 models in the mantle are mainly due to the strong a priori assumptions considered.
The results of the quake location are displayed in Figure 16, for the four approaches that are relocating source during the inversion process (M1a and M2 methods). Given that the true epicentral distance is relatively large (87.8°), the body waves are traveling in the mantle and so they mostly constrain the 600-900-km depth range. It is particularly enhanced in the M2b case: In order to increase the t S − t P differential time due to higher velocities in the mantle, the epicentral distance distribution is shifted toward the upper bound of the prior range (Figure 16a). M1a distributions of the epicentral distance, depth of the quake and origin time are spread, because the seismic velocities are explored in a wide parameter space (Figure 8a). The depth of the quake is well retrieved, and the epicentral distance is slightly overestimated (90 instead of 87.8°). The origin time distribution is large compared to the other ones, because no a priori assumption is made on it. The M2a method retrieve well all the location parameters, because the V S values in the mantle are very close to the input model (Figures 11d and 15) due to strong a priori conditions (Figure 8d). The M2c epicentral distance and depth of the quake distributions are slightly shifted toward higher values, because the V S distribution shows globally higher seismic velocities compared to the input model, in both crust and mantle.
We clearly observe on these blind test results the complementarity of the two parameterizations M1 and M2, to ensure sufficient coverage of the a priori model space. The relatively large prior bounds used by the M1 models, with no constraints from mineral physics and thermodynamics, allow more flexibility in the model sampling. They are powerful to constrain V S in the crust and also indicate regions of the parameter space in depth, where the data (mainly surface waves) are no more sensitive to the structure. The M2 models produce stable velocity models through the whole planet, in the sense that they satisfy current observations of mass, moment of inertia, and tidal response, considering a variety of starting bulk Mars compositions, and give a physical interpretation of the data. In the absence of surface wave measurements at large periods, allowing a better exploration of the Moho, they also explicitly output Moho depths only compatible with physical relations. However, they could show some limitations if Mars's structure reaches extreme bounds of what we expect, based on the knowledge on the interiors of Earth and Moon. The differences between the results from the six different methods are mainly due to the chosen prior assumptions. Comparisons with the prior distribution are essential to understand how the data improve our knowledge compared to the a priori conditions. Considering only one station, these results support the approach set up by the MSS of using different algorithms to constrain the deep interior structure of Mars.

Inversion Results Considering Only Body Waves
Seventy-two sols after landing, InSight began streaming continuous data to Earth. During the first 168 sols, until 31 July, InSight recorded 70 events, with at least 12 events of magnitude Mw3-4  and two larger events located near the Cerberus Fossae system at 25-30°distance. Only direct P and S phase identification is possible for S0173a and S0235b events, as scattering prevents positive identification for all other events Lognonné et al., 2020). Multidiffusion analysis have shown that the scattering is intermediate between that of the Earth and Moon . The first shear Q measurements provide nevertheless values, which are in accordance with the lithospheric Q for these two large events and are 2-3 times smaller than the lunar crustal Q .
For single station analysis, the absence of surface waves is however making the determination of both the quake location and the velocity structure more challenging.
To respond to these new observations and to show what could be investigated using body waves alone, algorithms that are using both surface and body waves are adjusted by removing the surface wave contribution. The posterior pdfs are displayed in Figure 17. The M1b and M1c methods, handling surface waves only, are not represented. The location and origin time are shown in Figure 19 and the data fits are displayed in Figure 18. The exploration of the whole sets of models allows to check the compatibility of each model with the recorded data. We observe that for all the methods, the posterior pdfs ( Figure 17) are very similar to the prior pdfs (Figure 8). The results of the M1a method show that the pdfs of the epicentral distance and quake location are almost identical to the prior distribution, which clearly demonstrate the trade-off between the V S profile and the quake location. The M2a and M2c distributions of the epicentral distance (Figure 19a) are broader compared to inversions made using both surface waves and body waves (Figure 16a). The M2b method always indicates a larger epicentral distance and lower tS − tP values. Constraining the velocity structure using body waves arrival times only seems difficult, in the sense that the algorithms will almost always find a combination of V S model and epicentral distance compatible with the data. In other words, the output epicentral distance distribution reflects the range of epicentral distances compatible with the a priori distribution of V S profiles.

RFs Calculation
Five different methods are used to computed P-to-S RFs. A comparison of results, both in the ZRT coordinate system and in the ray-centered LQT system that includes an additional rotation around the incidence angle so as to effectively separate P and SV contributions to the wavefield, is shown in Figure 20. A brief summary of the methods is given below. All of them have been applied to Earth data before, though in the context of having many tens to hundreds of usable events.
For application of the transdimensional hierarchical Bayesian deconvolution (THBD) method (Kolb & Lekić, 2014), the P wave coda of the event was first transformed into upgoing P and SV waveforms by applying a free-surface transfer matrix (Kennett, 1991). The elements of the transfer matrix depends on the ray parameter as well as the near-surface P and S wave velocities. Their values were estimated by minimizing the energy on the SV component during the first 2 s of the P wave arrival, following Abt et al. (2010). The best estimates for near-surface P and S wave velocities depend on the length of the deconvolution window and whether it includes the pP phase or not. The shorter window gives a ray parameter estimate of 3.7 s per degree, a P wave velocity of 4.3 km/s, and an S wave velocity of 2.5 km/s, whereas the longer window results in a ray parameter estimate of 4.5 s per degree, a P wave velocity of 4.2 km/s and an S wave velocity of 2.1 km/s. The THBD method then uses a reversible jump McMC algorithm (Green, 1995) to sample one million random realizations of RFs, varying the width, lag-time, amplitude, and number of Gaussian pulses that make up the RF waveform. Resulting RFs were band passed between 0.1 and 0.7 Hz. The result of the THBD method is an ensemble of RFs that are compatible with the data. Features that are common across the ensemble are considered as robust. The THBD results for the shorter deconvolution window are displayed in Figure 20b as gray shading for both P and SV, with black lines indicating the average RFs.  Table 1. The deterministic deconvolution methods applied include Wiener filter deconvolution in the time domain (Hannemann et al., 2017;Kind et al., 1995), water-level deconvolution (Clayton & Wiggins, 1976), iterative time domain deconvolution (Ligorria & Ammon, 1999), and extended-time multitaper frequency domain deconvolution (Helffrich, 2006).
For the Wiener filter deconvolution, the data were band-pass filtered between 20 s and 2 Hz. RFs were calculated both in the ZRT and in the LQT system. The incidence angle for rotation into the LQT system was estimated by polarization analysis of the P wave onset (Jurkevics, 1988), yielding an angle of 19°. A length of 80 s for the deconvolution window and a damping factor of 0.1 were used, which resulted in clear RF phases with little noise contamination. We also derived the frequency-dependent apparent P wave polarization by measuring the amplitudes of the resulting vertical and radial RF in different filter bands (see below).
Water level deconvolution was applied to obtain a radial RF, whereas the vertical RF was obtained by spectral whitening and autocorrelation (Tauzin et al., 2019). Gaussian low-pass filter were applied during deconvolution, where the width of the Gaussian was determined by the parameter a ¼ 2.5 rad/s for the radial and a ¼ 3.0 rad/s for the vertical component (Ammon, 1991). The smoothing width for the spectral whitening was set to W ¼ 0.05 Hz. To stabilize the long-period component of the vertical RF, a zero-phase second-order Butterworth high-pass was applied at 5-s period.

10.1029/2020EA001118
Earth and Space Science visible P wave arrival. Deconvolution is then performed using a 50-s window, a time-bandwidth product of 3 that translates to a frequency bandwidth of permissible spectral leakage of 0.2 Hz, and four tapers (Shibutani et al., 2008). Different time windows for the source function were tested and found to produce similar results.
The P wave coda of the synthetic event contains a strong secondary arrival at about 14 s after the onset, identified as pP (section 3 and Figure 1b), which is also visible in the RFs. When using the ZRT coordinate system, this results in a negative pulse on both the vertical (Z) and radial (R) component for some deconvolution methods (Figure 20a). Even though the amplitude of this pulse on the R component is comparable to that of other phases, for example, that at 8 s, this is not a P-to-S converted phase that we want to interpret, but rather a mapping of P energy onto the radial component. Rotation in the LQT or the P-SV coordinate system efficiently suppresses this arrival ( Figure 20b) and thus provides a better base for interpreting the RF in terms of P-to-S conversions.
The frequency content of the RFs from different deconvolution methods is not completely identical as the parameter that influence the frequency content were chosen individually for each method. The timing of the main phases in the radial and Q component RFs waveform, that is, at 1.8 s, between 5 and 10 s, and at 26 and 34 s, respectively, shows good agreement between most of the applied methods, though the amplitudes of the iterative time domain and extended-time multitaper deconvolution results are smaller, specifically for the early arrivals. The Q-RF waveforms of the various deterministic deconvolution methods generally lie within the uncertainty range provided by the THBD.
The results of the Wiener filter deconvolution were further used to measure the frequency-dependent apparent P wave polarization. This quantity is related to the apparent subsurface S wave velocity via a simple relation that also requires knowledge of the ray parameter (Hannemann et al., 2016;Svenningsen & Jacobsen, 2007;Wiechert, 1907). The derived apparent S wave velocities were jointly inverted with the RF waveform. Unfortunately, the signal-to-noise ratio of the radial component, as defined in Knapmeyer-Endrun et al. (2018), was only larger than 5, which is considered an indicator for a reliable measurement (Hannemann et al., 2016), for periods up to 4 s. Thus, the apparent velocity curve is only defined at the shortest periods and provides constraints on the S wave velocity in the near-surface layer only. The apparent P wave incidence angles range between 19°and 20°for these periods, in good agreement with the angle determined Figure 21. Amplitude stacking according to Zhu and Kanamori (2000). (a) Results for the average RF waveform from the THBD method using semblance weighting (e.g., Knapmeyer-Endrun et al., 2014), and a value of 5.5 km/s for the average crustal P wave velocity. The assumed ray parameter is 3.7 s per degree. (b) Results for the Wiener filter deconvolution using semblance weighting (e.g., Knapmeyer-Endrun et al., 2014) and averaging results for three different average crustal P wave velocities: 5.5, 6, and 6.5 km/s. Lower P wave velocities result in smaller values for the crustal thickness, and higher P wave velocities in larger values. The assumed ray parameter is 4 s per degree.

10.1029/2020EA001118
Earth and Space Science from polarization analysis of the actual waveform data. To derive apparent S wave velocities from these polarization angles, information on the ray parameter is needed. Based on the estimated epicentral distance of 85.7°and precalculated P wave ray parameters for a number of Mars models and event depths (Knapmeyer-Endrun et al., 2018), the estimated likely range is 3.6 to 4.0 s per degree.
In an initial interpretation of the RF waveforms, the working group unanimously identified the clear positive arrival at 8.0-8.3 s as the likely Moho conversion, with the following peak and trough around 26 and 34 s as the corresponding first (PpPs) and second (PsPs + PpSs) multiply reflected and converted phase, respectively. This interpretation is further substantiated by applying the stacking method of Zhu and Kanamori (2000) to the single RF trace. This methods stacks the amplitudes at the arrival times for hypothetical Ps conversions and their multiples, which are calculated for a given crustal thickness, V P /V S ratio and average crustal P wave velocity. The phases should stack coherently and thus produce an amplitude maximum for the true crustal thickness and V P /V S ratio. In general, the method benefits from the availability of data from multiple events with different ray parameters. However, it also provides a clear result for the single synthetic trace, pointing to a Moho depth of 52 to 66 km, depending on the assumed average crustal P wave velocity, and a V P /V S ratio, less well constrained than the Moho as usual, between 1.73 and 1.95 ( Figure 21). The assumed average crustal P wave velocity is found to have a larger influence on the actual resulting crustal thickness than the assumed ray parameter of the event, at least for variations within a sensible range (3.6-4.0 s per degree).
In contrast to the broad consensus on the identification of the Moho signal, there were different interpretations for intracrustal structure. The first interpretation considers the peak at 1.6-2.2 s as intracrustal conversion, with the peak and through at 5.5-6 and 7.1 s as the corresponding multiples, whereas an alternative interpretation saw the peak at 5.5-6 s as an additional intracrustal conversion. This is also mirrored in the stacking results that strongly point to a discontinuity around 10-km depth and also contain some constructively stacked energy around 40-km depth (Figure 21b). The waveform inversion accordingly also considers both two and three crustal layers (section 5.2).
We did not try to calculate S-to-P RFs for this event because the S arrival is rather weak, in keeping with some of the synthetics of the previous MQS blind test, which contained a number of models with broad S wave shadow zones in the mantle .

Inversion
A single algorithm was applied in order to invert the RFs for crustal structure. It was specifically developed for the application to small amounts of data, as, for example, InSight. Another algorithm, using a fully Bayesian approach, has been presented for application to InSight data from Mars before (Panning et al., 2017) but was not used during the blind test. Similar to the surface wave case, both algorithms produce an ensemble of models that can explain the data in order to also characterize uncertainty and trade-offs. The main characteristics of the method applied here are given in Table 3. All models adhere to a standard seismic parameterization as the crustal composition of Mars is expected to be variable and possibly, due to the low temperatures, not in thermodynamic equilibrium (e.g., Wood & Holloway, 1984). Additional complications arise from complex crustal lithologies and porosity. Thus, mineral physics constraints are weak and all considered models are parameterized in terms of seismic velocities and density as a function of depth.
We jointly invert the radial RF waveform resulting from Wiener filter deconvolution and the apparent S wave velocity curve derived from this waveform, as previous studies have shown that this combination helps to reduce the depth-velocity trade-off inherent to RFs as a traveltime method (Svenningsen & Jacobsen, 2007) and that the apparent velocity curves from just a few events should allow to recover the first-order crustal structure of Mars (Knapmeyer-Endrun et al., 2018).
The inversion is performed with the Conditional Neighbourhood Algorithm (NA) (Wathelet, 2008). The NA (Sambridge, 1999) is a direct search method that uses Voronoi cells to subdivide the parameter space and Note. For each layer, the depth, S wave velocity, and V P /V S ratio are abbreviated as D, V, and R, respectively. V h and R h represent the velocity and V P /V S ratio of the half-space.

10.1029/2020EA001118
Earth and Space Science preferentially looks for new models in the vicinity of the current best fitting models in a self-adaptive manner. In principle, it is able to identify several separate low-misfit regions during a single inversion run. The NA is a derivative-free method, so the complexity of the misfit function is the computationally limiting factor, as its derivatives are not required to guide the search (Sambridge, 1999). The modifications by Wathelet (2008) allow to define irregular limits to the searchable parameter space, for example, based on physical conditions or prior information, and include a dynamic scaling to keep the exploration of the parameter space as constant as possible over the iterations of an inversion run.
The forward calculation of RFs is based on the code by Shibutani et al. (1996) that uses a simple reflectivity matrix approach to provide the P-to-S response of a layer stack. The resulting synthetic vertical and radial RFs are convolved with the measured vertical RF to take into account source complexity, and, in this case, the pP conversion in the P wave coda. Tests with synthetic seismograms for Mars models have shown that this allows to obtain results comparable to full Instaseis synthetics based on an AxiSEM database (Nissen-Meyer et al., 2014) with a greatly reduced computation time (Knapmeyer-Endrun et al., 2018). Density was not used as an independent parameter during the inversion but calculated from P wave velocity values using Birch's law (Birch, 1961). V P /V S , on the other hand, was allowed to vary. For details on the parameterization, see Table 3.
Cases with both two and three constant-velocity layers over a half-space were investigated. Individual forward runs prior to the inversion process showed that misfit for the RF waveform was an order of magnitude less than that for the apparent velocity curve, so a weight factor of 10 was used for the RF misfit. As there was no standard deviation of the data available based on only a single event, the misfits were not weighed inversely with their standard errors in the objective function definition. The cost function is then defined using an L1 norm: where n iterates over the number of data points in the RF and m over the number of samples in the apparent velocity curve. Superscripts refer to observations (obs) and predicted data (pred). Due to the uncertainty in event location and velocity model, there is an uncertainty in the ray parameter which was estimated to lie between 3.6 and 4.0 s per degree. Both end-members were used to invert for the subsurface structure as they would span the range of all the possible ground models.
Two parameters that control the NA need to be tuned depending on the style of sampling needed. For a more explorative search that is robust against local minima, we perform 1,200 iterations in each inversion run with 500 models produced at each iteration and 300 cells resampled at each iteration, resulting in an ensemble of ∼400,000 models per run. Furthermore, each inversion was repeated several times to test the stability of the results. Figure 22 shows the resulting 1-D subsurface profiles for the two-layer and three-layer case, with the fit to the RF waveform and apparent velocity curve for the two-layer case given in Figures 23 and 24, respectively. The plot includes all models within a maximum misfit, ranked and color coded according to misfit. The maximum acceptable misfit is derived visually by considering the uncertainty in the RF waveform as derived by THBD ( Figure 20). The spread between the maximum and minimum RF amplitude at each point in time as given by THBD is added to the inverted RF waveform as a measure of uncertainty (dashed lines in Figure 23) and all models that produce RFs that mostly lie outside of this range are discarded.

Inversion Results and Discussion
It is evident from the velocity profiles that adding a third crustal layer to the model is not warranted by the data. The velocity contrast between the additional layer and the one adjacent to it is about 1.5% for both ray parameters considered, and the velocities obtained for this layer lie within the velocity range of the adjacent layer. The minimum misfit also just reduces insignificantly from 0.0065 to 0.0061 in going from two crustal layers (8 free parameters) to three layers (11 free parameters). Adding yet another fourth crustal layer again has no visible effect on the ground profiles, and therefore, a two-layer model is considered as the optimum for the inversion, in agreement with the true model. The V P /V S ratio could not be well constrained by the inversion, that is, models within the considered misfit range show V P /V S values across the whole allowed parameter space (Table 3) for all layers. Therefore, the V P /V S results are not shown here.
The S wave velocities and layer depths of the true model lie between the model families derived for the two end-member ray parameters. The uncertainty in the velocities gets larger at larger depth, though, and the S Figure 22. One-dimensional velocity profiles from joint inversion of receiver functions and apparent velocity curves for (a) two crustal layers and (b) three crustal layers; (i) and (ii) denote ray parameters 4.0 and 3.6 s per degree, respectively. The thin black lines represent the minimum and maximum parameter ranges, the thick black line represents the true model and the thick dotted brown line indicates the near-surface S wave velocity obtained from the free-surface transfer matrix for a ray parameter of 4.5 and 3.7 s per degree in (i) and (ii), respectively. wave velocity in the uppermost mantle is less well resolved than the crustal velocities. In fact, the true mantle S wave velocity is at the upper end of those resulting from the inversion. As there are no conversions from within the mantle, this velocity is not constrained by any traveltimes, like those of the intracrustal and Moho conversions, but only by the amplitudes of these phases. As amplitudes have a higher uncertainty than traveltimes in RF, this is one probable explanation. Another factor that contributes to the loss of resolution with depth is the limited period range for which apparent velocities could be measured. Nevertheless, crustal velocities, layering and thickness could be reasonably well determined form a single event and with a limited apparent velocity curve.
The S wave velocities derived from the free surface transfer matrix for the shallow subsurface show a clear trade-off with the associated ray parameters. For a ray parameter of 4.5 s per degree, which is larger than the range estimated based on epicentral distance and model uncertainty, the estimated velocities are too low, whereas they show a closer agreement with the inversion results and the true model for a ray parameter of 3.7 s per degree. For the actual marsquakes, MQS is providing probabilistic ray parameter distributions for each event, which was not the case here. This allows to better take the uncertainty in the ray parameter into account during the inversion and also provides a reference to compare the results from the free surface transfer matrix to. Thus, it helps to understand better how the velocities derived from both approaches might deviate from each other.   Table 4.

10.1029/2020EA001118
Earth and Space Science

Attenuation and Diffusion
Another seismic feature we would like to extract from the seismic signals is the attenuation, which is represented by the quality factor Q. Attenuation is a feature closely related to the thermal structure and the composition. One variable that increases the attenuation is the existence of the water (e.g., Karato, 2013). Observed seismic amplitude spectra A can be expressed as where ω is the angular frequency, t is the total traveltime, and A 0 is the original amplitude without attenuation (Aki & Richards, 2002). The quality factor Q can be obtained by fitting the model to the observed seismic spectra. Since the attenuation also depends on the traveltime, which varies with the seismic velocity model chosen, attenuation may also be expressed by t * , which is the ratio between Q and the traveltime t. It should also be noted that this Q should be regarded as an effective Q with multiple phenomena included. The quality factor Q will quantify the amount of energy loss during the wave propagation, including both elastic and inelastic effects. The elastic attenuation can be referred to as scattering due to the heterogeneity in the medium, while the inelastic attenuation can be referred to absorption of energy (e.g., Sato et al., 2012). For a layered structure, Equation 5 can be rewritten as and by using multiple events with different paths as inputs, we can invert for the 1-D structure of Q. In the blind test, however, we are focusing on what we can investigate with a single event and we will estimate the effective Q by fitting the spectral decay. With the real Martian data, we plan to investigate both the elastic and inelastic effect by using methods such as the one described in Gillet et al. (2017).
Here, we use the arrival picks provided by MQS and chose P, PP, and SH phases for Q estimation. Each phase arrival was cut into 5-to 10-s window chosen to include the complete signal pulse of each phase. Then seismic spectra were calculated using a multitaper method. Figure 25 and Table 4 show the results of the spectral fit and the summary of t * and Q obtained with the reference values used to create the synthetic data. Q was calculated using traveltimes calculated with the reference seismic velocity model. We found that our method was able to obtain correctly the t * and/or Q for the P phase while the estimation was less accurate for the S phase. These results show that the simple approach we have taken here is efficient to obtain a first estimate of the attenuation of Mars even with single event.

Conclusion
In this article, we describe the methodologies that have been proposed by the MSS, based on a single station and a single event, and we demonstrate the feasibility for locating Marsquake and for constraining the 1-D seismic velocity structure of Mars. This be achieved using RFs, and surface and body wave arrival times information. We also show how we could estimate the attenuation of Mars. The inversion of the RF waveforms and the apparent S wave velocity, taking into account the uncertainty of the ray parameter, provides the velocity profile and the Moho depth below the InSight station. Using six different methods, the combination of body waves and surface waves enable to determine V S in the crust integrated along the great circle path. Taking into account the 1σ uncertainties (Figure 15), the models range between −6% and +9% around the 1-D base model. Due to the period range of surface waves considered (and because high-frequency surface wave are contaminated by an extra energy, see Figure 6), the Moho depth and V S in the mantle are poorly constrained. Only the methods that are not inverting source and medium properties simultaneously are retrieving a fairly good estimation of Moho depth. We demonstrate the complementarity of using both models parameterized with seismic velocities and physical parameters, when handling such a nonlinear problem. Although the uncertainties on the structure are large, the approach chosen by the MSS, of combining the outputs of several methods, using different a priori conditions and different seismic data, helps to reveal biases of any single method and highlights the diversity of the models compatible with the data. However, complexities related to anisotropy and three-dimensional structure, particularly in the crust and lithosphere, undoubtedly complicate the interpretation as a 1-D radial model envisaged here. Lateral structural variations in the Martian crust are likely much stronger than on Earth due to the crustal dichotomy between the Southern and Northern Hemispheres (Zhong & Zuber, 2001), which complicates the wave propagation on Mars.
In the absence of surface waves, as is the case for the Marsquakes recorded by InSight so far, we demonstrated the difficulties to constrain independently the velocity structure and the source location. A future effort be made in order to create a joint inversion of both RFs and body wave arrival times. While the identification of seismic phases is up to now limited with the InSight data, if a number of additional phases could be identified, for example, reflections from the surface (SS and SSS) and Moho (PmP and SmS), they could add additional constraints on the structure and quake location. In the event of meteorite impacts, they could be located by one of several orbiting cameras, which could provide a known location. This would enable the direct inversion of all differential traveltimes with respect to the P arrival time and improve the constraints on the output 1-D model. Daubar et al. (2018) demonstrated that using a fixed location, the average seismic model could be constrained near the depths of the turning point of the ray paths. If more events with clear wave arrivals accumulate, joint inversion of all the events will be critical, enabling a continuous refinement of model parameter distributions. However, the feasibility of the present approach will strongly depend on the future data. Further work seems to be necessary to develop new techniques using, for instance, polarization and amplitude information, to fully investigate the recorded data.