Amir Khan
email: akhan@ethz.ch

Savas Ceylan

Martin Van Driel

Domenico Giardini

Philippe Lognonné

Henri Samuel

Nicholas C Schmerr

Simon C Stähler

Andrea C Duran

Quancheng Huang

Doyeon Kim

Adrien Broquet

Constantinos Charalambous

John F Clinton

Paul M Davis

Mélanie Drilleau

Foivos Karakostas

Vedran Lekic

Scott M Mclennan

Ross R Maguire

Chloé Michaut

Mark P Panning

William T Pike

Baptiste Pinot

Matthieu Plasman

John-Robert Scholz

Rudolf Widmer-Schnidrig

Tilman Spohn

Suzanne E Smrekar

William B Banerdt

Upper mantle structure of Mars from InSight seismic data

used observations of direct (P and S) and surface-reflected (PP, PPP, SS, and SSS) body wave phases from 8 low-frequency marsquakes to constrain the interior structure to 800 km depth. We found a structure compatible with a low-velocity zone associated with a thermal lithosphere much thicker than on Earth that is possibly related to a weak S-wave shadow zone at teleseismic distances. By combining the seismic constraints with geodynamic models, we predict that the crust is 13-20 times more enriched in heat-producing elements than the primitive mantle, greater than suggested by gamma-ray surface mapping, and has a moderate-to-elevated surface heat flow.

The Interior Exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) mission (1) touched down on Mars in Elysium Planitia (2) at the end of 2018, and has been acquiring high-quality seismic data with the Seismic Experiment for Interior Structure (SEIS) instrument (3) since early 2019 (4,5). SEIS's main instrument is a 3-component very broadband seismometer [START_REF]Materials and methods are available as supplementary materials[END_REF], and to date (Sol 676), >1000 distinct seismic events have been identified by the Marsquake Service (MQS) [START_REF] Clinton | [END_REF]. A primary goal of the InSight mission (1) is to image the interior structure of Mars from observations of seismic events and to use this to improve our understanding of its formation and evolution.

The compilation of travel time tables of seismic body waves by Jeffreys and Bullen (8) was of fundamental importance for establishing the first radial average structure of Earth's interior.

Earth's seismic velocity structure has been constrained by a plethora of seismic observations that includes body waves, surface waves, and normal modes [START_REF] Dziewonski | [END_REF](10)(11)(12). However, replicating this progress on the Moon has been more difficult. For the Moon, strong scattering in the shallow parts of the lunar crust have largely limited seismic observations to the main P-and S-wave arrivals (13), and core reflections (14,15), from which averaged radial seismic profiles of crust, mantle, and core could be deduced (16). Relative to the four-station seismic array that oper-ated on the Moon, travel time inversions on Mars are complicated because accurate epicentral distance and origin time measurements of marsquakes are difficult to obtain with only one seismic station. Only event distance can be estimated from arrival time measurements of direct Pand S-waves (5). To obtain information on interior structure, observations of additional seismic phases such as planet-circling surface waves (17,18), normal modes (19,20), or reflected and refracted body waves that have interacted with the surface or internal boundaries are needed. As

we have yet to positively identify planet-circling waves or normal modes on Mars, we employed observations of P-and S-wave differential travel time measurements together with surface reflected body wave phases PP, PPP, SS, and SSS (Fig. S1) to jointly invert for both epicentral distance and interior structure. This strategy allows us to construct a quantitative model of the seismic velocity structure of the Martian mantle.

Marsquakes have been divided into two main categories based on their frequency content.

The low-frequency (LF) events have energy dominantly below 1 Hz, and high-frequency (HF) events have energy dominantly greater above 2.4 Hz (5,[START_REF] Clinton | [END_REF]. InSight has recorded far more HF events, and these events often lack clearly identifiable P-and S-wave arrivals and polarizations.

Here, we focused on eight high-quality LF events of the 43 recorded to date by InSight [START_REF]InSight Marsquake Service, Mars Seismic Catalogue[END_REF], labelled by mission Sol of occurrence and sub-labeled alphabetically for Sols with more than 1 event: S0167b, S0173a, S0185a, S0189a, S0235b, S0325a, S0407a, and S0484b (Fig. S9).

Each of these events has a high signal-to-noise ratio with identifiable P-and S-waves, which allows for epicentral distance estimation [START_REF] Clinton | [END_REF]. We determined the back azimuth of three events (S0173a, S0235b, and S0325a) from polarization and found that all were located in the Cerberus Fossae region [START_REF] Clinton | [END_REF]. These 8 events occurred in the distance range ∼25 • -75 • and had moment magnitudes between M W =3.0-4.0 [START_REF] Clinton | [END_REF].

An example broad-band filtered seismogram from marsquake event S0235b is shown in Fig. 1A and includes our picks of the direct P-and S-wave arrivals. Marsquake waveforms are characterized by codas produced by scattering (4,5), which, like their lunar seismogram counterparts, complicates identification of seismic body wave phases. Consequently, we do not identify seismic arrivals directly in the time series, but instead use narrow-band filtered timedomain envelopes (Fig. 1B) called filter banks [START_REF]Materials and methods are available as supplementary materials[END_REF]. Fig. 1B shows a zoom-in around the main P-and S-wave arrivals for event S0235b and our identification of surface reflected body wave phases PP, PPP, and SS that were made on the filter banks (bold black lines). We made our picks on the envelope peaks rather than the phase onsets because arrival onsets are more difficult to discern for the later arriving seismic phases. Filter banks for the 8 LF events considered here are shown in Fig. S9.

To provide an independent verification of the filter bank based picks, we used complementary methods for the identification of seismic phases. These include: (i) polarisation filtering and vespagrams and (ii) waveform matching [START_REF]Materials and methods are available as supplementary materials[END_REF]. Depite the independent nature of the three analysis methods, our picks (Table S4) are consistent within the overlapping uncertainties as illustrated in Fig. S9. To ensure that the arrivals identified in the envelopes are not due to noise, we compared the observed waveforms to the pulse shape of the direct P-and S-wave arrivals.

As we expected for multiple reflected body waves, the waveforms we observed are consistent with shifts in phase of 90 • and 180 • for PP/SS and PPP/SSS, respectively [START_REF]Materials and methods are available as supplementary materials[END_REF]. Based on the absence of surface waves and dispersion in the LF seismograms, we assumed the events occurred below the crust-mantle interface (5,[START_REF] Brinkman | [END_REF] at a depth of 50 km. We found that varying assigned event depth in the range 50-70 km produced little difference in inversion results. We note that for some events, possible depth phases are discernable (Fig. 1B), but are not considered in the inversion.

We used the travel time picks based on the visual inspection of filter banks (bold numbers in Table S4) as our preliminary differential travel time data set to simultaneously invert for epicentral distance (∆) and elastic seismic wave velocity structure (18,23). We modeled Mars as a spherically symmetric planet, although we acknowledge that it varies laterally in crustal thickness (24,25) and possibly in seismic properties (20,26). Based on global and regional synthetic seismograms we computed for one-and three-dimensional (1D, 3D) Mars models (23,26), we found generally small travel time differences between 1D and 3D models for the direct P-and S-wave arrivals (<3 s). Thus picking uncertainty (5-10 s) should be larger than what we expect for 3D effects. This might not be the case for the multiple reflected waves. For instance, 3D effects could potentially be more pronounced for the surface reflected phases as they travel in the more heterogeneous shallow structure beneath Elysium Planitia (1,4).

We carried out two independent inversions to determine interior structure using a seismic and a geophysical parameterization [START_REF]Materials and methods are available as supplementary materials[END_REF]. The seismic parameterization is based on P-and S-wave velocities in a layered planet model (23), whereas the geophysical parameterization (18,27) relies on an unified description of phase equilibria, seismic properties, and thermo-chemical parameters. We assumed a homogeneous bulk composition derived from the Martian meteorites (28) and computed P-and S-wave velocities as a function of temperature, composition, and pressure using a Gibbs free-energy minimization (29). Based on stagnant-lid models that include a crust enriched in heat-producing elements (30)(31)(32)(33), the Martian geotherm was parameterized using variable conductive crustal and lithospheric geotherms for which the underlying mantle can be assumed to be adiabatic. The exact nature of the crustal geotherm is less important here because we relied on a seismic parameterization of the crust [START_REF]Materials and methods are available as supplementary materials[END_REF]. We computed mantle adiabats (isentropes) self-consistently from the entropy of the lithology at the pressure and temperature of the bottom of the thermal lithosphere [START_REF]Materials and methods are available as supplementary materials[END_REF]. To solve the inverse problem (6), we employed a stochastic sampling algorithm (34) that produces a suite of models that fit the differential travel times within uncertainties and enables estimates of model parameter uncertainty and variability to be estimated. In total, we sampled 10 5 models of which ∼10 4 were retained to ensure model independence.

We inverted for radial P-and S-wave velocity and geothermal profiles (Fig. 2A,B) using both the geophysical (red and blue models) and seismic (gray-shaded models) parameterizations. The seismic velocity profiles we obtained from both inversions are very similar, with slightly wider posterior model parameter ranges in the case of the seismic parameterization.

This reflects the inclusion of mineral physics information into the geophysical parameterization. Based on the ray path geometry and sensitivity (Fig. 2C), we found that the differential travel time data constrained structure to 800 km depth, with sensitivity below ∼500 km provided primarily by the two events located at ∆>60 • (S0185a and S0167b). Below 800 km depth, the distribution of inverted profiles starts to resemble the prior P-and S-wave velocity distributions (gray shaded areas in Fig. S10), signaling loss of resolution. Comparison of prior (Fig. S13) and posterior (Fig. 2A,B) distributions of geothermal and P-and S-wave velocity profiles based on the geophysical parameterization confirms that the differential travel time data constrain upper mantle structure. For both parameterizations, S-wave velocity is better constrained than P-wave velocity, which reflects the more abundant number of S-wave picks relative to P-wave picks.

We found good agreement between the inverted epicentral distances for both parameterizations (Fig. S11), with slightly wider ∆ distributions for the seismic parameterization because of the broader velocity distributions. We focus on the implications of the geophysical parameterization as it provides a good fit to the observations (Fig. 2D,E).

We grouped the geotherms, including lithospheric geothermal gradients (Fig. 2B, inset), and seismic wave velocity profiles into families based on lithospheric thickness (Z lit ): Z lit <400 km (Fig. S13, green models), Z lit =400-500 km (blue models), and Z lit =500-600 km (red models).

We found that only models that group in a Z lit range of 400-600 km (blue and red models), corresponding to a thermal gradient between 1.6-2.5 K/km (mantle potential temperatures between 1600-1700 K), fit the data. These thermal gradients are consistent with those obtained on the basis of radar observations of the faint lithospheric flexure beneath the north polar cap (1.3-3.2 K/km) that constrains the thermal state of the present-day lithosphere (35).

Our crustal velocity models (Fig. 2A, insets) include an upper-crustal low-velocity layer with S-and P-wave velocities in the ranges 2-3 and 3.5-5 km/s, respectively, separated from the lower crust by an intra-crustal discontinuity around 5-10 km depth, which had been observed previously (4), and a larger discontinuity around 30-50 km that could be the crust-mantle interface. Lateral variations in crustal structure can result in differences between near-station structure and profiles inverted from our differential travel time data set, which averages structure across a wide geographic region, yet may not be entirely representative of the average structure of Mars. The crustal discontinuities are discussed in more detail in (36).

At greater depth, the lithospheric structure is characterised by constant negative S-and neutral-to-positive P-wave velocity gradients, respectively, to depths of 400-600 km. The Swave velocity decrease, an inherent feature of the models and consistent with the data, is followed by an increase to 800 km depth. This behavior defines a S-wave low-velocity zone (LVZ) in the Martian upper mantle. The LVZ results from the large thermal gradient across the lithosphere that arises when a relatively thick stagnant conductive thermal boundary layer sits on top of a convective mantle, as previously postulated (37). A possible manifestation of an LVZ is the presence of a weak S-wave shadow zone at epicentral distances of ∼40ĉirc-60ĉirc. There is no equivalent LVZ for compressional waves as P-wave velocity is less affected by temperature.

The LVZ on Mars results from the decrease in seismic velocity associated with increasing temperature with depth dominating over the opposing increase in seismic velocity associated with the increase in pressure with depth. Other mechanisms, such as fluid and melt in the asthenosphere (38) and redox effects (39) are also thought to play a role in shaping Earth's LVZ (40), but we lack the resolution to address these issues.

To provide an observational constraint on the strength and extent of the LVZ, we considered the alignment of the seismic traces (Fig. S9A) to qualitatively estimate the S-wave amplitude behavior with distance. From the alignment, we make the following observations: i) for ∆<40 • , both P-and S-wave arrivals are identifiable (S0235b, S0407a, S0484b, S0173a, S0189a, and S0325a); ii) for 40 • <∆<59 • , only a P-wave arrival is visible above the background noise on the vertical component, although there is a signal consistent with a low-amplitude S-wave (S0183a); and iii) for ∆>59 • , P-and S-wave arrivals are again distinctly visible (S0185a and S0167b). We interpreted the recovery in S-wave amplitude as the bottom of the LVZ as the increase in S-wave amplitude is commensurate with velocity switching to an increase with depth. This sequence provides tentative observational evidence for a weak S-wave shadow zone in the ∼40 • -59 • epicentral distance range.

We compared the InSight observations to reflectivity synthetics (6) for a range of models covering the inverted blue and red families (Fig. 3). The synthetic relative S-to-P-wave amplitude ratios illustrate the effect of the LVZ and predicts that the S-to-P-wave amplitude ratio falls off with increasing epicentral distance, but stabilizes where the velocity gradient begins to increase with depth. The drop in synthetic relative amplitudes between ∼53 • -59 • is compatible with the low S-wave amplitude event S0183a (not included in our travel time inversions), and allows us to revise the event's distance and the models obtained from the travel time inversion. The colored box (beige and brown) shows the allowed range (∼40 • -59 • ) based on alignment (Fig. S9A), whereas the brown box corresponds to the location (∼54 • -59 • ) for which synthetically-predicted amplitudes drop as observed for S0183a. Comparison with InSight observations (Fig. 3, gray boxes) indicates that the blue model family predicts an amplitude behaviour that is in line with the location of the observed amplitude drop (S0183a) and increase (S0185a). The red model family would appear to put S0183a (drop) in slight contradiction with the inferred location of S0185a (increase).

The weak S-wave shadow zone is based on an absent or low-amplitude S-wave from a single event (S0183a), which could potentially arise from the source radiation pattern. However, the relatively long (several minutes) coda that is observed throughout suggests that near-source scattering is substantial and that S-wave energy is emitted in all directions. As the ray paths of the direct and coda S-waves are similar, the simultaneous absence of both is best explained by geometrical spreading such as that produced by a shadow zone. Attenuation may also result in additional complexity if the quality factors (Q) are strongly depth dependent, such that the effective attenuation along the path is different for the surface reflected phases. Our initial observations from the direct phases suggest an effective Q≈200-300 to distances of 45 • (4, 5) and a similar value is expected for the surface reflected phases. Alternatively, if we assume that attenuation increases with depth, the more deeply diving S-waves would simply weaken and would therefore be unable to replicate the amplitude effects of a LVZ as evidenced by S0183a and S0185a.

The presence of the LVZ has important implications for the thermal evolution of the Martian mantle, the formation of the crust, and for the planet's surface heat flow. To identify plausible geophysical parameters that pertain to the dynamical evolution of Mars, we computed presentday thermal profiles using a parameterised stagnant-lid mantle and core convection model that simulates 4.5 Gyr of planetary evolution (32,(41)(42)(43). We explored all relevant geodynamic model parameters in computing the thermal evolution of Mars [START_REF]Materials and methods are available as supplementary materials[END_REF] and compared the resulting present-day lithospheric thermal profiles and crust and lithospheric thicknesses with the models constrained seismically (Fig. 2B).

We found that the seismic results are reproducible by parameterized convection models for a restricted range of geodynamic model parameters. In particular, the models point toward an initially relatively cold state (mantle potential temperature in the range 1630-1720 K) and a moderately sluggish mantle (with a reference viscosity η 0 = 10 20.2 -10 21.8 Pa•s). The bulk heat-producing element (HPE) content of the primitive mantle could be comparable to estimates based on the Martian meteorites (28,44) or 25-50% more enriched (45). Moreover, to match the seismic results requires a crust that is 13 to 20 times more enriched in HPEs than the primitive mantle, leading to an average surface heat flow of 14-29 mW/m 2 , a value that is generally higher than previously estimated (27,(31)(32)(33)(46)(47)(48). These enrichment levels also call into question models of crustal composition [START_REF]Materials and methods are available as supplementary materials[END_REF]49,[START_REF] Taylor | Planetary Crusts: Their Composition, Origin and Evolution[END_REF]. Models estimating crustal HPE from orbital gamma-ray near-surface mapping [START_REF] Taylor | [END_REF] predict HPE enrichments ≤12 times primitive mantle [START_REF] Taylor | Planetary Crusts: Their Composition, Origin and Evolution[END_REF]. Accordingly, deeper crust may consist on average of petrologically more evolved lithologies, such as those observed at Gale crater and in Martian brecciated meteorites (52) and in TES Surface Type 2 (53).

We related mantle thermal structure to surface heat flow using radial models with an average crustal thickness, neglecting regional variations that could vary across the surface of Mars, perhaps by a factor of two or more (33,47). Fortuitously, the area sampled by the seismic data, between Cerberus Fossae and the InSight landing site, has been predicted to be relatively homogeneous in terms of surface heat flow and close to the average planetary value (33). Thus, our approach appears justified given the large uncertainty in current surface heat flow predictions.

As of Sol 676, the SEIS instrument has operated at optimal noise conditions during the Martian Northern spring and summer between Sols 100 and 520 (1,5). Since then, the average wind has increased to a level that would make all but the three largest events (S0173a, S0235b, and S0325a) unobservable. Seismic monitoring conditions have improved again since Sol 780 (February 2021) and we expect around ∼10 additional P-S travel time observations during the extended mission. Meanwhile, our preliminary radial velocity model for the upper mantle of Mars will help guide and inform search for other seismic arrivals particularly core-related phases, possibly aid in aligning currently unlocated LF events to complement our data set, constrain a variety of geophysical and geochemical models, and models of planetary evolution. results. A.K. wrote the main paper with contributions from H.S., N.S., S.C., M.v.D., P.L., D.G., S.M.M., C.M., T.S., V.L., S.E.S., A.B., and R.R.M. Supplementary section 1 was contributed by amplitude ratio computed using a range of models covering the inverted blue and red model families (Fig. 2) that correspond to models with lithospheric thicknesses in the ranges 400-500 km and 500-600 km, respectively. The gray boxes indicate the observed S/P amplitude ratio for the events that have been determined from the time-domain envelopes (Fig. S9A).

Figure captions

The height of the gray boxes indicates the uncertainty on the observed amplitudes, whereas the location and width of the light gray boxes are based on the inverted epicentral distance distributions (Fig. S11), except for S0183a. For S0183a the width of the beige and brown boxes show the ranges based on the alignment (Fig. S9A) and the synthetically-predicted amplitudes, respectively.

Figures

SM 2 Seismic body wave picks SM 2.1 Seismic ray paths and body wave nomenclature

Ray paths for all body wave phases considered in this study are shown in Fig. S1. The main P-wave is the first compressional wave that arrives at a receiver, whereas PP and PPP phases are compressional waves that have reflected once and twice, respectively, from underneath the surface at a point midway (PP) and at points one third and two thirds (PPP) between source and receiver (for a source located at the surface). The same arguments apply to all S-wave phases.

While ray paths and bottoming depths are model dependent and therefore differ in detail, the resolution is generally common to all models (cf. Fig. 2C).

SM 2.2 Body wave picks using filter banks

All low frequency family marsquakes observed to date are of magnitude <4 and occur at distances >1200 km and are generally only observable because of the extremely low noise floor of the VBB seismometer (4). None of the events considered here surpass ambient noise outside the frequency range 0.1-4 Hz. Moreover, clear surface wave arrivals have not been observed so far, i.e., fall below the noise for all events, indicating that the detected low frequency marsquakes are unlikely to be shallow. All the events discussed here have amplitude estimates that are sufficiently above the background seismic noise [START_REF] Clinton | [END_REF]. Given the moderate distances of all these events (25 • -75 • ), surface reflections (PP, SS, PPP, SSS) are predicted to be separated by less than ∼50 s (P-wave phases) and ∼80 s (S-wave phases) from the main phase arrival; a time window within which depth phases (pP, sS) and receiver-side conversions and scattering are also found (4,[START_REF] Brinkman | [END_REF]. For these reasons, the phases discussed here (PP, SS, PPP, SSS) are not identified on seismograms, as would be usual practice on Earth, but are picked on smoothed time domain envelopes of the seismograms filtered to a set of narrow frequency bands (hereinafter filter banks). We used frequency bands that we specifically tailored for each event. We computed filter banks that are 1/2 octave wide on each side for frequencies between 1/16 Hz and 1 Hz, and determined the bands where the event is most visible and contains no glitch-or windrelated artefacts (54). For most of the events, this frequency band is 1/2.8-1 Hz. The arrivals show an increase in epicentral distance with S-P travel time (move-out). We only considered as robust body wave phase arrivals that persist across at least two frequency bands and show a consistent move-out between events. We observed main P-wave arrivals between 1/2 and 1 Hz for all events, most prominently on the vertical component. The PP and PPP body wave phases are mostly visible between 1/2 and 1/1.4 Hz. Main S-wave arrivals are, as we expected, most prominent on the horizontal components at longer periods, between 1/10 Hz and 1/2 Hz. We observed the multiples (SS, SSS) between 1/2.8 Hz and 1/2 Hz. Because arrival onsets are generally less clear for the later phases, we made our picks on the envelope peaks rather than phase onsets.

Our strategy, which relies on the use of filter banks, of which two examples are shown in Figs. S2-S3 (for each filter bank we also show the unfiltered time domain seismograms), can be summarised as follows:

1. We commence with the marsquake service-based picks (21) of the main P-and S-wave arrivals and use a filter bank to determine the highest frequency in which these phases are clearly observed. This frequency is controlled mainly by attenuation, since all events are small enough that their source spectrum can be assumed to be flat up to 1 Hz. For P-waves, this frequency is typically 1 Hz or 1.4 Hz, whereas for S-waves, it is 1/2 Hz or 1/1.4 Hz in ground displacement.

2. In bands at or directly below the main frequency, we identify and pick clear secondary peaks in the P-and S-wave coda that are consistent across at least two frequency bands.

Since the phase onset is usually masked by the coda of preceding phases, we pick the time of the maximum peak. For consistency, we repick the main phases on their peaks. The P-and S-wave peaks can therefore show slight differences to the values in the marsquake service catalog, which are picked on onsets; yet differences are within the uncertainties assigned to the picks.

3. From comparison of multiple events, we find that most main arrivals are followed by up to three peaks within 50-80 s that are most prominent at frequencies between 1/2.8 and 1/1.4 Hz. Of the three peaks, the two later-arriving peaks show a move-out with distance, whereas this is not the case for the first peak. Consequently, we interpret the first peak as a depth phase (either pP/sP or sS/pS).

4. Any arrival that coincides with a glitch (marked as gray vertical bands in the filter banks) is excluded (e.g., the SSS arrival for S0173a shown in Fig. S3 is picked on the North and not the East component).

While there are other peaks that at first appear likely, closer examination using the above procedure explains why such peaks are not retained in the final data set. We observed that arrival times of peaks can vary across a filter bank. This variation forms part of the uncertainty estimate on individual picks, in addition to the width of the entire pulse. The phase picks made using time-domain envelopes are listed in Table S1. By iterative application of this procedure, we converged on the set of picks shown in Fig. S9 and summarised in Table S4.

Table S1. Summary of body wave phase picks using filter banks. All travel times are relative to the main P-wave arrivals. The instrument response was removed from the data to obtain seismograms in velocity, and the three SEIS VBB components were rotated to the vertical (Z), North (N), and East (E) coordinate frame. We bandpass filtered the data between 0.3-0.9 Hz to remove long-period noise and avoid the 1 Hz tick-noise from cross-talk between the seismometer and the temperature sensors in the SEIS instrument package (55,56). We picked the P-wave arrival of each event on the vertical component (BHZ) based on the MQS-reported P-wave arrival time [START_REF] Clinton | [END_REF]. To obtain event backazimuth, we computed the average principal axis of P-wave particle motions, which is expressed as P = (P E , P N , P Z ), using a 5 s time window centered around the P wave. The back-azimuth of the event is then calculated using the following equation ( 54)

Event PP PPP S SS SSS (s) (s) (s) (s) (s) S0235b 17±5 36±5 167±5 195±5 - S0407a 22±5 38±5 168±5 194±5 - S0484b 22±5 40±5 172±5 - - S0173a - - 173±5 197±5 212±5 S0189a - - 191±5 213±5 239±5 S0325a 29±8 - 230±5 257±5 280±5 S0185a - - 356±5 397±5 - S0167b 55±5 - 418±5 466±5 502±5
BAZ = arctan P E P N + π (1)
We then rotated the horizontal components (BHN and BHE) of the event to radial (BHR) and transverse (BHT) components based on the P-wave derived back-azimuth.

The waveforms of Marsquakes are often characterized by strong coda due to scattering (4), therefore complicating the identifications of body waves. Polarization filtering is a method used to enhance linearly polarized waves such as body waves while suppressing other non-polarized waves including surface waves, background noise, or scattered waves. A polarization filter was used to detect core-reflected body waves in the Apollo seismic data (14). Here we followed the original method by (57) to design a polarization filter for improving the signal-to-noise ratio (SNR) of Martian body waves. First, we computed the covariance matrix of the threecomponent data (R, T, Z) in a given time window (5 s) around t 0 , and the eigenvectors and eigenvalues of this matrix. The rectilinearity at time t 0 is defined as the following equation

RL(t 0 ) = 1 - λ 2 λ 1 n J , (2)
where λ 1 and λ 2 are the largest and second-largest eigenvalues respectively and n and J are empirical exponents. The eigenvector associated with the principal axis is Ē = (e 1 , e 2 , e 3 ) with respect to the RTZ coordinate system. The weights on each direction at time t 0 is given by

D i (t 0 ) = e K i , (3)
where i = 1, 2, 3 (R, T, Z) and K is an empirical exponent. We chose these empirical values for the exponents: n=0.5, J=1, and K=2. The polarization filter for each component is defined as the product of rectilinearity and the corresponding direction weight. The filtered three-component seismograms at any time t are given by

R f (t) = R(t) • RL(t) • D R (t) (4)
T f (t) = T (t) • RL(t) • D T (t) (5)
Z f (t) = Z(t) • RL(t) • D Z (t) (6)
We applied the polarization filter to all three components (BHZ, BHR, BHT), followed by an amplitude envelope of the data. We then repicked the P-wave arrival on the vertical component, and the S-wave arrival on the transverse component, with guidance from the MQS picks [START_REF] Clinton | [END_REF].

Here we chose to pick P-and S-wave arrivals at the maximum amplitude of the envelope. We avoided known glitches identified by MQS by excluding glitched time windows as potential picks (Fig. S9, (54)). A subsequent event distance alignment was made based upon our picks for the main P-and S-wave arrivals, and the corresponding waveforms are ordered by S-P time (Fig. S4).

SM 2.3.2 Vespagram analysis

To detect seismic phases that may arrive coherently after the direct P-and S-waves, we applied a velocity spectral analysis (vespagram) approach to stack the envelopes across a range of slownesses for common body waves. The vespagram approach was used to identify multiple reflected phases (e.g., SS and PP) and core phases (e.g., ScS and PcP) on Earth using seismic arrays (58,59) and has often been used to study low amplitude phases that arrive near a larger amplitude parent, (e.g., the SS and PP precursors ( 60)). The vespagram approach improves the signal to noise of coherent phases with similar move-outs and also distinguishes phases with distinct slowness (e.g., core phases). Here, we adapted the vespagram approach to a source array configuration on Mars, in which we leverage the relative timings of seismic waves arriving from a range of epicentral distances at a single station, in lieu of the more traditional station array method. The epicentral distance of each event was derived from the S-P time based on a model constructed only for the purpose of alignment (as described in the main manuscript), and the average distance was then used as the reference distance for vespagram analysis. Note that no vespagram predictions for events S0185a and S0167b are included because the larger distances of these events compared to the rest would violate the plane wave approximation used in the vespagram analysis. Glitches were removed by zeroing the amplitudes in the time windows of glitch arrivals. We created vespagrams for the polarization-filtered LF envelopes on two components: (1) BHZ component aligned on the P-wave arrival; (2) BHT component aligned on the S-wave arrival. The goal is to identify any coherent body waves that come in with consistent timing and slownesses after the initial P-and S-wave.

The BHZ component vespagram (Fig. S5a) reveals focused P-wave energy at zero relative slowness and a subsequent packet of coherent energy with slightly larger relative slowness (0.8±0.3 s/deg, 23±5 s in Fig. S5c). We interpret this second arrival as the seismic phase PP. A third more ambiguous arrival is tentatively interpreted as the PPP phase (1.2±0.4 s/deg, 42±5 s in Fig. S5c), although there are a series of arrivals near PPP, which could result from structural complexity near the bounce point. The vespagram for the BHT component has low energy at the expected arrival times of P-wave energy, suggesting that the back-azimuth used to rotate the records was appropriate (Fig. S5b). The BHT component vespagrams are aligned on the S-wave pick, so the S-wave energy focuses at zero slowness and time. A subsequent packet of energy at larger slowness (0.5±0.2 s/deg, 25±5 s in Fig. S5d) is interpreted as the SS phase.

There is a series of more complex packets of energy following the SS phase, which we interpret as the SSS phase (1.0±0.2 s/deg, 43±5 s in Fig. S5d).

We used our vespagram results to guide the identification of PP and SS arrivals in the individual polarization-filtered event waveforms, as well as, where warranted, PPP and SSS arrivals.

The observed slowness of each phase from the vespagram and reference distance (31.59 • ) was used to compute the predicted arrival time of PP and SS, as well as PPP and SSS for individual events (Fig. S9). The predicted arrival times then guided picks on the individual body waves in each event, with selection based upon the most coherent and high amplitude peak within ±10 s of the vespagram-derived predicted travel times. These picked times are summarized in Table S2.

SM 2.4 Body wave picks using waveform matching

Surface-reflected body wave phases (PP, PPP, SS, SSS) go through a caustic after each reflection at the surface and two-way transmission through the crust. For SS (and PP), this leads to a pulse distortion corresponding to a Hilbert transform of the emitted S-wave phase (61,62), in addition to a change in polarity related to the reflection and a further possible polarity change at the source, when the differences between the take-off angles are large between S, SS, and SSS. As shown in Fig. S6, this phase distortion and polarity change is clearly identified in the events with the best signal-to-noise ratio, such as S0235b, S0407a, S0484b, and S0173a, for which an alignment of the Hilbert transformed S-and SS-waveforms can be performed. The alignment remains stable for various filter settings (see below) and provides an independent means to the filter bank method (section SM 2.2) described in the main text for identifying phase arrivals. A similar alignment can be performed for PP versus P (63). This results in SS-P and PP-P differential travel times with a root-mean-square (rms) error ≤4 s in comparison to the measurements made using the filter bank method and ≤2 s with the vespagram method (section SM 2.3).

For SS and S, we search for these phases by waveform matching of the S-wave on the transverse component. This is done by computing, for each time window, the best least-squares fit of 20 s of signal with the first 15 s of the S-wave (the S waveform therefore starts 5 s before the MQS reported arrival time [START_REF] Clinton | [END_REF]). Prior to analysis, all seismograms have been deglitched (4, 54) and the instrument response has been removed to produce final seismograms in ground displacement (with a high-pass cutoff frequency of 0.1 Hz). The least-squares fit is made using the displacement seismograms that have been further band-pass filtered from 0.3 Hz to f c , where f c assumes values between 0.4 Hz and 0.9 Hz, in steps of 0.1 Hz. For each time window and f c value, the fit is computed by scaling the S-wave amplitude, which provides an estimation of the reflection coefficient of the phase with respect to S. The resultant signal is either the Hilbert transform of the S-wave coda (for SS) or the S-wave coda without the Hilbert transform (for SSS, crustal reflections, or sS). The same alignment is also made for SSS (with respect to S)

and PPP (with respect to P) but without the Hilbert transform due to the double π/2 phase shift.

These signals are generally more affected by the coda of the primary wave and have a larger probability than the SS to have a polarity change related to the take-off angle.

Fig. S7 shows results for S0235b and includes SS-and SSS-wave arrival times (relative to S) obtained by the filter bank and vespagram methods, respectively. Fig. S7A,C,E show the case when the Hilbert transform of the coda is used for the waveform matching, while Fig. S7B,D,F show the waveform match without application of the Hilbert transform. SS phases are suggested when the three following conditions hold: (i) decrease in variance residual (Fig. S7A);

(ii) small variations of the reflection coefficient for the different low pass filters (Fig. S7C), as a consequence of the low frequency dependence of the crustal reflection coefficients; and (iii)

negative phase/S ratio on the Hilbert-transformed signal (Fig. S7E). SSS phases follow the same criteria, but on the non-Hilbert transformed signal (Fig. S7B,D,F). If we have polarity inversion of SS or SSS due to the source geometry, nearby secondary minima must be considered and included in the error assessment. For S0235b, the best matches are achieved for the SS phase with a 22.5 s delay with respect to the S, which is close to the SS pick proposed in the vespagram analysis (22.8 s, see Table S2). For such a delay, the variance residual is 35%, while the relative amplitude of the SS phase relative to the S (i.e., scaling coefficient) is found to be -0.2. Uncertainties in the polarity are also integrated into the error bar of ±2.5 s (green arrow), which includes matches with comparable amplitudes in the range 0.15-0.20. For the SSS, we observe the best match at about 38 s with a scaling stability of 50% and a positive amplitude of 0.15. A negative polarity, however, is found at 37 s and leads us to propose that SSS-S is 38 s.

The corresponding SSS-wave arrival proposed using the vespagram method is 34.5 s. Note also that in Fig. S7, at about 10 s after the S-wave arrival, a reflected S-wave is identified in the S coda with a negative reflection coefficient of about -0.15. This coefficient can be obtained by a discontinuous increase in S-wave velocity from 1.9 km/s to 2.6 km/s for the 9-km discontinuity proposed in ( 4). An equivalent reverberation is found for SS. The waveform matches are shown in Fig. S8 for the SS-wave arrival of S0235b. This confirms the waveform fit after Hilbert transformation and also includes the inner-crustal reflected phase of the SS about 10 s after the main SS-wave arrival. The same analysis was applied to SS and SSS for the other events and to PP and PPP (63).

Results are summarized in Table S3 for the four events for which a good or reasonable waveform match is obtained. The mean squared SS-S differences for all common picks compared with those made using the filter banks (Table S1) and the vespagrams (Table S2) are <4 s and <2 s, respectively. For SSS-S, differences amount to 0.5 s for S0173 (filter banks). Relative to the vespagrams, the rms for all 4 events and picks is <2.5 s. For PP-P and PPP-P, rms differences are also <4 s in the case of the filter banks (for S0235b, S0407a, and S0484b) and <2.5 s in the case of the vespagrams for all events and picks. These differences suggest that the 5 s error assumed in the inversion is likely to be on the conservative side.

Table S3. Summary of the secondary phase delays (PP-P, PPP-P, SS-S, SSS-S) for the four events with the largest signal-to-noise ratio using waveform matching. Both the P-wave and the S-P differential times are from the MQS catalog (7, 21) except for S0484b, for which the start of the P-wave is 2020-04-07 08:52.34 UTC while the start of the S-wave is based on the MQS picks. D show the relative differences between the scaling coefficients, when computed for the 6 different bandwidth filters (see Fig. S8). The scaling coefficient can be interpreted as the relative amplitude ratio of SS or SSS phase with respect to S. Panels E and F indicate average values of the scaling coefficient, again for the 6 bandwidth filters and with respect to the S-wave packet. For optimally matching a phase, the scaling coefficient can be interpreted as the relative amplitude ratio of this phase with respect to S. For S0235b, the best match is obtained for 22.5±2.5 s, as indicated by the green arrow (panels A, C, and E). For the SSS phase of S0235b, we propose an arrival at 38±3 s, as indicated by the green arrow in panels B, D, and F. Time after S refers to the MQS-picked S-wave arrival time. In panel A, the S-wave (gray) is compared to the Hilbert transform of the S-wave, delayed by 22.5 s and scaled by -0.2 (blue). A delay of 22.5 s provides a good match to the waveforms between 185 s and 210 s, suggesting that all these phases have gone through a surface reflection and a caustic. The window used for the least-squares waveform fitting is given by the two vertical red dashed lines, while the S waveform is given by the two green dashed lines. The green solid line is the MQS S-wave arrival time, and the red solid line is the proposed SS-wave arrival time, using the inverted SS-S delay with respect to S (Fig. S6). The waveform matching remains good for all bandwidths from 0.3 Hz to 0.9 Hz. In panel B, the same waveform matching is attempted without performing the Hilbert transform of the S-wave pulse. Even if some parts can be matched with a positive reflection coefficient, the matching is not achieved at longer periods beyond the 20 s window, illustrating that this phase cannot be explained with a Hilbert transform.

Event PP-P PPP-P SS-S SSS-S S-P

(s) (s) (s) (s) (s) S0235b

SM 3 Stochastic inversion

To solve the inverse problem d = g(m), where d is a data vector consisting of observations and g is an operator that maps from the model space into the data space, we employ a Bayesian approach as outlined in ( 34)

σ(m) = k • h(m)L(m), (7)
where h(m) is the prior model parameter probability distribution, i.e., the information on model parameters procured independently of data, L(m) is the likelihood function, which measures the misfit between observed and predicted data, k is a normalization constant, and σ(m) is the posterior model parameter distribution. σ(m) represents the solution to the inverse problem above. L(m) is determined from the observed data, data uncertainties, and the manner in which the latter are used in modeling data noise (to be described in the following).

The Metropolis algorithm is employed to sample the posterior distribution (Eq. 7) in the model space (34). This algorithm, which samples the model space in a random fashion, is an importance sampling algorithm, i.e., it ensures that models that fit data well and are simultaneously consistent with prior information are sampled more frequently. The Metropolis algorithm samples the model space with a sampling density that is proportional to the (target) posterior probability density and thus ensures that low-probability areas are sampled less intensively.

This is an important feature of any algorithm that wishes to randomly sample high-dimensional model spaces where the probability density over large proportions of the volume are near-zero.

SM 3.1 Sampling the posterior

Assuming that the errors can be modeled using an exponential probability density (L 1 -norm), the likelihood function can be written as

L(m) ∝ exp - i |d i obs -d i cal (m)| σ i (8)
where i runs over the number of travel time picks (27), d obs and d cal (m) denote observed and calculated data, respectively, and σ differential travel time uncertainty (bold numbers in Table S4). The reason we resort to the L 1 -norm is because the posterior probability density function is less prone to be influenced by outliers in the data set [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF][START_REF] Khan | [END_REF].

SM 3.2 Seismic model parameterization

In order to provide an independent check of the inversion results presented in the main manuscript, we inverted the differential travel time data using a "standard" seismic parameterization, i.e., a layered model with variable P-and S-wave velocities. In contrast to the methodology applied in the main part of the manuscript, this approach provides increased flexibility in determining the velocity models, since it allows inversion for both P-and S-wave velocities independently of any thermodynamic or petrologic property. Yet, herein also lies its disadvantage. To retrieve the general features of the velocity structure, we consider a simplified model that is divided into two crustal and six fixed mantle layers. The crustal layers are variable in thickness, whereas the mantle layers are of fixed thickness. Each layer is parameterised in terms of variable P-and S-wave velocities. Crustal velocities are assumed to increase as a function of depth. The a priori model parameter information, which consists of constraints on absolute Pand S-wave velocities, is summarized in Table S5. As part of the inversion, we also determine epicentral distance.

To infer velocity structure, we invert the differential observed body wave phase picks based on the time-domaine envelope picks that are tabulated in Table S4 (bold numbers) using the same Bayesian framework described in section SM 3 and employed in the inversion in the main part of the manuscript. To sample the posterior distribution (Eq. 7), we employ the Metropolis-Hastings algorithm, assuming that error can be modeled using an exponential probability density (Eq. 8). Sampling is performed using ten independent chains with a total length of 10 4

iterations, characterized by identical initial models but different randomly chosen initial perturbation. Finally, the inverted radial S-and P-wave velocity profiles for the seismic model parameterization are shown in Fig. S10 and are based on the 3×10 4 best-fitting candidates.

The current set of differential travel time data allows us to constrain structure down to 800 km depth. In agreement with the features described in the main text, the S-wave velocity models retrieved here show slightly negative gradients toward the bottom of the thermal lithosphere around ∼400 km depth. Sampled epicentral distances are shown in Fig. S11 and are seen to be in good agreement with those obtained from the geophysical parameterization (see next section). Fig. S11. Inverted epicentral distances for the 8 events considered here using the seismic (light shaded bins) and the geophysical parameterization (dark shaded bin).

SM 3.3 Geophysical model parameterization

The seismic velocity and thermal models discussed in the main manuscript and shown in Fig. 2A, B are based on the geophysical parameterization that we describe here. To construct models of the internal structure of Mars, we follow the approach outlined in (27), where mantle seismic properties (P-and S-wave speeds) are computed using petrologic phase-equilibrium computations in the NCFMAS model chemical system comprising the oxides CaO-FeO-MgO-

Al 2 O 3 -SiO 2 -Na 2 O.
For this, we employ Gibbs free-energy minimization (29), which computes stable mantle mineralogy and physical properties as a function of temperature, composition, and pressure based on the thermodynamic formulation and parameters of (66,67). The virtue of this parameterization is that it relies on a unified description of phase equilibria, seismic properties, and thermo-chemical parameters. Major sources of uncertainty in the thermodynamic calculations are absence of experimental constraints on the parameters relevant for the thermodynamic formalism and parameterization of Stixrude and Lithgow-Bertelloni (67). ( 68) estimated the accuracy of the elastic moduli and density to be ∼0.5 and ∼1-2%, respectively.

Tests conducted as part of this study have shown that changing bulk mantle composition has little impact on seismic structure (see also (48)), since differences among the major model Martian bulk compositions of e.g., Dreibus and Wänke (44), Lodders and Fegley (69), Sanloup

, Khan and Connolly (71), Taylor (28) are relatively small. Accordingly, we fixed the mantle composition to that of (28). For the crust, we rely on a "standard" seismic parameterization (Pand S-wave velocity) and consider a 3-layer crust with variable depth nodes and variable P-to-S-wave velocity scaling (Fig. 2A, insets).

Martian geotherms are parameterized by a conductive crust and lithosphere underneath which the mantle is assumed to be adiabatic (Fig. S12). The conductive part is defined by two linear thermal gradients (blue lines in Fig. S12): the first is between the surface and the bottom of the crust-mantle interface (of variable thickness Z and temperatures T ) and the sec-ond between the crust-mantle interface and the base of the lithosphere (of variable thickness Z lit and temperature T lit ), respectively. Because of enrichment in heat-producing elements, a moderate curvature may develop in the crustal thermal gradient. However, the exact nature of the crustal geotherm is unimportant for the inversion since we resorted to a seismic parameterization of the crust as described above. This assumption of crustal seismic structure decoupled from temperature originates from the fact that the crustal seismic structure of Mars may be dominated by various distinct processes that do not directly or primarily relate to temperature such as alteration, porosity, compositional and structural heterogeneity, and fracturing/damage (48).

The thermal structure of the mantle is assumed to be adiabatic and the adiabats (isentropes) are computed self-consistently from the entropy of the lithology at the pressure and temperature equivalent of the bottom of the lithosphere (29). Model parameter values and ranges are summarised in Table S6.

T surf T lit Z lit

Temperature Depth

Lithosphere Asthenosphere -C). Insets in panels (A), (B), and (C) show the distribution of sampled prior lithospheric geothermal gradients (dT/dz) and crustal S-and Pwave velocity structure, respectively. Profiles are colour-coded according to lithospheric thickness: 300-400 km (green), 400-500 km (blue), and 500-600 km (red). The lithospheric thermal gradient is determined from the temperature at the crust-mantle interface and the bottom of the lithosphere and the difference in depth between the two points. no significant changes in the characteristics of the posterior probability densities. To obtain independent samples, we introduced an "elapse time" (number of iterations) between retention of samples, which was found to be around 10 by analysing the autocorrelation function of the fluctuations of the likelihood function. In all 10 5 models were sampled with an acceptance rate of 32%. Inverted ("relocated") epicentral distances for the eight events considered in the main text are shown in Fig. S11. In comparison to the distances determined using the seismic parameterization, the distributions here are narrower, which reflects the fact that the velocity distributions are narrower because of added mineral physics information in the present parameterization. Finally, we show a zoomed-in version of the data fit shown in the main text (Fig. 2D,E) in Fig. S14.

T' Z' Moho

Differential travel time [s]

SM 4 Martian seismic velocity models and amplitude computations

Fig. S15 shows two seismic velocity models that are characteristic of the inverted families of models shown in Fig. 2A in the main manuscript. The blue and red models correspond to models with lithospheric thicknesses in the range 400-500 km and 500-600 km, respectively.

The models are available in digital format.

SM 4.1 Computing synthetic waveforms and amplitudes

The seismic S-/P-wave amplitudes show in Fig. 3 are obtained from reflectivity synthetics (72) calculated for a 50 km deep double-couple source (strike=0 • , dip=45 • , rake=90 • ) in a range of models that covers the blue and red inverted models shown in Fig. 2. The source mechanism is obtained from waveform inversion of the event S0235b [START_REF] Brinkman | [END_REF]. This implies the same source mechanism for all events, which is a working assumption that will be tested in a forthcoming study. The S-and P-wave attenuation models employed for this purpose are shown in Fig. S16 and have been constructed so as to be in overall agreement with InSight observations that indicate an effective Q around 300 based on spectral fits of the closer events and a higher effective Q for the more distant events (4,5). Synthetic data are in acceleration and band-pass filtered from S4).

Since a vast majority of Martian events do not exhibit clear P-and S-wave onsets in the time series, we based our phase picking practice on identifying seismic phases using the timedomain envelope peaks. This practice, by definition, causes some level of uncertainty when determining the observed S-to-P amplitude ratios (Fig. 3) depending on whether the amplitudes are measured where the envelope starts to emerge, or at the peak of the envelope. In order to compare observational amplitude ratios with those from 1D synthetics, we consider the spread in S-to-P values, based on the minimum and maximum P-and S-wave envelope amplitudes within the uncertainty windows shown in Fig. S9, respectively, as representative of measure- To establish a connection to the thermal history of the planet, we performed a set of geodynamic calculations using a parameterized stagnant-lid mantle convection model coupled to a convective liquid core that simulates 4.5 Gyr of thermo-chemical evolution (32), in which we updated the melting model that feeds into the crustal production using melting curves that are more representative of a Martian composition (74). We conducted an exploration of all relevant geodynamic model parameters (Table S7) in computing the thermal evolution of Mars.

These parameters include initial mantle temperature beneath the lithosphere (T m0 ), initial core temperature at the core-mantle boundary (T c0 ), mantle rheology (reference viscosity η 0 ), effective activation energy (E * ) and volume (V * ), crustal enrichment factor in heat-producing elements (HPE) relative to the primitive mantle (Λ), and initial bulk heat-producing element contents (C U , C T h , and C K ). Parameter search ranges are given in Table S7. We considered only initial core temperatures hotter than mantle temperatures to allow for the generation of an early Martian dynamo (75). Within the explored parameter space, we retained only models for which the mantle remained in a convective state (supercritical Rayleigh number) as suggested by the observed traces of recent volcanism on the surface of Mars (76,77). Note also that the upper bound for the range considered for the crustal enrichment factor (Table S7) could have been larger i.e., 30 instead of 20. However, we found no thermal evolution run compatible with crustal enrichment factors greater than about 20. Beyond this threshold, the HPE content in the crust would be so large that the mantle would be entirely depleted in HPEs, or would require a negative mantle HPE concentration, which is unphysical. Consequently, we set the upper bound for the crustal enrichment to 20. Each model parameter was randomly sampled within the prior ranges and checked to ensure that sampling was sufficiently dense for meaningful statistics.

Fig. S17c), moderately thick present-day crustal thicknesses (D cr =40-90 km, fig. S17g), moderate to large HPE content (C U > 7 ppb, Fig. S17e) and a present-day surface heat flow > 12 mW/m 2 (Fig. S17h). Sampled values of other quantities (activation energy and volume, crustal enrichment factor relative to the primitive mantle and initial thermal state) for this family are similar to the prior range, indicating that these quantities, strictly speaking, cannot be constrained even if the distributions favor, in a statistical sense, moderate E * values (150-350 kJ/mol, Fig. S17a), relatively small values for V * (<∼ 4 cm 3 /mol, Fig. S17b), and moderate to large HPE crustal enrichment (Λ ∼10-17, Fig. S17f).

The blue family generally yields considerably more constrained thermo-chemical histories and governing parameters. These profiles are associated with a large to strong bulk HPE content (C U > 17 ppb, Fig. S17m), a moderate to large crustal enrichment of HPE (Λ > 8-20, Fig. S17n), resulting in moderate to large values of present-day surface heat flow (22-48 mW/m 2 , Fig. S17p). The present-day crustal thickness associated with this family is relatively large (D cr =70-90 km, Fig. S17o). Similarly to what was obtained for the red family, the corresponding values for mantle rheological parameters span most of the explored range. However, the obtained statistical distributions favor even narrower ranges than those corresponding to the red family. Most likely values are E * = 150 -280 kJ/mol for the effective activation energy (Fig. S17i), V * = 0 -3 cm 3 /mol for the effective activation volume (Fig. S17j), and 10 21.2 -10 22 Pa s for the reference viscosity (Fig. S17k). However, the associated initial thermal state remains poorly constrained (Fig. S17l).

The most constrained quantities and their corresponding ranges for the two families of lithospheric geotherms are listed in Table S8. The difficulty in constraining the thermal state and mantle rheology reflects the loose constraints on the present-day geothermal lithosphere that are incapable of reducing the trade-offs between temperature and rheology for the thermal evolution of Mars (hotter initial temperatures and lower intrinsic mantle viscosity can both enhance of the lithosphere is not rigid, and can therefore be penetrated by hot uppwellings, which in conjunction with other processes (e.g., the presence of water), would trigger decompression melting at depths less than 500 km (79), and allow to trigger recent volcanic events as observed at the surface (76,77). Applying the same constraints on crustal thickness to the blue family of geotherms results in no thermo-chemically compatible solutions because the associated crustal thickness range is above 70 km (Fig. S17o).

Overall, combining considerations on present-day lithospheric geothermal gradient, lithospheric thickness, and crust-mantle interface depth obtained from seismic observations point towards an initially relatively cold thermal state, a relatively sluggish mantle whose bulk HPE content challenges the most accepted current estimates that are more enriched in HPEs. In addition, matching the seismically induced lithospheric and crustal constraints requires a crust that is 13 to 20 times more enriched in HPEs than the primitive mantle, and a corresponding average surface heat flow of 14-29 mW/m 2 , which is at the higher end or above pre-mission estimates.

SM 5.2 Compositional models of the Martian crust

The level of crustal HPE enrichment implied by our thermal modeling also has implications for constraining the chemical composition of the Martian crust. Models for Martian crustal composition all rely on near-surface measurements (microns to meters), such as chemical deconvolution of surface mineralogy determined by thermal emission (TES) and visible to near-infrared (VNIR) spectroscopy (49,53,80), direct chemical mapping of the near-surface by gamma-ray spectroscopy (GRS) (49,[START_REF] Taylor | Planetary Crusts: Their Composition, Origin and Evolution[END_REF] or average compositions of surficial regolith and sediment [START_REF] Taylor | Planetary Crusts: Their Composition, Origin and Evolution[END_REF].

In the absence of direct observations of lower crustal compositions, these near-surface mea-

Fig. S17. Distributions of physical quantities associated with the parameterized thermo-chemical evolution of Mars and corresponding to two different families of lithospheric geotherms obtained in the inversion of the seismic data: (a-h) relatively hot and thick lithosphere (red models in Fig. 2B) and (i-p) relatively hot and moderately thick lithosphere (blue models in Fig. 2B). The initial heat-producing elements (HPE) concentration is randomly sampled. The prior ranges are displayed in light blue. The red and the blue set corresponds to 10 4 and to 3307 randomly sampled thermo-chemical histories, respectively, among a prior set (transparent light blue) composed of 6 • 10 5 model evolutions. See text for further details. surements necessarily are assumed to be representative of the entire Maritan crust, leading to a major uncertainty.

The most comprehensive crustal model to date, based on Martian soils and GRS, corrected for near-surface volatiles and a small meteoritic component, proposes HPE abundances of K=3740 ppm, Th=0.70 ppm and U=0.18 ppm [START_REF] Taylor | Planetary Crusts: Their Composition, Origin and Evolution[END_REF]. For the primitive mantle HPE values in Table S7, based on (44), a crustal HPE enrichment factor of 12 is obtained. This is likely an upper limit since other models of Martian primitive mantle, based on Martian meteorites (28,45) or mixtures of classes of chondritic meteorites (69,70) predict even higher HPE abundances and would thus lead to lower crustal HPE enrichment factors.

Why do GRS-based crustal composition models appear to underestimate HPE abundances?

One possibility is that the Martian crust has internally differentiated into petrologically distinct upper and lower crust, analogous to the terrestrial crust. This would also be consistent with our seismic evidence for crustal layering (Fig. 2a; also see ref. [START_REF] Brinkman | [END_REF]). However, on Earth, both continental and oceanic crustal differentiation result in an upward enrichment of incompatible HPE, leaving the upper crust more, rather than less, enriched. In the case of continental crust, this is due largely to intracrustal partial melting processes giving rise to granitic and other evolved magmatic rocks in the upper crust [START_REF] Taylor | The continental crust: Its composition and evolution[END_REF] and in the much thinner oceanic crust, results from fractionation of primary mantle melts into MORB-like melts preserved in sheeted dikes and pillow lavas (Oceanic Crust Layer 2) overlying cumulate-rich gabbroic rocks (Layer 3) with lower HPE abundances (White and Klein, 2014). If analogous intracrustal differentiation processes were to explain low HPE in the Martian upper crust, they would have to be fundamentally different to result in downward enrichments in incompatible elements.

There are at least three (non-mutually exclusive) possibilities to explain why GRS-based crustal composition models may underestimate HPE abundances. One possibility is that the crust has been internally differentiated into petrologically distinct upper and lower crust, analo-gous to terrestrial continental crust. However, intracrustal partial melting would lead to upward enrichment of the incompatible HPE [START_REF] Taylor | The continental crust: Its composition and evolution[END_REF], leaving the upper crust more, rather than less, enriched. A second possibility is that the crust has been underplated by later, more HPE-enriched magmas. However, fractionation of ascending mantle melts, leading to underplating, more likely results in accumulation of ultramafic to mafic cumulates at the base of the crust which would be more depleted in HPE [START_REF] Mccubbin | [END_REF]. We consider the more likely possibility is that the uppermost crust has preferentially incorporated components of HPE-depleted mafic material (e.g., SNC basaltic meteorite compositions) during younger magmatic re-surfacing that is known to have occurred, thus effectively "diluting" near the surface an older, more HPE-enriched bulk crust.

Recent observations from rover in-situ measurements, Martian meteorites and orbital spectroscopy have identified considerable magmatic diversity on Mars including numerous occurrences of petrologically evolved igneous rocks that could form part of an older HPE-enriched crust (see reviews in (52,83)). These include ancient (Noachian to early Hesperian) igneous rocks in Gale crater that may be part alkaline igneous differentiation series or mantle melts, Noachian evolved igneous clasts in Martian brecciated meteorites (NWA 7034 and its numerous pairs) and regions of possible felsic materials in the southern highlands and near the base of Valles Marineris. In addition, it has long been recognized on the basis of orbital thermal emission spectroscopy (TES) that much of the younger northern lowland surface is characterized by "andesite"-like compositions termed TES Surface Type 2 (84), possessing K abundances in excess of 1% (53). Whether such compositions reflect igneous precursors or relatively young surficial alteration processes remains controversial (53) but, if reflecting igneous processes, could also represent a significant reservoir of felsic crust.

  S.C., D.G., A.K., N.S., Q.H., P.L., D.K., C.C., J.F.C., A.C.D., P.M.D., M.D., F.K., V.L., R.R.M., W.T.P., M.P., J.-R.S., and R.W.-S. Supplementary section 2 was contributed by A.K. and A.C.D. Supplementary section 3 was contributed by A.C.D., A.K., N.S., S.C., Q.H., and D.K. Supplementary sections 4-5 were contributed by H.S, S.M.M., and A.B. The InSight mission is managed by W.B.B., S.E.S., and M.P.P. Competing interests: The authors declare no competing interests. Data and materials availability: The InSight seismic event catalogue (comprising all events including phase picks, until July 2020) and waveform data are available from the IRIS-DMC and SEIS-InSight data portal https://www.seis-insight.eu/en/science. The catalogue and waveform data have the identifiers https://doi.org/10.12686/a9 and https://doi.org/10.18715/SEIS.INSIGHT.XB_2016, respectively. All In-Sight seismic data used in this study are also available from the IPGP data center: http: //dx.doi.org/10.18715/SEIS.INSIGHT.XB_2016.

Fig. 1 .

 1 Fig. 1. Marsquake recorded by the InSight seismometer. (A) Three-component broad-band filtered (1.5-8 s) seismogram of event S0235b. The main direct P-and S-wave arrivals are indicated by the vertical arrows. (B) Close-up view of the P-wave (left) and S-wave (right) arrivals, showing the direct and surface-reflected body waves (PP, PPP, and SS). Because arrival onsets are less palpable for the surface-reflected phases, arrivals are picked on the peak rather than the onset (arrows). The arrival of a possible depth phase is also marked. Bold black lines indicate the time-domain envelopes that are used for picking phase arrivals. The P-and S-wave envelopes are filtered at 2 s and 1.4 s, respectively, and are 0.5 octave wide on either side. Horizontal gray bars indicate the measurement uncertainty on the picks.

Fig. 2 .

 2 Fig.2. Summary of Mars's upper mantle structure. (A-B) Inverted S-and P-wave velocity and geothermal profiles. Coloured (red and blue) and gray-shaded models are obtained from the geophysical and seismic inversion, respectively. Insets in panels (A) and (B) show the distribution of sampled crustal S-and P-wave velocity structure and lithospheric geothermal gradients (dT/dz), respectively. Profiles are colour-coded according to lithospheric thickness: 400-500 km (blue) and 500-600 km (red). The lithospheric thermal gradient is determined from the temperature at the crust-mantle interface and the bottom of the lithosphere and the difference in depth between the two points. For comparison, prior sampled models are shown in Fig.S13. Gray-shaded contours in panel A indicate the 50, 75, and 90%-credible intervals computed from the distribution of models inverted using a purely seismic parameterization[START_REF]Materials and methods are available as supplementary materials[END_REF].

( C )

 C Fig. S14.

Fig. 3 .

 3 Fig. 3. Seismic amplitude behavior with distance. Blue and red bands show the predicted S/P

Figure 3 SM 1

 31 Figure 1

Fig. S1 .

 S1 Fig.S1. Ray paths for the seismic body wave phases (direct P-and S-wave and surface reflections PP, PPP, SS, and SSS) considered in this study. We assume a spherically symmetric medium and use a source located at 50 km depth.

Fig. S2 .

 S2 Fig. S2. Filter bank for event S0235b centered around P-wave (vertical) and S-wave (north/south and east/west) arrivals. Filter bands are half an octave wide on each side, centered on the frequencies marked on the vertical axis. Window length used for smoothing is 8 s. Underlying time domain data is in acceleration after instrument correction. Bold coloured filter banks are the frequencies and components on which the picks were based. The blue dashed lines mark the picks used in this paper. Vertical arrow on horizontal component indicates a possible depth phase. Vertical and horizontal gray bars indicate the presence of glitches and the time period used for normalizing envelopes, respectively. Black traces at the bottom show the unfiltered seismograms from which the filter banks are constructed.

Fig. S4 .

 S4 Fig. S4. Alignment of polarization filtered low-frequency (LF) and broadband (BB) events. (a) P-wave alignment on the vertical component (BHZ). (b) S-wave alignment on the transverse component (BHT), except for S0167b event, which is shown on the radial component (BHR). The blue lines indicate the envelopes of time-series data.The purple lines denote the PP, SS, PPP, and SSS travel time picks, and the purple shaded regions are the picking uncertainties. The red lines represent the predicted travel times from vespagram analysis. Note that no vespagram predictions for S0185a and S0167b are included because the distances of these events are too large compared to other events. Letter "G" represents the instrumental glitches which are also highlighted as faded gray sections of the waveforms(54).

Fig. S5 .

 S5 Fig. S5. Vespagrams of LF and BB event envelopes on (a) vertical (BHZ) and (b) transverse (BHT) components. The predicted slowness and arrival times of body waves are shown in red crosses. These predictions are calculated using the reference model in this study at the reference distance (31.59 • , S-P time=185.6 s). The vespagram stacked envelopes of (c) P waves on BHZ component and (d) S waves on BHT component, these plots are from in the red boxes highlighted in (a) and (b) respectively. The gray curves are the stacked envelopes at different slowness. The black curves highlight the envelopes of multiple reflected waves at the corresponding measured slowness. The red lines denote the predicted travel times at the reference distance using the reference model in this study.

Fig. S7 .

 S7 Fig.S7. Results of waveform matching for event S0235b for a 15 s window length. The solid red lines indicate the proposed SS arrival obtained with the filter bank method (section SM 2.2). Magenta and red dashed lines are the SS and SSS arrivals picked with the vespagram method (section SM 2.3). Panels A and B show the normalised least-squares variance on the transverse horizontal component, which is obtained by fitting in the coda the scaled first 15 s of the S-wave after Hilbert transformation (A) and no Hilbert transformation (B). The black and blue lines represent the best fit for the different low pass frequencies (see text) and the average fit, respectively. Panels C and D show the relative differences between the scaling coefficients, when computed for the 6 different bandwidth filters (see Fig.S8). The scaling coefficient can be interpreted as the relative amplitude ratio of SS or SSS phase with respect to S. Panels E and F indicate average values of the scaling coefficient, again for the 6 bandwidth filters and with respect to the S-wave packet. For optimally matching a phase, the scaling coefficient can be interpreted as the relative amplitude ratio of this phase with respect to S. For S0235b, the best match is obtained for 22.5±2.5 s, as indicated by the green arrow (panels A, C, and E). For the SSS phase of S0235b, we propose an arrival at 38±3 s, as indicated by the green arrow in panels B, D, and F. Time after S refers to the MQS-picked S-wave arrival time.

Fig. S8 .

 S8 Fig.S8. Low-pass filtered (from 0.4 Hz to 0.9 Hz) transverse-component waveforms for event S0235b. All traces are in displacement and have been high-pass filtered at 0.3 Hz and normalized to 0.1 (in order to fit with the 0.1 Hz frequency increment on the figure). In panel A, the S-wave (gray) is compared to the Hilbert transform of the S-wave, delayed by 22.5 s and scaled by -0.2 (blue). A delay of 22.5 s provides a good match to the waveforms between 185 s and 210 s, suggesting that all these phases have gone through a surface reflection and a caustic. The window used for the least-squares waveform fitting is given by the two vertical red dashed lines, while the S waveform is given by the two green dashed lines. The green solid line is the MQS S-wave arrival time, and the red solid line is the proposed SS-wave arrival time, using the inverted SS-S delay with respect to S (Fig.S6). The waveform matching remains good for all bandwidths from 0.3 Hz to 0.9 Hz. In panel B, the same waveform matching is attempted without performing the Hilbert transform of the S-wave pulse. Even if some parts can be matched with a positive reflection coefficient, the matching is not achieved at longer periods beyond the 20 s window, illustrating that this phase cannot be explained with a Hilbert transform.

Fig. S9 .

 S9 Fig.S9. Summary of the seismic observations and body wave arrivals as a function of S-P travel time difference (epicentral distance). (A) Vertical-component spectral envelopes and (B)-(C) three-component (vertical, north, and east) time-domain envelopes of acceleration seismograms after correction for instrument response. The alignment in (A) employs P-and S-wave travel time profiles that are computed from characteristic models (Fig.S15) for a 50-km-deep source, and assuming visible main energy packages are direct P-and S-waves. Events are aligned by matching their P-and S-wave envelopes with the theoretical travel times, with the exception of S0183a for which only an impulsive P-wave envelope is visible. Instead, S0183a is aligned based on the S-to-P amplitudes (see main text and Fig.3). This event is shown in gray to indicate that it is not considered in the inversion. (B)-(C) Bold vertical red and blue lines indicate P, PP, PPP, and S, SS, SSS picks, respectively, whereas vertical red and blue shaded areas denote the pick uncertainties. The coloured vertical dashed gray and dotted red and blue lines indicate body wave picks obtained using independent methods: dotted (polarisation filter and vespagram) and dashed (waveform matching). Downward pointing arrows indicate possible depth-related seismic phases (pP, pS, sS and sP) that are not used in the inversion. In all panels, occurrence of glitches or wind-related artefacts are indicated by light gray-shaded portions. Red and purple lines above the filter banks in (B)-(C) show their sum.

8 CFig. S10 .

 8S10 Fig. S10. Inverted radial profiles of S-wave (A) and P-wave (B) velocity for the seismic model parameterization. Gray-shaded areas represent the a priori distributions on P-and S-wave velocities. (C) sampled posterior P-and S-wave velocity distributions in the various layers (layer numbers are indicated on the upper part of each panel).

Fig. S12 .

 S12 Fig. S12. Parameterization of the Martian crust and mantle geotherm. T surf , T , and T lit indicate temperatures at the surface and bottom of the crust and lithosphere, respectively, whereas Z and Z lit refer to crustal and lithospheric thickness. The blue and red lines indicate conductive and adiabatic parts of the geotherm, respectively.

Fig. S13 .

 S13 Fig.S13. Sampled prior geotherms (A) and S-and P-wave velocity profiles (B-C). Insets in panels (A), (B), and (C) show the distribution of sampled prior lithospheric geothermal gradients (dT/dz) and crustal S-and Pwave velocity structure, respectively. Profiles are colour-coded according to lithospheric thickness: 300-400 km (green), 400-500 km (blue), and 500-600 km (red). The lithospheric thermal gradient is determined from the temperature at the crust-mantle interface and the bottom of the lithosphere and the difference in depth between the two points.

Fig. S14 .

 S14 Fig.S14. Differential body wave travel time misfits for all sampled models. Green and purple lines denote differential travel times computed using the inverted models and squares and circles indicate the observations including errorbars. For the travel time calculations performed here, we always pick the first arrival.

Fig. S15 .

 S15 Fig. S15. Seismic velocity profiles characteristic of models with thick (red) and thin (blue) lithospheres. The models are representative of the inverted models shown in Fig. 2A.

0. 1

 1 to 1.0 Hz. S-and P-wave amplitudes are measured on the vertical component of motion by finding the peak envelope amplitude within ±7.5 s of the predicted S-or P-wave arrival time from the TauP Toolkit (73) for each model respectively. Results are smoothed with a 7-point horizontal smoothing operator that smooths over 7 points in distance, which corresponds to 7 • around the the center point (3 • on either side of the center point). The S-P travel times and corresponding epicentral distances for marsquakes are from the time-domain envelope picks (bold numbers in Table

SM 5

 5 Fig.S16. S-and P-wave attenuation models used for the synthetic amplitude computations.

Fig. S18 .

 S18 Fig.S18. Distributions of physical quantities associated with the parameterized thermo-chemical evolution of Mars and corresponding to a subset (286 random samples) of the red family displayed in Fig.S17. This subset corresponds to thermo-chemical histories leading to present-day crustal thicknesses thinner than 55 km in line with the seismic observations (see insets in Fig.2A). The prior range displayed in Fig.S17is shown with the same color. See text for further details.

Table S2 .

 S2 Summary of body wave phase picks using polarisation filters and vespagrams. All travel times are relative to the main P-wave arrivals.

	Event	PP	PPP	S	SS	SSS
		(s)	(s)	(s)	(s)	(s)
	S0235b 21.5±5 39.5±5 171.1±2.5 193.9±5 205.6±5
	S0407a	23.1±5 47.2±5 171.0±2.5 193.3±5 205.5±5
	S0484b 19.9±5 40.2±5 173.4±2.5 194.6±5 208.6±5
	S0173a	19.2±5 37.2±5 172.8±2.5 197.5±5 213.6±5
	S0189a	20.8±5 36.9±5 195.1±2.5 228.3±5 253.2±5
	S0325a	30.1±5 52.8±5 230.2±2.5 260.2±5 282.3±5
	S0185a* 54.0±5 95.0±5 360.4±2.5 403.1±5 434.0±5
	S0167b* 53.0±5 78.6±5 414.5±2.5 465.9±5 493.1±5
	*not used in the vespagram analysis (see text).

Table S5 .

 S5 Seismic model parameterization and prior model parameter information.

	Description	Quantity Parameter Value/Range	Distribution
	Crustal properties 2	V S	2-4.2 km/s	uniform
		-	V P	α•V S	V P and V S cannot
					decrease with depth
	V P /V S scaling	1	α	1.65-1.85	uniform
	Mantle properties	6	V S	3-5.5 km/s	uniform
		6	V P	6-10 km/s	uniform
	Crustal layers	-	Z	0-15 km	fixed
		-		15-50 km
	Depth nodes	-	Z	50, 100, 200,	fixed
		-		400, 600, 1000,
		-		1588 km
	Epicentral distance 8	∆	0 • -180 •	uniform

Table S7 .

 S7 Summary of the main physical parameters used for modeling the thermo-chemical evolution of a Mars-like planet. Brackets in the fourth column indicate explored range of values. From(32).

	Symbol Meaning	Reference value Explored range
	T m 0	Initial Mantle temperature	1700 K	[1700-2000] K
	T c 0	Initial Core temperature	2100 K	T m 0 + [100-500] K
	T s	Surface temperature	220 K	-
	ρ m	Mantle density	3500 kg/m 3	-
	ρ c	Core density	6500 kg/m 3	-
	ρ cr	Crust density	2900 kg/m 3	-
	g	Surface gravity	3.7 m/s 2	-
	g c	Core surface gravity	3.1 m/s 2	-
	C m	Mantle specific heat at constant pressure	1142.0 J/kg•K	-
	C c	Core specific heat at constant pressure	840 J/kg•K	-
	C cr	Crust specific heat at constant pressure	1000 J/kg•K	-
	k m	Mantle thermal conductivity	4.0 W/m•K	-
	k cr	Crust thermal conductivity	2.5 W/m•K	[2-6] W/m•K
	k c	Core thermal conductivity	40 W/m•K	-
	α m	Silicates thermal expansion coefficient	2.0 10 -5 K -1	-
	α c	Core thermal expansion coefficient	1.5 10 -5 K -1	-
	R p	Mars radius	3.3895 10 6 m	-
	R c	Mars core radius	1.800 10 6 m	-
	c	Ratio of the mean and upper core temperature 1.05	-
	u 0	Mantle convective velocity scale	2 10 -12 m/s	-
	Λ	Crustal enrichment factor	10	[5-20]
	C U	Initial U bulk concentration in the mantle	16 ppb (ng/g)	[1-40] ppb
	C T h	Initial Th bulk concentration in the mantle	56 ppb (ng/g)	4×U
	C K	Initial K bulk concentration in the mantle	305 ppm (µg/g) 2 10 4 ×U
	Ra c	Critical value of the Rayleigh number	450	-
	L m	Mantle latent heat of fusion/crystallization	6 10 5 J/kg	-
	R	Gas constant	8.31 J/K•mol	-
	E *	Effective activation energy	300 kJ/mol	[100-500] kJ/mol
	V *	Effective activation volume	5 cm 3 /mol	[0-10] cm 3 /mol
	P ref	Reference pressure	3 10 9 Ps	-
	T ref	Reference temperature	1600 K	-
	η 0	Reference viscosity	10 21 Pa s	[10 20 -10 23 ] Pa s
	D cr	Crustal thickness	65 km	-

SM 1 SEIS instrument . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . SM 2 Seismic body wave picks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . SM 2.1 Seismic ray paths and body wave nomenclature . . . . . . . . . . . . . . . SM 2.2 Body wave picks using filter banks . . . . . . . . . . . . . . . . . . . . . . SM 2.3 Body Wave Picks using polarization filters and vespagrams . . . . . . . . . SM 2.3.1 Signal processing steps . . . . . . . . . . . . . . . . . . . . . . . . . SM 2.3.2 Vespagram analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . SM 2.4 Body wave picks using waveform matching . . . . . . . . . . . . . . . . . SM 2.5 Summary of body wave picks . . . . . . . . . . . . . . . . . . . . . . . . . SM 3 Stochastic inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . SM 3.1 Sampling the posterior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . SM 3.2 Seismic model parameterization . . . . . . . . . . . . . . . . . . . . . . . . SM 3.3 Geophysical model parameterization . . . . . . . . . . . . . . . . . . . . . SM 3.3.1 Prior information . . . . . . . . . . . . . . . . . . . . . . . . . . . . SM 3.3.2 Sampling the posterior . . . . . . . . . . . . . . . . . . . . . . . . . SM 4 Martian seismic velocity models and amplitude computations . . . . . . . . . . . SM 4.1 Computing synthetic waveforms and amplitudes . . . . . . . . . . . . . . . SM 5 Constraints on thermo-chemical history . . . . . . . . . . . . . . . . . . . . . . . SM 5.1 Thermo-chemical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . SM 5.2 Compositional models of the Martian crust . . . . . . . . . . . . . . . . . . SM 5.3 Thermal gradients from lithospheric flexure . . . . . . . . . . . . . . . . . SM 6 Geodynamic model parameters . . . . . . . . . . . . . . . . . . . . . . . . . . .

Acknowledgments

We acknowledge NASA, CNES, partner agencies and Institutions (UKSA, SSO, DLR, JPL, IPGP-CNRS, ETHZ, IC, MPS-MPG) and the operators of JPL, SISMOC, MSDS, IRIS-DMC and PDS for providing SEED SEIS data. We are grateful to four reviewers for constructive comments that led to an improved manuscript. Funding: A.C.D., A.K., D.G., J.C., M.v.D., and S.C.S. acknowledge support from ETHZ through the ETH+ funding scheme (ETH+02 19-1:"Planet Mars"). M.P.P. was supported by the NASA InSight mission and funds from the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. D.K., F.K., N.S., Q.H., and R.R.M. were supported by NASA-Grant 80NSSC18K1628. V.L. was funded by a Packard Foundation Fellowship. B.P., C.M., H.S., M.D., M.P., and P.L. acknowledge the support of CNES for the SEIS analysis, in addition to an ANR grant (MAGIS, ANR-19-CE31-0008-08). W.T.P. and C.C. acknowledge support from the UK Space Agency. Numerical computations were partly performed on the S-CAPAD platform, IPGP, France. S.M.M. was supported by NASA Grant 80NSSC18K1622. This work was further supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID s922. This is InSight contribution number 129. Author contributions: S.C., D.G., A.K., P.L., N.S., J.F.C., S.C.S., M.v.D., and Q.H. analysed the seismic data and made arrival time picks. A.K. and A.C.D. performed the inversions. A.K., S.C., D.G., M.v.D., N.S., H.S., C.M., and S.E.S. participated and contributed to the interpretation of the

Supplementary materials

Fig. S6. Illustration of waveform matching for the two largest events and their normalized displacement seismograms in the 0.3-0.9 Hz bandwidth. The blue traces correspond to the Hilbert transform of the transverse component, while the gray traces are their opposites (sign-flipped) with specific scaling for each event (0.5 for S0173a and 0.25 for S0235b). Delays, here ∆(SS-S), are indicated in the figure legend. Time after S refers to the MQS-picked S-wave arrival.

SM 2.5 Summary of body wave picks

The 8 low-frequency events considered in this study are shown as vertical-component spectral envelopes (Fig. S9A) and as three-component time-domain envelopes (Fig. S9B,C). All body wave picks made with the three different methods are compiled in Table S4. Only the TDE picks that have been confirmed by at least one other method are considered in the inversion and in table S4.

The time-domain envelopes are a compact way of representing seismic time series, whereas the spectral envelopes provide an overview of the distance-travel time behaviour of the seismic observations. For computing the spectral envelopes, we first calculate the acceleration spectrograms from detrended and instrument response-removed seismic waveforms, subsequently sum the spectral amplitudes along the frequency axis between two frequencies. This allows us to obtain total energy as a function of time. The frequency bands employed are specifically tailored for each event (0.2-0.9 Hz for most of the events), determined from visual inspection of the spectra of the events where these are most clear. The prior model parameter information described above is summarized in table S6 below. The chosen prior ranges represent the information acquired from data and results from experimental and numerical studies as discussed in the foregoing sections. Sampled prior geotherms and Sand P-wave velocity profiles are shown in Fig. S13.

Sampling the posterior

To ensure convergence of the McMC algorithm, we monitored the time series of all model parameters from the algorithm to verify that these remained stationary during the iterations. We also checked for similarity of model parameter values across independent chains. To further assess the performance of the McMC algorithm, we monitored the changes to the probability densities and deemed an adequate number of samples as having been sampled, once there were

The initial HPE content was varied by randomly sampling C U between the bounds indicated in Table S7, and the concentrations of Th and K were set using the following ratios:

and C K /C U = 2 • 10 4 , respectively. The resulting variations in initial bulk heat production rates ranged from 1.6 pW/kg to 10 pW/kg at present. For each set of the aforementioned parameters, we computed the thermo-chemical evolution of Mars and compared the resultant present-day areotherm with the characteristics of the geotherms corresponding to the two different families displayed in Fig. 2B: the blue family characterized by hotter and thinner lithospheric geotherms and the red family that corresponds to colder and thicker lithospheric geotherms.

Upon sampling the model space, we retained the thermo-chemical history model solutions that led to the desired lithospheric geotherm characteristics displayed in Fig. 2B. We rely on the following criteria for present-day quantities:

1. A lithospheric thickness (Z lit ) ranging between 400 km and 500 km for the blue family, and between 500 km and 600 km for the red family.

2. A lithospheric temperature gradient ranging between 1.75 K/km and 2.4 K/km for the blue family, and between 1.6 K/km and 2.4 K/km for the red family.

3. Uppermost mantle temperatures ranging between 1680 K and 1790 K for the blue family and between 1680 K and 1785 K for the red family.

4. Following application of the above criteria, we consider the remaining models with crustal thicknesses <55 km.

The corresponding distributions of physical quantities that satisfy the above criteria are displayed in Fig. S17. Panels a-h and i-p are dedicated to the red and blue family of models shown in Fig. 2B, respectively, along with the prior range (displayed in transparent light blue). The red family is associated with small to moderate reference viscosities (η 0 = 10 20.5 -10 22.3 Pa s, early planetary cooling, and vice versa).

As shown in Fig. 2B (see histograms in inset), the red family is statistically more likely than the blue family. While the geodynamic quantities associated with red geotherms are generally not strongly constrained (Fig. S17a-h), further insight into the value of the geodynamic quantities can be gained by considering constraints on the present-day crustal thicknesses, as suggested in Fig. 2A, where the crust-mantle interface depth appears to be shallower than 55 km.

The corresponding subset of thermo-chemical histories is shown in Fig. S18. It represents a considerably smaller fraction of the possible range associated with the red or blue geothermal families, with values that span a continuous but smaller prior range. The mantle rheology is constrained to E * > 200 kJ/mol (Fig. S18a), V * < 7 cm 3 /mol (Fig. S18b), η 0 = 10 20.2 -10 21.8 Pa s (Fig. S18c). The bulk mantle HPE content (C U = 10 -25 ppb, Fig. S18e) is generally relatively high but still compatible with cosmochemical estimates (28,44). These relatively high values for bulk HPE content would tend to trigger larger amounts of shallow melting, which could produce a relatively thick present-day crust. However, these are compensated by the associated strong crustal enrichment (Λ=13-20, Fig. S18f) that implies a strong HPE depletion of the mantle. Contrary to the cases displayed in Fig. S17, the initial thermal state is constrained to relatively cold values (T m 0 < 1785 K, Fig. S18d). These cold initial temperatures also contribute to smaller crustal production, resulting in a relatively thin crust at present (D cr = 41 -55 km, Fig. S18g)). The relatively high bulk HPE content also implies relatively large values of presentday surface heat flow (14-29 mW/m 2 , Fig. S18h) that are in the higher range or even above pre-mission estimates (27,(31)(32)(33)78). Note that most of the aforementioned studies considered smaller HPE contents, which explains such differences. The values associated with this subset are summarized in Table S9. The computed present-day lithospheric geotherms include a deformable thermal boundary layer about 100 km thick (the exact value varies among the models because it is a function of the evolving mantle convective vigour (32)). This lowermost region SM 5.3 Thermal gradients from lithospheric flexure

In the main text, our constraints on the seismically-determined lithospheric thermal gradients (1.6--2.5 K/km) were found to be consistent with the thermal gradients obtained by (35) (1.3-3.2 K/km). The latter range in thermal gradients was determined by converting their upper (450 km) and lower bounds (330 km) on the elastic thickness using an equating bending moment approach (85). In the conversion, we accounted for crustal heat production (∼8.5 mW/m 2 for a 35-km thick crust, ( 86)) and assumed a wet diabase rheology for the crust and a mantle thermal conductivity of 4 W/m K. The strain rate was varied from 10 -14 s -1 to 10 -24 s -1 and the mantle rheology was set to either wet or dry olivine (87).

SM 6 Geodynamic model parameters

The values of the rheological parameters mentioned in the main text along with the main governing thermo-chemical evolution parameters are listed below in Table S7. Tables S8, and S9 summarise the set of physical quantities that could be constrained from comparison of the thermo-chemical history models with the geophysical models (Section SM 5).

Table S8. Range of quantities associated with the two families of lithospheric geothermal gradients inferred from the differential travel time inversion (Fig. 2B), which are compatible with thermo-chemical evolution of a Mars-like planet (see Fig. S17).

Lithospheric Table S9. Range of quantities associated with the lithospheric geothermal gradients inferred from the differential travel time inversion (red family, see Fig. 2B) with the additional constraint on present-day crustal thickness D cr < 55 km, which are compatible with thermo-chemical evolutions of a Mars-like planet (Fig. S18).

Quantity


Céline Eiden

Véronique Pinzani

Marion Laureau

Nicolas Chapet

Juline Beringer

François Roubille

Mustapha Sebbane

Hélène Peyriere

Clomiphene misuse and risk of severe cardiovascular events

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.