Evidence for a liquid silicate layer atop the Martian core

Seismic recordings made during the InSight mission1 suggested that Mars’s liquid core would need to be approximately 27% lighter than pure liquid iron2,3, implying a considerable complement of light elements. Core compositions based on seismic and bulk geophysical constraints, however, require larger quantities of the volatile elements hydrogen, carbon and sulfur than those that were cosmochemically available in the likely building blocks of Mars4. Here we show that multiply diffracted P waves along a stratified core–mantle boundary region of Mars in combination with first-principles computations of the thermoelastic properties of liquid iron-rich alloys3 require the presence of a fully molten silicate layer overlying a smaller, denser liquid core. Inverting differential body wave travel time data with particular sensitivity to the core–mantle boundary region suggests a decreased core radius of 1,675 ± 30 km associated with an increased density of 6.65 ± 0.1 g cm−3, relative to previous models2,4–8, while the thickness and density of the molten silicate layer are 150 ± 15 km and 4.05 ± 0.05 g cm−3, respectively. The core properties inferred here reconcile bulk geophysical and cosmochemical requirements, consistent with a core containing 85–91 wt% iron–nickel and 9–15 wt% light elements, chiefly sulfur, carbon, oxygen and hydrogen. The chemical characteristics of a molten silicate layer above the core may be revealed by products of Martian magmatism.


S1000a
Since a comprehensive description of the identification of PP-and SS-wave arrivals is given in [1], here we focus mainly on the identification of the SKSwave arrival.Through the application of filter banks, polarized waveforms, and polarization analysis [2,3], we refined the initial picks and considerably reduced the uncertainties of [4], allowing us to properly align the event and better constrain its location.Envelopes of the event with the PP-and SS-wave picks and their uncertainties represented as vertical lines and bars, respectively, are shown in Figure S1.
As the location of the event is well known [5], waveforms are rotated to and envelopes are shown in the Vertical-Radial-Transverse (ZRT) system.Gray curves in Figures S1A-C represent the time-domain envelopes of the bandpass filtered waveforms (0.2-0.5 Hz), while blue curves represent envelopes of the waveforms after applying a time-domain polarization filter [6].Because body waves exhibit a high degree of linear polarization, the application of the polarization filter increases the signal-to-noise ratio of the linearly polarized part of the signal.
Figure S1D shows a three-component scalogram of the event computed as the absolute value of the continuous wavelet transform of the signal.This scalogram, built as the sum of the squared scalograms of each component, illustrates the change in frequency content of the signal with time.While the main PPand SS-wave arrivals are discernible on account of the broad frequency content of the envelopes and scalogram, the identification of additional phases in the PP-and SS-wave coda is challenging, because of possible interference of phases and scattering signal after the PP-wave arrival, as illustrated in Figure S1.In addition, non-seismic signals in the form of large glitches are visible in Figure S1D, (smaller glitches are indicated by the orange box).Glitches, which are transient one-sided pulses often accompanied by high-frequency spikes [7], are mainly attributed to SEIS-internal stress relaxations or by tilts of the instrument [8], and contaminate the signal making the identification of seismic phases difficult.
To narrow the search for the SKS phase, we predict the arrival time of SKS by considering the Martian models of [3] and perturb the P-wave velocity in the core by ±15% to account for the lack of information on seismic core properties.For S1000a, the so-predicted time range is represented by the green shaded area in Figures S1A-C.While the time-domain envelopes of the nonpolarized waveforms do not exhibit any particularly distinctive arrival within the expected time range, a peak is discernible around 340 s in the R component of the polarized waveform envelopes, where the SKS arrival is expected.This large-amplitude arrival, is strongly linearly polarized and consistent with synthetically-computed SKS waveform amplitudes that indicate a relatively large-amplitude SKS phase (Figure S2).No other large-amplitude arrival compatible with the characteristics of SKS is discernible within the expected time window.As a means of validating the selected seismic phase, we filter the velocity traces in narrow frequency bands, allowing us to better visualize the partitioning of energy across different frequencies.We use a bandwidth half an octave wide around each central frequency, ranging from 1/5.7 to 1/1.4 Hz.The selected filter banks cover the frequency bands where the low-frequency energy of the event is most visible and is mainly devoid of glitch-, donkor atmospheric-related artefacts [8].Donks, interpreted as stress release from temperature cycling on the lander, are pulses of energy with a typical duration of a few seconds, and are generally observed at frequencies higher than 10 Hz [7]. Figure S1E shows the envelopes of the polarized and non-polarized R-component waveforms in a time window that includes the SKS-wave arrival.In the filter banks, we observe a strongly polarized arrival at frequencies ranging from 1/2.8 to 1/1.4 Hz that conforms with our established criteria for identifying seismic phases [3] and which we therefore pick as the SKS phase.

S0976a
Analogous to the analysis presented for event S1000a, Figure S3 shows envelopes of polarized (blue) and non-polarized (gray) waveforms in the Vertical-North-East (ZNE) system (Figure S3A-C) and three-component scalogram (Figure S3D) for event S0976a.While an estimate of the backazimuth for event S0976a based on the observed particle motion of the PP-wave arrival is available [9], the uncertainty is considerable (±25 o ).
In contrast to S01000a, S0976a has a much more impulsive PP-wave arrival, allowing us to easily pick it in the band-pass filtered (0.1-0.5 Hz) polarized and non-polarized waveforms.In addition, its arrival is clearly discernible in the scalogram, characterized by a considerable increase of energy.For the SSwave arrival, due to the presence of a small glitch which affects the waveforms, we mainly rely on the scalogram to identify the energy onset.Both, PP-and SS-wave picks are represented as vertical lines in Figure S3.While the increase of energy is clear in the scalogram for both PP-and SS-waves, the picks also coincide with the break of the slope in the Z-component envelopes for PP and the horizontal components for SS.However, due to the presence of several glitches close to the PP-and SS-wave arrivals and the windy conditions that are present throughout the event [9], we assign larger uncertainties to both picks (represented by the vertical bars in Figure S3) in comparison to event S1000a.Phase picks for S0976a and their uncertainties [1] are shown as vertical solid lines and orange bars, respectively.Envelopes are masked where glitches occurred to avoid misinterpretation of phases.Green band represents range for expected arrival of SKS based on inverted models where the core velocities were perturbed ±15%.
Because S0976a has not been reliably located, the prediction of a potential SKS phase relies on the alignment of the event through PP-and SS-wave arrival times.We follow the procedure for S1000a and compute a range of arrival times for SKS to narrow down the possible arrival time interval of SKS.In comparison to S1000a, where a single strong SKS candidate is apparent, for S0976a a series of similar-amplitude and strongly linearly-polarized arrivals in the predicted time window for SKS are present.The lack of a clear candidate for SKS could be due to 1) lack of an accurate location that would enable a rotation of the traces into the ZRT system and 2) the generally noisy (wind) conditions that characterise the event [9].The lack of a precise location requires analysis of both North and East components.Time-domain envelopes across different frequency bands are shown in Figure S4 for both horizontal components.The time window during which SKS is expected is indicated by the green-shaded area in Figure S4, and several linearly-polarized peaks are present that are labeled from I-IV.Of these, we can rule out I because it would correspond to an SKS arrival that is too early relative to what we would expect based on S1000a (this range is represented by dotted vertical lines).This leaves arrivals II, III, and IV.Of these, arrival II (∼285 s) is only seen on the North component but present across the same frequencies (1/2.8-1/1.4Hz) as S1000a, whereas arrival III (∼310 s) appears on both components and at the same frequencies (1/2 and 1/1.4 Hz).Because of this fundamental uncertainty in the SKS arrival, we abstain from picking the latter for event S0976a.
In [10], several independent methods were employed to pick SKS for S0976a based on either waveform or cross correlation methods.The SKS-PP travel time spread for the five methods covers the range from 297.4 s to 311 s with two picks in the lower range (297.4 s and 298.3 s) and three picks toward the upper range (309 s, 310 s, and 311 s).Irving et al. resorted to averaging all the picks as a means of producing a representative SKS pick (303.9 s) and used the standard deviation of the five picks as indicative of pick uncertainty CONTENTS (7.1 s).This final pick, however, is not consistent with all of the individual picks.Clearly some of the picks represent outliers.In contrast, here we rely on a single processing/picking method for consistency, yet we refrain from picking SKS for S0976a because of the inherent difficulty in making a consistent pick.

Prior model parameter information
Model parameters are illustrated in Figure S5, while model parameters and prior model parameter ranges are summarized in Table S1 below.The prior ranges follow our earlier studies [e.g., 2,3,11,12] and cover wide model parameter ranges.Fig. S5 Thermo-chemical and seismic model parameterisation (crust).An example areotherm is shown in red and defined using the parameters surface temperature T surf (fixed), temperature at the base of the lithosphere T lit , and lithospheric thickness Z lit .Prior ranges on T lit and Z lit are indicated by the grey box.Temperatures below the lithosphere are computed from the entropy of the lithology at the base of the lithosphere.Mantle composition is assumed to be homoegenous, indicated by the single-coloured region from the base of the crust to the top of the liquid silicate layer (LSL).All model parameters and prior ranges are defined in Table S1.For illustration, a self-consistently-computed shear-wave velocity profile for the red areotherm is shown in black (for the composition of [13]).The LSL and liquid core are described using separate 3rd-order Burch-Murnaghan (BM) EoSs and parameterised in terms of adiabatic bulk modulus K 0S , its pressure derivative K ′ 0S , and molar density ρ 0 , referenced at the conditions of the solid mantle-LSL interface for the LSL and LSL-liquid core interface for the core, respectively.The thickness of the LSL is determined by ∆Z, whereas core radius is defined by Rcore.
3 Inversion results for a different mantle composition  We use the spectral element code AxiSEM [14] to model high-frequency waveforms that match the frequency content of the data (up to 1 Hz) assuming an inclined (45 • ) source.Figures S7-S8 show synthetic waveform sections for a melt-layer model obtained here (Figure 2A) and a Martian model from [11].[14] assuming an inclined (45 • ) source.Main body-wave arrivals (P and PP) and arrivals of phases that interact with the liquid silicate layer (LSL) are indicated by dotted and solid curves across the section.Inset at the top provides a zoom into the predicted diffracted Pwaves from the mantle-LSL interface (P LSL diff ) and the LSL-liquid core interface (P CMB diff ) for epicentral distances similar to that of S1000a (∼ 126 o ).Inset at the bottom provides a zoom into the core-reflected P-waves from the top of the LSL (PdP) and the liquid core (PDcDP) for epicentral distances similar to that of S1094b (∼ 59 o ).Black dashed lines demarcate both areas for which zoom in insets are produced.[11].Waveforms are computed employing the spectral element code AxiSEM [14] assuming an inclined (45 • ) source.Main body-wave arrivals (P and PP) and arrivals of phases that interact with the core-mantle boundary (CMB) are indicated by solid curves across the section.Inset at the top provides a zoom into the predicted diffracted P-waves (Pdiff) from the CMB interface for epicentral distances similar to that of S1000a (∼ 126 o ).Inset at the bottom provides a zoom into the core-reflected P-waves from the top of the core (PcP) for epicentral distances similar to that of S1094b (∼ 59 o ).Black dashed lines demarcate both areas for which zoom in insets are produced.
6 Observational evidence for a liquid silicate layer Generally, marsquake seismograms are characterized by scattering coda that complicates identification of seismic body wave phases needed for structural interpretation [2,3,15,16].As a consequence, seismic arrivals, with the exception of the direct P-and S-waves, are not immediately identifiable in the time series.Instead, a number of processing steps have been developed specifically for the purpose of detecting seismic phase arrivals, revolving around the use of narrow-band-filtered time-domain polarized and unpolarized envelopes and waveforms labelled filter banks [2,3].This increases the signal-to-ratio of the linearly-polarized part of the signal that involves body-wave phases [6].Filtered polarized vertical-component traces and their time-domain envelopes for events S1000a and S1094b are shown in Figure S9.The polarized traces, in comparison to the original waveforms (not shown) exhibit packages of energy that are associated with body wave arrivals and allows for improved understanding of frequency-dependent energy partitioning.

Filter banks
We filtered the waveform traces (velocity) in narrow frequency bands using a bandwidth half an octave wide around each central frequency, ranging from 1/5 to 1 Hz, and computed their envelopes.The selected filter banks cover the frequency bands where the low-frequency energy of the event is most visible and is mostly devoid of glitch-, donk-or atmospheric-related artefacts [17].For event S1000a, the main diffracted P-wave arrival (P CMB diff ) is discernible across four filter banks from 1/1.4 Hz to 1/4 Hz, whereas the lower-amplitude Pwave that is diffracted along the mantle-liquid silicate layer interface (P LSL diff ) and the reverberation (P diff ˆLSL P diff ) are visible across the frequency bands from 1/1.4 Hz to 1/2.8 Hz (Figure S9C).The onsets of the individual energy packages in the filter banks coincide with the main arrivals that are visible in the polarized waveform and envelope shown in Figure S9B.Based on the identification procedure of body-wave arrivals described in [3], we identify discrete energy onsets as seismic phases from the filter banks and the polarized and non-polarized traces/envelopes.To qualify as seismic phases, arrivals must be present across different frequency bands and in both, non-polarized and polarized traces/envelopes.None of the other arrivals exhibit these characteristics.Following this, the identified arrivals are selected as possible P LSL diff and P diff ˆLSL P diff picks as indicated in Figure S9.Note that P CMB diff has been identified by [1].The vertical red bar in Figure S9B,C indicate the differential travel time ranges relative to the main PP-wave arrival (for S1000a) computed using the inverted models shown in the main manuscript in Figure 2A.The identified P LSL diff arrival and the predictions based on the inverted models are seen to be in good agreement.

Polarisation analysis
To provide an independent verification of the selected arrivals, we follow [2,3] and consider the particle motion (polarization).A scalogram of the event, computed as the absolute value of the continuous wavelet transform of the signal, is shown in Figure S10A and illustrates how the frequency content of the signal changes with time.Energy is seen to increase over the background noise in the frequency range and starting in the time window for which the inverted models predict the arrival of P LSL diff (delimited by dashed lines).As rectilinear polarization with consistent particle motion is expected for the diffracted P-waves, we compute the polarization attributes of the three-component seismic data and infer the instantaneous back-azimuth following [18].To represent changes in azimuth with the arrival of seismic phases, we compute normalized azimuthal probability density curves in the frequency band 0.2-0.7 Hz considering 10-s long time windows.These are shown in Figure S10B.Until the arrival of the P LSL diff wave, no consistent azimuth is observed.With the arrival of P LSL diff , however, particle motion appears to be consistent with the azimuth based on the imaged location (indicated by the horizontal red line).Following the arrival of P LSL diff , the azimuth generally exhibits a consistent pattern around 30 o across the Pdiff wavetrain.

Waveform matching
Since P LSL diff , P CMB diff , and P diff ˆLSL P diff , follow approximately similar paths, their waveforms are expected to be similar.This provides another means of verifying our selected arrivals.We filter the vertical-and radial-component waveforms (0.2--0.7 Hz) and consider a 10-s long window that includes the observed P CMB diff arrival to build a template trace.The vertical-component event and template trace are shown in Figure S10C in dark gray and red, respectively.Following this, we cross-correlate both traces considering lags every 0.05 s, and compute at each lag a correlation coefficient CC Z and CC R for the vertical and radial component, respectively, defined by the mean of the cross-correlation function.The similarity-coefficient employed to measure resemblance between the two diffracted P-wave arrivals is computed as the root-sum-square of the cross-correlation coefficients for both components ( CC Z 2 + CC R 2 ). Figure S10D shows the temporal evolution of the similarity coefficient (gray) and the threshold (red horizontal line) employed to detect waveforms that match the template trace.The analysis indicates three prominent peaks in the similarity coefficient that are above our defined threshold that are indicated by the vertical coloured lines in Figure S10C and are seen to coincide with the arrival of P LSL diff and P diff ˆLSL P diff , in conformity with the synthetic waveform predictions (Figure S2).Based on the similarity in moveout, waveform, and differential travel time between P CMB diff and P LSL diff , on the one hand, and P diff ˆLSL P diff and P CMB diff , on the other hand, the third phase represents a reverberating P CMB diff in the liquid silicate layer.

Core-interacting P-waves
For event S1094b (Figure S9C), we identify packages of energy in the envelopes and waveforms across the main frequency bands (1/1.4Hz and 1/2 Hz), where seismic arrivals from low-frequency events are expected.Note that for this particular event the main P-wave arrival can only be identified above 1 Hz.The vertical orange/blue bars indicate the differential travel time ranges relative to the main P-wave arrival computed using the inverted models shown in the main manuscript in Figure 2A.However, because the impact event S1094b is overall noisy and given the location of PdP and PDcDP in the PP-wave coda, identifying these arrivals is very challenging even with the filter banks.For these reasons polarization and waveform matching analysis as performed for S1000a is equally difficult.
Finally, Figure S11 shows a waveform section with the travel time curves for PdP and PDcDP superposed.This figure indicates that for S-P differential travel times below 200 s (∆ <35 • ), PdP and PDcDP will be located in the Swave coda and thus largely unidentifiable.According to Figure S11, the optimal epicentral distance range for observation of PdP and PDcDP is ∼35 • -60 • , i.e., S-P differential travel times around 225-350 s.

C
Fig. S9 Body-wave ray path geometry of seismic phases that interact with the liquid silicate layer (LSL) (A) and analysis of vertical-component waveforms of events S1094b (B-C) and S1000a (D-E).Panel (A) shows the predicted ray paths of seismic phases that provide additional constraints on the LSL: diffracted P-waves from the mantle-LSL interface (P LSL diff ) and the LSL-liquid core interface (P CMB diff ), including a diffracted P wave reverberating in the LSL (P diff ˆLSL P diff ); and core-reflected P-waves from the top of the LSL (PdP) and the liquid core (PDcDP).The source locations of the high-quality events S1094b (∆ ∼ 59 o ) and S1000a (∆ ∼ 126 o ) employed to predict the ray paths are represented by black stars.Panel (B) shows the vertical-component waveform (dark gray) and envelope (light gray) for event S1000a (filtered between 0.2-0.7 Hz).The arrivals of P LSL diff , P CMB diff , and P diff ˆLSL P diff are indicated.Also shown is the time window where the inverted models (Figure 2A 2A.Each trace (envelope) is normalised to the S-wave.CF designates the differential travel time distance range for the low frequency events that cluster around Cerberus Fossae, while the remaining events appear to be unrelated to any surface tectonic expression [18].

Analysis of Sdiff
As in the case of the diffracted P-wave (SM section 6.1), the inverted models presented here should give rise to diffracted S-waves along the mantle-LSL and LSL-liquid core interfaces.Following the analysis of [1], we first compute synthetic seismograms to analyze the characteristics that can be expected for Sdiff waves.As before, we employ the spectral element code AxiSEM [14] to model high-frequency waveforms that match the frequency content of the data (up to 1 Hz) assuming, as before, an inclined (45 • ) force.From the computed seismogram, we calculate the time-domain envelope by taking the square root of the summed squared amplitudes of the trace and its Hilbert-transform in 5 s long time windows to better depict the arrivals of energy packages.The time-domain envelopes for the vertical, radial, and transverse components are shown in gray in Figure S12A.While a robust SKS-wave arrival is observed in the radial component, Sdiff, which is expected to arrive with transverse polarization, is barely visible.In line herewith, Sdiff is unindentifiable in the observed seismograms (Figure S5C).8 Geo-and cosmochemical constraints Experimental and cosmochemical arguments can place further constraints on the composition of the Martian core.The solubility of C is 4-4.5 wt% in S-rich (10 wt%) Fe alloy near the conditions of Martian core formation, ∼13 GPa and 2073-2273 K [19].As such, model compositions that exceed 4.5 wt% C are not realistic candidates for the Martian core.The light element content may be further restricted by considering the compositional range of potential Martian precursors.Of the extant set of meteorites, the chondritic meteorites, by virtue of their high volatile contents relative to their achondritic counterparts [20], are the most likely sources of S, C, O, and H. Mass-independent stable isotope anomalies of refractory lithophile elements and O in Mars indicate a strong affinity for ordinary chondrites with a lesser contribution from carbonaceous chondrites [21][22][23].Given that S is less volatile than C during cosmochemical processes [e.g.24], the S/C ratios in any potential chondritic building blocks should be minimum estimates for that of Mars.Both S and C are moderatelyto highly siderophile during Martian core formation [19] such that >99.5% of the total Martian budget of S is expected to reside in the core, increasing to >99.9% for C. Thus, the S/C ratio of the Martian core is comparable to that of bulk Mars, and hence can be used to compare with those of chondritic meteorites.The data of [20] indicate that the S/C ratio ranges from ∼2 (CI) to ∼20 (H, L, LL).Using the abundances of moderately volatile lithophile elements to fix the bulk S content of Mars [12], the maximum S content of the Martian core, given its smaller mass fraction (21±0.5%), is 12 wt%.Chondritic constraints then place C at between 0.6 wt% and 6 wt%.The upper end of this range is unfeasible because i) it exceeds C saturation and ii) CI chondrites alone are likely not viable precursors for Mars.Based on previously-constrained OC:CC ratios of 9:1 for the moderately volatile element Zn [25], which has a similar volatility to S, we consider the S/C ratio of bulk Mars to be ∼6, but with considerable uncertainty.This would translate into no more than 2 wt% C in the Martian core.However, C-rich compositions around 3-4 wt% C cannot be excluded cosmochemically.Provided that an estimate for S near 9-12 wt% holds for the Martian core, then, together with the temperature of equilibration and the FeO content of the mantle, the O content of the Martian core is readily calculated.The addition of S to Fe-rich alloy causes a strong decrease in the activity coefficient of O, stabilising O in the metal as S increases [26,27].At 2100 K and a mantle FeO content of 13.7 wt% from [12], the expression of [27] returns 1-2 wt% O in the Martian core.

Estimating the Martian mantle liquidus
In order to estimate the liquidus temperature of an hypothetical liquid layer at the base of the Martian mantle, the crysoscopic equation [28] is employed: The form shown above comes from integration of the van't Hoff equation from a reference temperature T 0 and composition x 0 (here, x = Mg/[Mg+Fe 2+ ]), to the temperature, T and composition, x of interest (R is the gas constant).For the simple binary system approximating that of the Martian mantle, which, like that of the Earth, is dominated by olivine [e.g., 12,29], ∆H f us can be estimated by that of forsterite/fayalite; that is, the liquidus temperature relates to that of the reaction: olivine(s) = olivine(l).The value of ∆H f us has been determined at 1 bar for fayalite at 1490 K [30, 89.3±1.1 kJ/mol] and that estimated for forsterite at 2163 K is 114±20 kJ/mol [31].Based on these data, three simplifications are made to estimate the ∆H f us at 20 GPa for the complex Martian mantle mineralogies.The first is that ∆H f us is taken as 110 kJ/mol, which is the value that reproduces the binary forsterite-fayalite liquidus temperature at 1 bar using equation 1 (Figure S14).This value is intermediate between experimentally determined ∆H f us for both endmembers.The second is that we assume ∆H f us does not change significantly as a function of pressure.The suitability of this assumption may be assessed by examining the K Fe−Mg D between olivine and melt as a function of pressure.Experimental data [32,33] indicate only a marginal effect of pressure on the value of K Fe−Mg D , which is counteracted by that of temperature.For our  purposes, the assumption of constant ∆H f us , both as a function of composition and of pressure, is sufficient.The third assumption is that the heat of fusion of olivine is representative of that at the liquidus at the temperature and pressure of interest, which is substantiated by experiments that indicate olivine is a liquidus phase up to ∼15 GPa [34,35].
In order to estimate T 0 , the liquidus temperature determined by [34] for a nominally anhydrous peridotite composition in the simplified CaO-MgO-Al 2 O 3 -SiO 2 (CMAS) system is used.These authors report a liquidus temperature around 2598 K at 20 GPa, for x 0 = 1 (i.e., iron-free).A T 0 of 2600 K is therefore taken.The resultant liquidus depression as a function of Mg# is shown in Figure S15.At a nominal pressure of 20 GPa, it shows a ∼300 K drop in T from Mg# = 1 to Mg# = 0.5 (2288 K), with a further 300 K drop to Mg# = 0.2.The cryoscopic estimate is consistent with the liquidus temperatures of the Mars-like compositions investigated by [35], who report a liquidus temperature around 2473±100 K at 23 GPa.
Fig. S1 Time-domain envelopes of band-pass filtered (0.2-0.5 Hz) polarized (blue) and non polarized (gray) waveforms of (A-C), three-component scalogram (D), and radial-component polarised (colour) and non-polarised time-domain waveforms and envelopes across different frequency bands (E) for event S1000a.Phase picks for S1000a and their uncertainties [1] are shown as vertical dashed lines and orange bars, respectively.Envelopes are masked where glitches occurred to avoid misinterpretation of phases.The vertical green band represents the expected arrival time range of SKS.Orange box in (D) indicates the presence of a small glitch.Envelopes are normalized for visual purposes.Time is relative to the arrival of the PP-wave.

Fig. S2
Fig. S2 Time-domain envelopes (A-C) and three-component scalograms (D) for a synthetic event with a differential travel time T SS -T PP matching the observation for S1000a.To represent an impact, we compute synthetic waveforms simulating an inclined force at the surface for a range of angles, including vertical force.Their time-domain envelopes are shown in gray, while a mean representative envelope is shown in black.Arrival times for PP and SS computed using ray theory are shown as dotted vertical lines, while the arrival time for SKS is shown by a solid vertical line.

Fig. S3
Fig. S3 Time domain envelopes of band-pass filtered (0.2-0.5 Hz) polarized (blue) and non polarized (gray) waveforms of (A-C) and three-component scalogram (D) of event S0976a.Phase picks for S0976a and their uncertainties[1] are shown as vertical solid lines and orange bars, respectively.Envelopes are masked where glitches occurred to avoid misinterpretation of phases.Green band represents range for expected arrival of SKS based on inverted models where the core velocities were perturbed ±15%.
Fig. S4 Time-domain envelopes across different frequency bands (filter banks) of polarized (color) and non polarized (gray) waveforms in the North (A) and East (B) components.Green-shaded area represents predicted range for SKS arrival based on inverted models and alignment of the differential travel time T SS − T PP .Vertical dotted lines indicate predicted range for SKS arrival based on inverted models [3] that fit the SKS travel time of S1000a.Boxes indicate potential SKS-wave arrivals for each component.Amplitudes of traces and envelopes are normalized by their maxima and scaled for better visualisation.Time is relative to arrival of PP seismic energy in the time-domain waveforms.

Fig. S6
Fig. S6 Summary of Mars's interior structure.(A) Inverted S-and P-wave velocity and density profiles.For comparison, black solid and dashed lines in (A) represent the range of core profiles determined previously from seismic core-transiting (SKS) data [10].(B) Body wave ray path geometry for all events (labeled with stars) considered in this study.Colorbar denotes ray path density, i.e., the number of rays passing through a given area, based on the inverted models shown in panel (A), which explains the diffuseness of the ray paths and source locations.The horizontal column underneath "InSight" is the radial sensitivity and computed as the integrated ray path density with epicentral distance.(C) Inverted liquid silicate layer (LSL) properties in the form of mean density (ρ LSL ), mean P-wave velocity ( VLSL P ), and thickness (∆Z).(D) Inverted core properties in the form of mean density (ρcore) and core radius (Rcore).

Fig. S7
Fig.S7Synthetic section of vertical component waveforms for epicentral distances ranging from 50 to 130 o .Waveforms are computed employing the spectral element code AxiSEM[14] assuming an inclined (45 • ) source.Main body-wave arrivals (P and PP) and arrivals of phases that interact with the liquid silicate layer (LSL) are indicated by dotted and solid curves across the section.Inset at the top provides a zoom into the predicted diffracted Pwaves from the mantle-LSL interface (P LSL diff ) and the LSL-liquid core interface (P CMB diff ) for epicentral distances similar to that of S1000a (∼ 126 o ).Inset at the bottom provides a zoom into the core-reflected P-waves from the top of the LSL (PdP) and the liquid core (PDcDP) for epicentral distances similar to that of S1094b (∼ 59 o ).Black dashed lines demarcate both areas for which zoom in insets are produced.

Fig
Fig.S8Synthetic section of vertical component waveforms for epicentral distances ranging from 50 to 130 o employing an earlier Martian model[11].Waveforms are computed employing the spectral element code AxiSEM[14] assuming an inclined (45 • ) source.Main body-wave arrivals (P and PP) and arrivals of phases that interact with the core-mantle boundary (CMB) are indicated by solid curves across the section.Inset at the top provides a zoom into the predicted diffracted P-waves (Pdiff) from the CMB interface for epicentral distances similar to that of S1000a (∼ 126 o ).Inset at the bottom provides a zoom into the core-reflected P-waves from the top of the core (PcP) for epicentral distances similar to that of S1094b (∼ 59 o ).Black dashed lines demarcate both areas for which zoom in insets are produced.
Fig.S9Body-wave ray path geometry of seismic phases that interact with the liquid silicate layer (LSL) (A) and analysis of vertical-component waveforms of events S1094b (B-C) and S1000a (D-E).Panel (A) shows the predicted ray paths of seismic phases that provide additional constraints on the LSL: diffracted P-waves from the mantle-LSL interface (P LSL diff ) and the LSL-liquid core interface (P CMB diff ), including a diffracted P wave reverberating in the LSL (P diff ˆLSL P diff ); and core-reflected P-waves from the top of the LSL (PdP) and the liquid core (PDcDP).The source locations of the high-quality events S1094b (∆ ∼ 59 o ) and S1000a (∆ ∼ 126 o ) employed to predict the ray paths are represented by black stars.Panel (B) shows the vertical-component waveform (dark gray) and envelope (light gray) for event S1000a (filtered between 0.2-0.7 Hz).The arrivals of P LSL diff , P CMB diff , and P diff ˆLSL P diff are indicated.Also shown is the time window where the inverted models (Figure 2A) predict the arrival of P LSL diff (red vertical bar).(C) Waveform traces (dark color) and their time-domain envelopes (light color) filtered across narrow frequency bands (filter banks).(D) polarized vertical-component waveform (dark gray) and envelope (light gray) for event S1094b (filtered between 0.2-0.7 Hz) in the time window where the inverted models predict the arrivals of PdP and PDcDP (orange and blue vertical bars).Accompanying filter banks are illustrated in panel (E).Each trace is filtered using a bandwidth half an octave wide around each indicated central frequency.Amplitudes of traces and envelopes in (C) and (E) are normalized by their maxima for better visualization.
Fig. S10 Analysis of event S1000a in a time window showing the arrivals of the diffracted P-waves.The three-component scalogram of event S1000a and the temporal change in azimuthal density in the 0.2-0.7 Hz frequency band are shown in panels (A) and (B), respectively.For comparison, the horizontal red line in (B) represents the azimuth based on the imaged location.Panel (C) shows the vertical-component waveform (dark gray) and envelope (light gray) for event S1000a (filtered between 0.2-0.7 Hz).The vertical-component template trace employed for waveform matching is superimposed in red.Solid coloured lines in(A-C) indicate the times at which the similarity between the template and event trace is higher than 1.The similarity between the event trace and template trace is computed as the root sum squared of the cross-correlation coefficients obtained for the vertical (CC Z ) and radial (CC R ) components, and is shown in panel (D).Horizontal red line in (D) indicates the threshold employed for the waveform matching detections.Horizontal red bar on top and vertical red dotted lines in panels A-C represent the predicted arrivals of P LSL diff based on the inverted models.
Fig.S11Waveform section showing bandpass filtered (0.2-0.8 Hz) vertical-component envelopes (color) and predicted travel time curves for the direct P-and S-waves and interfaceand core-reflected P-waves (PdP and PDcDP) based on the inverted seismic models shown in Figure2A.Each trace (envelope) is normalised to the S-wave.CF designates the differential travel time distance range for the low frequency events that cluster around Cerberus Fossae, while the remaining events appear to be unrelated to any surface tectonic expression[18].

Fig. S12
Fig. S12 Time-domain envelopes (A) and three-component scalograms (B) for a synthetic event with a differential travel time of T SS − T PP matching the observation of S1000a.Synthetic waveforms are computed employing an inclined (45 • ) force as source.Arrival times for PP, SKS, Sdiff, and SS computed using ray theory are shown as dotted vertical lines.Time is relative to the theoretical PP-wave arrival.
Fig.S13Mars's core composition and light element budget.Each gray-shaded circle represents a Martian core composition in the liquid Fe-Ni-X (X=S,C,O,H) senary system, which has been obtained by matching density of the Fe-Ni-X mixtures based on the ab initio molecular dynamics (AIMD) simulations with the inverted seismic core profiles (corresponding to the light blue profiles in Figure4in the main text).The inset shows the subset of core compositions that remain after application of geo-and cosmochemical constraints (see Methods and SM section 8).
Fig.S14The forsterite-fayalite binary loop at 1 bar pressure.The thermodynamic properties determined are shown in black, and the approximation of the liquidus curve using eq. 1 for T 0 = 2163 K and ∆H f us = 110 kJ/mol shown at 0.1 Mg# intervals are shown in red.
Fig.S15The estimated liquidus of the Martian mantle with a peridotitic composition at 20 GPa as a function of its Mg# using the cryoscopic equation (eq. 1) at constant T 0 = 2600 K and ∆H f us = 110 kJ/mol.

Table S1
Overview of model parameters, their quantity, prior model ranges and probability distributions employed in the geophysical parameterization