# Seismic Velocity Variations in a 3D Martian Mantle: Implications for the InSight Measurements

A.-c. Plesa, E. Bozdağ, A. Rivoldini, M. Knapmeyer, S. Mclennan, S.

Padovan, N. Tosi, D. Breuer, D. Peter, S. Stähler, et al.

► To cite this version:

A.-c. Plesa, E. Bozdağ, A. Rivoldini, M. Knapmeyer, S. Mclennan, et al. Seismic Velocity Variations in a 3D Martian Mantle: Implications for the InSight Measurements. Journal of Geophysical Research. Planets, 2021, 126 (6), 10.1029/2020JE006755. [https://doi.org/10.1029/2020JE006755](https://doi.org/10.1029/2020JE006755)

HAL Id: hal-03365699

https://hal.science/hal-03365699

Submitted on 9 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Planets

## RESEARCH ARTICLE

10.1029/2020JE006755

**Key Points:**

* Models show up to 10% difference between maximum and minimum S wave velocity either below the crust or at the depth of phase transitions

* The seismic velocity pattern in the lithosphere correlates with the crustal thickness dichotomy and can extend to depths >400 km

* Models with a crust containing <20% of the total heat production show shadow zones that are incompatible with current seismic observations

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

A.-C. Plesa,
ana.plesa@dlr.de

**Citation:**

Plesa, A.-C., Bozdağ, E., Rivoldini, A., Knapmeyer, M., McLennan, S. M., Padovan, S., et al. (2021). Seismic velocity variations in a 3D Martian mantle: Implications for the InSight measurements. *Journal of Geophysical Research: Planets*, 126, e2020JE006755. https://doi.org/10.1029/2020JE006755

Received 30 OCT 2020
Accepted 18 APR 2021

**Author Contributions:**

**Conceptualization:** A.-C. Plesa, E. Bozdağ, A. Rivoldini, M. Knapmeyer, S. M. McLennan, S. Padovan, N. Tosi, D. Breuer, D. Peter, S. Stähler, M. A. Wieczorek, M. van Driel, A. Khan, T. Spohn

**Funding acquisition:** T. Spohn

**Investigation:** A.-C. Plesa, E. Bozdağ

**Methodology:** A.-C. Plesa, E. Bozdağ, A. Rivoldini, M. Knapmeyer, S. M. McLennan

**Software:** A.-C. Plesa, A. Rivoldini, M. Knapmeyer

**Validation:** A.-C. Plesa, E. Bozdağ

© 2021. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

# Seismic Velocity Variations in a 3D Martian Mantle: Implications for the InSight Measurements

A.-C. Plesa¹ ID, E. Bozdağ² ID, A. Rivoldini³ ID, M. Knapmeyer¹ ID, S. M. McLennan⁴ ID, S. Padovan¹ ID, N. Tosi¹ ID, D. Breuer¹ ID, D. Peter⁵ ID, S. Stähler⁶ ID, M. A. Wieczorek⁷ ID, M. van Driel⁶ ID, A. Khan⁶ ID, and T. Spohn¹,⁸ ID

¹German Aerospace Center (DLR), Berlin, Germany, ²Colorado School of Mines, Golden, CO, USA, ³Royal Observatory of Belgium, Brussels, Belgium, ⁴Department of Geosciences, Stony Brook University, Stony Brook, NY, USA, ⁵Physical Science and Engineering Division, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia, ⁶Institute of Geophysics, ETH Zurich, Zurich, Switzerland, ⁷Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, France, ⁸International Space Science Institute (ISSI), Bern, Switzerland

**Abstract** We use a large data set of 3D thermal evolution models to predict the distribution of present-day seismic velocities in the Martian interior. Our models show a difference between maximum and minimum S wave velocity of up to 10% either below the crust, where thermal variations are largest, or at the depth of the olivine to wadsleyite phase transition, located at around 1,000–1,200 km depth. Models with thick lithospheres on average have weak low-velocity zones that extend deeper than 400 km and seismic velocity variations in the uppermost 400–600 km that closely follow the crustal thickness pattern. For these cases, the crust contains more than half of the total amount of heat-producing elements. Models with limited crustal heat production have thinner lithospheres and shallower but prominent low-velocity zones that are incompatible with Interior exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) observations. Seismic events suggested to originate in Cerberus Fossae indicate the absence of S wave shadow zones in 25°–30° epicentral distance. This result is compatible with previous best fit models that require a large average lithospheric thickness and a crust containing more than half of the bulk amount of heat-producing elements to be compatible with geological and geophysical constraints. Ongoing and future InSight measurements that will determine the existence of a weak low-velocity zone will directly bear on the crustal heat production.

**Plain Language Summary** The crustal thickness variations and the crustal enrichment in heat-producing elements directly affect the thermal state of the lithosphere and in turn the distribution of seismic velocities in the interior of Mars. Thermal evolution models in a 3D geometry with a crust that contains more than half of the total radioactive heat production show large variations of the seismic velocities in the lithosphere. These models are characterized by a weak low-velocity zone that extends locally to depths larger than 400 km and a seismic velocity pattern similar to the crustal thickness pattern down to 600 km depth. Models, with limited crustal heat production, and hence higher mantle heat production, lead to a thinner lithosphere that results in shallower but more prominent low-velocity zones. The latter produce S wave shadow zones that are incompatible with clear S-phase arrivals for events located close to Cerberus Fossae. The absence of S wave shadow zones between the Interior exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) landing site and Cerberus Fossae is in line with other geological and geophysical constraints that require a large fraction of heat-producing elements to be located in the Martian crust. Future InSight measurements will put further constraints on the distribution of heat-producing elements in the Martian interior.

## 1. Introduction

The Interior exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) mission (Bauerdt et al., 2020) landed on November 2018 in the Elysium Planitia region (Golombek et al., 2020) bringing a suite of geophysical instruments to Mars. Equipped with a seismometer (Lognonné et al., 2019), an X-band radio system that together with the NASA Deep Space Network tracking stations determines the rotation of Mars (Folkner et al., 2018), and a heat flow probe (Spohn et al., 2018), as well as with payloads for atmospheric science (Spiga et al., 2018) and magnetic field investigations (Johnson et al., 2020), InSight
**Writing – original draft:** A.-C. Plesa,
E. Bozdağ, A. Rivoldini, M. Knapmeyer,
S. M. McLennan, S. Padovan, N. Tosi,
D. Breuer, D. Peter, S. Stähler,
M. A. Wieczorek, M. van Driel,
A. Khan, T. Spohn

aims at understanding the interior structure and present-day thermal state of Mars and can provide valua-
ble constraints for the planet's thermochemical history (Smrekar et al., 2019).

Since February 2019, the seismometer has continuously recorded Mars' seismic activity (Banerdt et al., 2020; Lognonné et al., 2020), and a total of 465 events have been identified up to Sol 478/March 31, 2020 (Clinton et al., 2021; InSight Marsquake Service, 2020). Among them, two high-quality events on Sols 173 and 235 (S0173a and S0235b, respectively) have distinguishable *P* and *S* wave arrivals and are suggested to originate in Cerberus Fossae (Giardini et al., 2020), a potentially active fault system (J. Taylor et al., 2013). Additional events on Sols 183, 205, and 325 (S0183a, S0205a, and S0325a, respectively) show a dominant *P*-phase and only a weak *S*-phase, but the location of these events has large uncertainties (Giardini et al., 2020). In addition, the presence of glitches at the expected *S*-phase arrivals for two of these events (i.e., S0205a and S0325a) makes their interpretation difficult.

The presence of clear *P* and *S* wave arrivals for the two events localized in Cerberus Fossae (i.e., S0173a, and S0235b) can be used to constrain the lithospheric temperature, as it excludes shadow zones at epicentral distances between 25° and 30°. Shadow zones may be produced by pronounced low-velocity zones, where seismic velocities decrease with depth. This causes a downward refraction of seismic waves and thus prevents waves from arriving at certain distances. The presence of low-velocity zones depends on the lithospheric temperature on Mars (Mocquet & Menvielle, 2000; Zheng et al., 2015).

In this study, we analyze the spatial distribution of seismic velocities for a large set of 3D thermal evolution
models of the interior of Mars (Plesa, Padovan, et al., 2018). Additionally, we perform ray-tracing compu-
tations to exclude models with *S* wave shadow zones at epicentral distances between 25° and 30°, with a
special focus on paths to the S0173a and S0235b events in Cerberus Fossae (Giardini et al., 2020). To this
end, we use the entire 3D temperature field and compute local profiles between the InSight location and
Cerberus Fossae to investigate the *P*- and *S*-raypaths.

Previous studies have investigated seismic velocities in the Martian mantle based on 1D parametrizations (e.g., Bagheri et al., 2019; Franck & Kowalle, 1994; Gudkova & Zharkov, 2004; Khan et al., 2018; Mocquet et al., 1996; Sohl & Spohn, 1997; Zheng et al., 2015). So far, no systematic study of the effects of 3D thermal variations in the interior of Mars on the seismic velocities has been undertaken. Since thermal evolution models suggest that more than 50% of the total heat-producing elements (HPEs) need to be located in the crust to be able to explain available geophysical, petrological, and geological constraints (Plesa, Padovan, et al., 2018; Thiriet et al., 2018), regional crustal thickness variations are expected to directly affect subsurface temperatures and, hence, seismic velocities in the lithosphere. Moreover, the presence and spatial distribution (i.e., local or global) of a LVZ strongly depends on the temperature variations in the interior. While the effects of a global LVZ can be modeled in 1D, a 3D investigation is required to address the consequences of lateral variations in the depth of the LVZ, in particular between the northern and southern hemispheres.

We describe the thermal evolution models, the thermodynamic calculations, and the ray-tracing analysis
used in this study in Section 2. We present the results of the entire 3D thermal evolution data set with
special focus on several end-member and intermediate cases in Section 3. A discussion of the constraints
provided by the absence of shadow zones between InSight and Cerberus Fossae is presented in Section 4,
followed by conclusions in Section 5. We emphasize that the 3D seismic velocity distribution as derived in
this study can be used to estimate the effects of a 3D thermal structure on seismic wave propagation models
(Bissig et al., 2018; Bozdağ et al., 2017) and provide a framework for the interpretation of the InSight seismic
measurements.

**2. Methods**

Here, we employ the entire suite of 3D thermal evolution models previously presented in Plesa, Padovan, et al. (2018) to investigate the effects of the distribution of HPEs and mantle composition on seismic velocity variations in the interior of Mars and to select models that are compatible with the absence of an S-shadow zone between the InSight landing site and Cerberus Fossae.
<table><caption>**Table 1**
Summary of Crustal Thickness and Crustal Amount of HPEs Used in This Study</caption><thead><tr><th>Case id</th><th>Crustal thickness (km)</th><th>Crustal amount of HPEs (% of total bulk HPE)</th><th>Case number (Plesa, Padovan, et al., 2018)</th></tr></thead><tbody><tr><td>lowE<sup>a</sup></td><td>87.1</td><td>19.5</td><td>68</td></tr><tr><td>mediumE<sup>b</sup></td><td>62</td><td>67.8</td><td>85</td></tr><tr><td>highE<sup>b</sup></td><td>87.1</td><td>97.6</td><td>65</td></tr></tbody></table>

Notes. For all cases, the core-mantle boundary is located at 1,550 km depth. For each model, we also list the corresponding case number as presented in Plesa, Padovan, et al. (2018), where details on all additional parameters can be found. Cases lowE and highE represent end-members in terms of their crustal amount of HPEs and temperature profiles. Cases that contain the description mediumE have a moderate amount of crustal HPEs and show an intermediate temperature profile.

Abbreviation: HPE, heat-producing element.

<sup>a</sup>The crustal heat production rate (per kg) is five times lower than the GRS surface estimate. <sup>b</sup>The model assumes that the GRS surface data are representative for the entire crust.

## 2.1. Thermal Evolution Models

The thermal evolution models used in this study investigate a wide range of parameters, including the reference viscosity and pressure dependence of the viscosity, initial mantle temperature, initial core-mantle boundary (CMB) temperature, core radius, constant or temperature- and pressure-dependent thermal expansivity, a thermal-conductivity contrast between the mantle and crust, various crustal thickness models, and distribution of HPEs between the mantle and crust (Plesa, Padovan, et al., 2018).

All models use a crust, whose thickness is constant in time but varies spatially. This is compatible with the fact that the crust was built very early during Mars' history (Nimmo & Tanaka, 2005). The crustal thickness variations were derived from gravity and topography data (Wieczorek et al., 2019) and were used previously in geodynamical models (Plesa et al., 2016).

The final temperature distribution in our models is the result of the 4.5 Gyr thermal evolution of the planet, which takes into account the cooling history of Mars and the generation of heat by radioactive elements. To test the effects of temperature variations on the seismic velocities, we choose end-member and intermediate cases out of a pool of 130

model runs that have been presented in Plesa, Padovan, et al. (2018). Two cases that contain a low (lowE) and high (highE) amount of HPEs in the crust lead to hot and cold end-member temperature profiles, respectively. The CMB is located at 1,550 km depth and is compatible with the range of core sizes as suggested by estimates of the tidal Love number $k_2$ (Khan et al., 2018; Plesa, Padovan, et al., 2018; Rivoldini et al., 2011). Both cases use a thick crust with an average crustal thickness of 87.1 km. While the highE case contains most of the HPEs in the crust (97.6% of the total bulk amount of HPEs), the crust for the lowE case contains less than 20% HPEs of the total bulk amount. Additionally, we use an intermediate case (mediumE) that is one of the best fit cases of Plesa, Padovan, et al. (2018), for which the average crustal thickness is 62 km and 67.8% of the bulk amount of HPEs is located in the crust. For cases highE and medium, the concentration of K, Th, and U in the crust is $K = 3,651$ ppm, $Th = 670$ ppb, and $U = 192$ ppb. These values are a factor of about 12 higher than the bulk concentration of Wänke and Dreibus (1994) for the primitive mantle (i.e., mantle + crust) and give a crustal heat production rate of 49 pW/kg at the present day and 275 pW/kg at 4.5 Ga, which is consistent with the near-surface heat production derived from gamma-ray measurements (GRS, Hahn et al., 2011; G. J. Taylor et al., 2006). Case lowE uses lower crustal HPEs values of $K = 729$ ppm, $Th = 134$ ppb, and $U = 38$ ppb, which are a factor of about 2.4 higher than the bulk concentration of Wänke and Dreibus (1994) and gives a five times lower crustal heat production than the value derived from GRS measurements. The cases highE, mediumE, and lowE correspond to cases 65, 85, and 68, respectively, of Plesa, Padovan, et al. (2018). All models, for which a detailed analysis has been performed here, are summarized in Table 1. Additional models using different crustal thickness and crustal enrichment of HPEs are discussed in the supplementary information (SI, Section S3).

## 2.2. Thermodynamical Models

The seismic velocities for a given Mars mantle composition are calculated at each location in the 3D model from the temperature field using the thermodynamic formulation and thermodynamic database of Stixrude and Lithgow-Bertelloni (2011) with the Perple_X software (Connolly, 2009). The Perple_X readout tables used in this study are included in the SI (Data Sets S1 and S2). To study how mineralogical models affect the distribution of seismic velocities in the mantle, we use the compositional model of G. J. Taylor (2013), which we abbreviate as TAY13, and that of Yoshizaki and McDonough (2020), which will be further referred to as YOS20. The TAY13 model is an update of the Wänke and Dreibus (1994) compositional model, taking into account recent chemical data from surface analyses and Martian meteorites. With no significant differences to the Wänke and Dreibus (1994) compositional model, TAY13 is currently considered as the standard compositional model for Mars. Unlike for Wänke and Dreibus (1994) and TAY13, the YOS20 model does not
assume that the refractory major element composition is the same as that of CI carbonaceous chondrites. In the YOS20 model, the bulk silicate Mars (BSM) is 2.26 times more enriched in refractory lithophile elements than the CI carbonaceous chondrites and is systematically depleted in moderately volatile lithophile elements, including the HPE potassium, although not to the same degree as is the TAY13 model (TAY13, K = 309 ppm; YOS20, K = 360 ppm). The YOS20 model has a higher Mg# for the BSM resulting in an FeO content of 14.7 ± 1.0 wt%. This value is lower than the FeO content of the TAY13 model (18.1 ± 2.2 wt%). We note that the temperature data calculated in the study of Plesa, Padovan, et al. (2018) used the bulk amount of HPEs of the Wänke and Dreibus (1994) compositional model. While the TAY13 has a nearly identical bulk amount of HPEs, the HPE abundances in the YOS20 compositional model result in a 17% larger heat production, which would lead to slightly higher present-day mantle temperatures. Nevertheless, we note that the range of temperatures that result by varying the distribution of HPEs between the mantle and crust covers the range of temperature variations resulting from the two compositional models.

For the crust, we employ a mineralogical model based on S. R. Taylor and McLennan (2009) and calculate the seismic velocities using the mineral physics toolbox BurnMan (Cottaar et al., 2014, 2016). The minerals and their proportions, which have been used to calculate the density and seismic velocities of the crust, have been computed using a CIPW mineralogical norm and the wt% oxides listed in S. R. Taylor and McLennan (2009) excluding MnO, Cr₂O₃, and P₂O₅, since minerals based on these oxides are not present in the mineral databases available for BurnMan. Our calculations use the mineral database of Stixrude and Lithgow-Bertelloni (2011) with the equation of state from Stixrude and Lithgow-Bertelloni (2005). However, we use the ilmenite values from the Holland and Powell (2011) mineral database, since in the mineral database of Stixrude and Lithgow-Bertelloni (2011), only a high-pressure phase was available (akimotoite). We have also assumed that all iron is present as Fe²⁺. Even with these simplifications, the resulting average density of the crust is 3,150 kg m⁻³ in perfect agreement with the crustal density estimate of S. R. Taylor and McLennan (2009). The crustal composition model of S. R. Taylor and McLennan (2009) is based on a combination of Martian regolith and GRS data that was corrected for near-surface volatiles and a 2% meteoritic (CI) component. In detail, K and Th were taken from the GRS global average but increased by a factor of 1.12 to account for H, S, and Cl in near-surface deposits; U was calculated assuming Th/U = 3.8. In this study, we focus on seismic velocity variations in the mantle, and the details of the crustal mineralogy play a secondary role. Therefore, we have kept the crustal mineralogy constant for all models discussed here, and we did not consider the effects of crustal porosity.

## 2.3. Ray-Tracing Models

For the 1D ray tracing, we use the TTBOX package, which was validated and verified by Knapmeyer (2004, 2005). This package automatically takes into account multipathing and returns all geometrically possible arrivals for a given phase at a given distance. It also finds rays, which arrive at the target distance Δ along the long arc 360°−Δ, or after one complete revolution 360° + Δ.

To determine the presence of shadow zones, we compute 1D seismic velocity profiles by averaging our data along a great circle segment between the InSight landing site and Cerberus Fossae. We test several source depths of 50, 100, and 150 km since neither of the events shows a discernible surface wave train.

In a 1D seismic velocity model, a shadow zone exists only if the Flat Earth Transformed velocity (Müller, 1976) has a negative gradient. We analyze the average models and search for shadow zones in four steps: (a) look for negative velocity gradients, (b) compute *P*- and *S*-raypaths for a distance range from 5° to 65°, increasing the takeoff angle at the source in steps of 0.02° to obtain an overview of the focusing/de-focusing of the wavefield, (c) try and shoot *P*- and *S*-rays to the mean epicentral distances of events S0173a (29°) and S0235b (26°), which will not be possible if the model has a shadow zone at these distances, and (d) compute *P* and *S* travel time curves for distances from 0° to 45°, in steps of 0.1°, to determine the extent of possible shadow zones.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Results of model calculations comparing the mantle composition of G. J. Taylor (2013) and Yoshizaki and McDonough (2020) for the cases lowE, mediumE, and highE: (a) average temperature profiles used as input for the two mantle compositions, (b) average density profiles, (c) average S wave velocity, and (d) average P wave velocity profile. The horizontal dashed lines in panel (a) show the average lithosphere thickness of each model. For the crust, we used the composition discussed in S. R. Taylor and McLennan (2009) and in Section 2.2.</figcaption></figure>

**3. Results**

In the following, we present seismic velocities and ray-tracing calculations for the 130 thermal evolution models discussed in Plesa, Padovan, et al. (2018) and the compositional models of G. J. Taylor (2013) and Yoshizaki and McDonough (2020). A detailed analysis is shown for three models described in Table 1.

**3.1. The Effect of Composition on Seismic Velocities**

In Figure 1, we show the effect of the two compositional models on the density and seismic velocities for the hot and cold end-members and an intermediate case. The lower FeO content of the YOS20 model (14.7 ± 1.0 wt%) compared to the value of TAY13 (18.1 ± 2.2 wt%) results in lower densities and slightly lower seismic velocities in the Martian mantle (Figure 1). We note, however, that the differences between the seismic velocities obtained with the two compositional models are relatively minor compared to the variations of these quantities obtained by using the different temperature profiles of this study.

The temperature profiles for cases lowE, mediumE, and highE are shown in Figure 1a. The thickness of the lithosphere, which is shown by horizontal lines in Figure 1a, has been set at the depth where the temperature gradient drops below 1 K/km and marks the transition between the conductive upper part (stagnant lid and thermal boundary layer) and the convective mantle below. The effect of the temperature on the density
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Phase proportions computed below the crust for the mantle composition of G. J. Taylor (2013) (top row) and Yoshizaki and McDonough (2020) (bottom row) for the temperature profiles shown in Figure 1. The mineral assemblage is composed of olivine (o), wadsleyite (wad), ringwoodite (ring), ferropericlase (fp), garnet (gt), orthopyroxene (opx), clinopyroxene (cpx), and C2/c clinopyroxene (C2/c). The left panels (a and d) show the phase proportions for the hot end-member, low amount of crustal HPEs case (lowE). The middle panels (b and e) show the phase proportions for the intermediate thermal profile (mediumE), while the right panels (c and f) show the stable mineral phases obtained for the cold end-member, high amount of crustal HPEs case (highE).</figcaption></figure>

profile is shown in Figure 1b. The densities and seismic velocities rapidly increase with depth in the uppermost ~150 km. While cases lowE and highE have a thicker crust (average crustal thickness of 87.1 km), the average crustal thickness in case mediumE is only 62 km. Hence, the depth at which crustal material is no longer present (not even locally) is shallower compared to the lowE and highE cases. The temperature gradient in the lithosphere can lead to a decrease of density and seismic velocities with increasing depth. This effect is most pronounced for the lowE case, which has the highest plotted temperature gradient. The highE case, in contrast, only shows a minor decrease in the S wave velocity with depth and has a continuously increasing density profile. We note that the regions with where the density decreases with depth for the lowE models in Figure 1b are located entirely within the lithosphere and should thus remain stable over billions of years.

In Figure 2, we illustrate the effect of temperature on the phase proportions. Additional figures for a more detailed comparison are presented in the supplementary information (SI, Section S1). First, the TAY13 compositional model has a higher amount of olivine and its high-pressure phases compared to the YOS20 composition. In contrast to TAY13, the YOS20 composition has a higher amount of clinopyroxene in the upper part of the mantle. In both compositions, the pressure at which the phase transition from olivine to wadsleyite occurs is around 13 GPa (about 1,000 km depth) in the cold temperature case (highE) and increases with increasing temperature for the other two cases (i.e., 1,070 km depth for mediumE and 1,120 km depth
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Results of model calculations comparing three end-member models in terms of crustal heat-producing elements (lowE, mediumE, and highE). In each panel, the full lines show the average profiles of (a) temperature, (b) density, (c) *S* wave velocity, and (d) *P* wave velocity, while the shaded envelopes illustrate their spatial variations at each depth. The average profiles computed between the Interior exploration using Seismic Investigations, Geodesy and Heat Transport landing site and Cerberus Fossae are shown as dashed lines. The average thickness of the lithosphere is plotted as horizontal lines in panel (a).</figcaption></figure>

for lowE, respectively). This phase transition is the most prominent one in our models and is observed both in density and the seismic velocities.

Ringwoodite is stable only for the coldest temperature profile (highE) at depths greater than about 1,300 km, although it may appear locally for the intermediate temperature case (mediumE, see additional discussion in Section S1). Conversely, ferropericlase is more stable at these depths for the warmer thermal profiles (cases lowE and mediumE, respectively) compared to the cold temperature profile, where most heat sources are located in the crust (highE). We observe a negative velocity gradient less than 100 km above the CMB for the cold end-member model (highE) and TAY13 composition, as larger proportions of garnet and ferropericlase become stable at the expense of ringwoodite compared to the YOS20 composition. This effect is visible only locally for the intermediate thermal profile (mediumE) and entirely absent for the hot end-member, for which ringwoodite is not stable in the lower part of the mantle at all (see also Figure 3).

We note that perovskite is not stable for the three end-member cases discussed here, due to the large core that has been used for these models. Perovskite can exist at the base of the Martian mantle, if the pressure and temperature are sufficiently high (Rivoldini et al., 2011; Verhoeven et al., 2005).

The phase transition from orthopyroxene to C2/c clinopyroxene (Finkelstein et al., 2015) takes place between 875 and 1,100 km in the intermediate case (mediumE) for the YOS20 composition, although locally
this phase transition also takes place for the TAY13 model for the same thermal profile. Its effects on the seismic velocities are rather weak as the phase proportion is small. This phase transition does not occur for the hot end-member (low crustal HPEs, lowE) and TAY13 composition, while for the same case and YOS20 composition it is barely noticeable on local scale (see Section S1). However, for the cold end-member, where most of HPEs are located in the crust, this phase transition is clearly visible and takes place at shallower depths (it extends from 600 to 900 km). Its effect on the seismic velocities starts to be visible at around 600 km but is considerably smaller than the effect of the phase transition from olivine to wadsleyite, as the proportion of C2/c is small.

## 3.2. The Effect of Temperature Variations on Seismic Velocities

In Figure 3, we show how the temperature profile (computed after 4.5 Gyr of evolution) affects the computed seismic velocities using the TAY13 mantle composition. A similar comparison for the YOS20 composition is shown in Section S2 of the SI. In addition, to showing 1D profiles averaged over the entire mantle, we also show profiles averaged over the great circle between the InSight landing site and Cerberus Fossae, where the two high-quality events have been located. The average profiles and the profiles between InSight and Cerberus Fossae are similar with the largest difference observed for the cold thermal profile (highE case). However, even for that case, the differences in the seismic velocities are minor, and mostly observed in the shallow part, reflecting the thinner crust between InSight and Cerberus Fossae compared to the average value. Indeed, previous modeling studies have also suggested that the InSight landing site is representative of the global average in terms of surface heat flow (Plesa et al., 2016).

The crustal thickness and the distribution of HPEs between the mantle and crust has a first-order effect on the temperature distribution and lithospheric thickness (Figure 3a). The variations of the lithosphere thickness are shown in Figures S4 and S6 of the SI. The effects of the mantle temperature profile are visible on the density and seismic velocities profiles (Figures 3b-3d). Compared to the highE model, the hot temperature profile of the lowE model results in a pronounced LVZ below the crust and to a 100 km deeper phase transition from olivine to wadsleyite, which lies for all cases in Figure 3 between 1,000 and 1,200 km.

We have computed the 3D seismic velocities distribution based on the present-day temperature variations and according to the prescribed mineralogical models of TAY13 and YOS20. In Figure 4, we show maps of the shear wave velocity at several depths in the mantle for the lowE, mediumE, and highE models with the TAY13 composition. A similar figure for the YOS20 composition is shown in the SI (Section S2). The maps in Figure 4 show depth slices where the seismic velocity pattern changes from one dominated by the crustal thickness pattern to a more homogeneous distribution. The gradient $dV_s/dz$ is shown in the lower half of this plot at the same depths. These maps demonstrate that with increasing depth, the seismic velocity gradient first changes from negative to positive values in the southern hemisphere.

For a low concentration of HPEs in the crust (which implies a high concentration in the mantle), the seismic velocity variations in the upper part of the mantle are limited to a depth of about 300 km (lowE model, Figure 4a). A prominent low-velocity zone is visible in this case as a global feature to depths of about 300 km and extends locally in limited areas to 400 km depth (Figure 4d). For the case where nearly all HPEs are located in the crust (highE model), the crustal thickness pattern is visible to depths of 700 km and deeper. We note that in this case the phase transition in the pyroxene system (i.e., orthopyronexe to C2/c clinopyroxene) starts locally at a depth between 600 and 700 km (cf. detailed maps in Section 2 of the SI), leading to larger variations in the S wave velocity at that depth. This is also shown in the corresponding phase diagram (Figure 2c) and in the envelope around the averaged S wave velocity profile for the highE model (Figure 3c). The effect of this phase transition is also observed in the S wave velocity gradient (Figure 4f). At 700 km depth, the location of the phase transition is shown by the region of highest S wave velocity gradient around the northern lowlands. Compared to the lowE model, the highE case shows a less prominent LVZ in the upper part of the mantle that extends over a larger range of depths (starting from the base of the crust and extending up to about 900 km depth).

The S wave velocity variations of the mediumE case lie between those of the lowE and highE end-members, as shown in Figures 4a-4c. They follow the crustal thickness pattern up to 500 km depth, and the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Results of model calculations comparing three end-member models in terms of crustal heat-producing elements. Top panels: seismic velocity variations with respect to the mean value at three depth slices for the (a) lowE, (b) mediumE, and (c) highE models. Bottom panels: corresponding S wave velocity gradient with depth. The color scale is saturated for the largest S wave velocity gradient in panel (f) at 700 km depth. Negative values of the S wave velocity gradient show the locations of the low-velocity zones. Note that the chosen depths for each of the three models are different and were chosen to highlight the transition from where temperatures are controlled by crustal thickness variations to more uniform values in a convecting mantle. The mantle compositional model used here is TAY13 (G. J. Taylor, 2013). For orientation, a shaded-relief map of the surface topography is overlain on all maps. The location of the Interior exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) landing site and the fault system in Cerberus Fossae are indicated by the yellow triangle and the red circle, respectively, on the seismic velocity gradient maps.</figcaption></figure>

low-velocity zone extends locally up to depths of about 500–600 km. The LVZ is less pronounced than in the lowE model but more prominent than in the highE case, as observed in Figures 4d–4f.

In Figure 5, we plot 2D histograms of the S wave velocity variations for each of the three models (lowE, mediumE, and highE, respectively) as a function of depth using the TAY13 mantle composition. A similar comparison for the YOS20 composition is shown in Section 2 of the SI. Since the differences between the mantle and crustal seismic velocities are considerably larger than for regions located solely in the mantle, we plot the S wave variations starting at 250 km depth. Regions located at depths larger than 250 km are well within the mantle, and the differences shown in Figure 5 are entirely related to the mantle seismic velocities.

For the lowE case, where most of the HPEs are located in the mantle, lateral variations in seismic velocity below the crust are restricted to the uppermost 400 km (Figure 5a). On the other hand, for case highE, where a large amount of HPEs is in the crust, which, in turn, leads to a cold and thick lithosphere,
<figure><img src="image_7.png" /><figcaption>**Figure 5.** S wave velocity variations shown as 2D histograms as a function of depth for the (a) lowE, (b) mediumE, and (c) highE cases. The histograms were computed by equal-area binning the data at each depth between 250 km and core-mantle boundary, and the colors are proportional to the frequency of the S wave velocity values at each depth. The mantle compositional model used here is TAY13 (G. J. Taylor, 2013).</figcaption></figure>

prominent lateral variations in seismic velocity extend to deeper depths. In this case, a bimodal distribution of seismic velocities that is caused by crustal thickness variations can be observed to depths of about 700 km (Figure 5c). Starting at about 600 km depth, the phase transition from orthopyroxene to C2/c takes place. This phase transition leads also to a bimodal distribution of S wave velocity for the highE case that is associated with the fact that the phase transition extends over a range of about 300 km (cf. Figure S3 of the SI). The mediumE case lies between the lowE and highE end-members. For all cases, large seismic velocity variations are visible where the olivine to wadsleyite phase transition occurs at depths near 1,000–1,200 km.

### 3.3. Low-Velocity Zones

In Figure 6, we show a comparison of seismic velocity results for all 130 thermal evolution models using both the TAY13 and YOS20 compositions (260 models in total). For each case, the maximum seismic velocity variations (difference computed between maximum and minimum S wave velocity at each depth) can reach up to 10%, as shown in panel (a). Panel (b) shows that these large variations are located either at the depth where the transition from olivine to wadsleyite occurs or below the crust where temperature variations are the largest.

Decreasing the amount of HPEs in the crust increases the heat production in the mantle and hence gives rise to a more prominent low-velocity zone. This is the case for models that employ a thin crust and/or have significantly lower concentrations of HPEs than that deduced from surface GRS measurements (Hahn et al., 2011; G. J. Taylor et al., 2006). On the other hand, if most of the HPEs reside in the crust (thick crust and/or a higher amount of crustal heat production), the colder mantle leads to a less pronounced global LVZ (Figure 6c). For most models, the depth of a global LVZ (minimum S wave velocity gradient) is located in the uppermost 300 km (Figure 6d). Some of the models that use the TAY13 composition and employ a thick crust (average crustal thickness >62 km) containing more than 60% of the bulk HPEs, have the minimum $V_s$ gradient located close to the CMB. This LVZ located close to the CMB is observed in particular for cold temperatures (e.g., highE case in Figure 1). While a significant number of models that use the TAY13 composition present LVZs located close to the CMB, only a limited number of models using the YOS20 composition shows this behavior (cf. Figure S11 in the SI). Thus, the LVZ close to the CMB is affected by both temperature and composition, specifically by how Fe influences $V_s$ of each stable phase (see additional discussion in the SI, Section S1).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Summary of seismic velocity variations for the 3D thermal evolution models of Plesa, Padovan, et al. (2018):</figcaption></figure>

(a) maximum $V_s$ difference at each depth, calculated as max($\frac{V_{max} - V_{min}}{V_{avg}}$). Two compositional models have been used: G. J. Taylor (2013) (circles) and Yoshizaki and McDonough (2020) (stars).

(b) The depth at which the maximum difference shown in panel (a) was attained; (c) minimum $V_s$ gradient computed using the average profile; (d) the depth at which the minimum $V_s$ gradient shown in panel (c) was attained. Note the broken y-axis in panels (b and d) that shows the relevant depths.

### 3.4. S Wave Shadow Zones

Since many models show negative velocity gradients with depth, we test all models for the existence of seismic shadow zones between the InSight landing site and the locations of Marsquakes, with special consideration of the epicentral distances of events S0173a and S0235b (29° ± 3° and 26° ± 3°, respectively, Giardini et al., 2020).

Figure 7 shows ray-tracing calculations for the velocity models that are presented in Figures 3 and 4. A similar analysis for the YOS20 composition is shown in Figure S12 of the SI. For simplicity, our ray-tracing calculations make use of 1D seismic velocity profiles that represent averages between the InSight landing site and the presumed source in Cerberus Fossae. For the mediumE and highE models, no LVZ is observed for epicentral distances of up to 45°. The lowE case, however, shows two LVZs. The first LVZ is located at about 100 km depth and is responsible for the shadow zone observed for epicentral distance between 25° and 35°. The second low-velocity zone is located below 220 km depth and is only relevant for epicentral distances larger than 50°, which are not relevant for the seismic events discussed here.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Ray tracing and shadow zone analysis between the Interior exploration using Seismic Investigations, Geodesy and Heat Transport landing site and the direction of Cerberus Fossae. Three cases are presented for lowE (left column), mediumE (middle column), and highE (right column) models using the TAY13 composition. Panels (a-c): Density (black), *S* wave velocity (red), and *P* wave velocity (blue). Solid lines show depth profiles in radial coordinates, whereas the dashed lines are after the Flat Earth Transform (Müller, 1976). Horizontal black lines mark the Moho depth. Panels (d-f): Raypaths for *P* (blue) and *S* (red), for a source depth of 50 km. Pale thin red and blue lines show the *S*- and *P*-raypaths, respectively, over the entire distance range between 0° and 45°. Thick red and blue lines are *S*- and *P*-rays, respectively, to the epicentral distances of events S0173a (29°) and S0235b (26°) that are marked by black horizontal lines. Panels (g-i): Travel time curves computed every 0.1° of epicentral distance, minus a linear trend of 6 s per degree distance. Blue and red circles mark *P* and *S* travel times, respectively, computed at the epicentral distances of events S0173a (29°) and S0235b (26°). Gaps in the travel time curves denote shadow zones. Small negative velocity gradients in the crust cause multipathing and thus additional *S* wave arrivals in panels (g and i). Panels (j-l): Ray parameter (the invariant of Snell's law) required to arrive at given distances, for *S* (red) and *P* (blue).</figcaption></figure>

We perform the ray-tracing analysis for all seismic velocity models and discard those with shadow zones at the epicentral distances of events S0173a and S0235b (i.e., 25°-35°). From a total of 260 models (130 thermal evolution models using two compositions, i.e., TAY13 and YOS20), we find that eight models using the TAY13 and YOS20 compositions had *S* wave shadow zones in the relevant distance range (i.e., 25°-35°) when the source depth is located at 50 km. These models contained less than 20% of the total amount of
HPEs in the crust, employed a low mantle reference viscosity (i.e., η<sub>ref</sub> < 5 × 10<sup>20</sup> Pa s), or had a lower ther-
mal conductivity of the crust compared to other cases (2 W m<sup>-1</sup> K<sup>-1</sup> vs. 3 W m<sup>-1</sup> K<sup>-1</sup>). For source depths
located at 100 and 150 km, the models that produced shadow zones showed similar parameters as for a
source depth of 50 km. However, as the source depth increased, the number of models presenting shadow
zones was found to decrease. Two cases, for which shadow zones are present, had unrelated parameters
with thicker mantles and hence smaller cores (mantle thickness *D* ≥ 1,700 km). A detailed analysis of all
models is presented in Section S6 and Tables S9–S11 of the SI. Interestingly, all models that show shadow
zones have a relatively thin lithosphere (lithosphere thickness thinner than 450 km on average) and high
average thermal gradients in the lithosphere (larger than 3.5 K/km). This is shown in Figure 8 by the black
symbols. The average thermal gradient has been calculated by computing a volume-averaged gradient up to
a depth where the thermal gradient drops below 1 K/km. As a comparison, Figure 8 shows also the best fit
models of Plesa, Padovan, et al. (2018) that match a large number of geological and geophysical constraints
as colored symbols.

We note that additional models, many of them employing the YOS20 composition, do not have distinct shadow zones but defocus the wavefield, such that it is unlikely to observe clear arrivals. These models have similar parameters to the models that have distinct shadow zones (i.e., low amount of HPEs in the crust and a low reference viscosity). In addition, models with a low activation volume (≤6 cm³/mol) and hence a low pressure-dependence of the viscosity, and/or a thick mantle (D = 1,900 km) are more prone to defocus the wavefield. To be able to provide a quantitative analysis of our models, we calculate the number of S-rays that arrive at epicentral distances between 25° and 30° relative to the number of S-rays arriving at the same epicentral distance range in a constant velocity model. The distance range for ray takeoff was chosen between 5° and 65°, with a ray takeoff angle resolution of 0.02°. In Figure 9, we show the percentage of S-rays arriving between 25° and 30° epicentral distance relative to a constant velocity model for TAY13 (Figures 9a, 9c, and 9e) and YOS20 (Figures 9b, 9d, and 9f) for source depths of 50, 100, and 150 km. The color of each symbol indicates the amount of HPEs in the crust relative to the bulk. We also mark the best fit models of Plesa, Padovan, et al. (2018) by magenta and blue colors for TAY13 and YOS20, respectively. A complete list of models and their analysis is available in the supplementary material.

While the number of *S*-rays arriving at 25°-30° epicentral distance range increases with increasing source depth, models with large average lithospheric temperature gradients and small amount of crustal HPEs stand out showing less than 20% *S* waves arrivals compared to a constant velocity model, even for a source depth located at 150 km.

4. Discussion

The distribution of seismic velocities, in particular in the lithosphere, is sensitive to the amount of HPEs located in the crust. Additionally, the magnitude (i.e., how well pronounced a LVZ is) and extent (i.e., how deep a LVZ reaches and whether its distribution is global or local) of a LVZ depends on the thermal structure of the lithosphere and hence on the abundances of HPEs in the mantle and crust and the 3D variations in crustal thickness. Additional factors such as the presence of partial melts and volatiles, which were not considered here, could further affect the distribution of seismic velocities, as it is the case on the Earth (e.g., Golos et al., 2020). However, we note that in contrast to Earth such effects may be more limited on Mars. Although volcanic activity as recent as several Myr has been recorded in Tharsis and Elysium (Hauber et al., 2011; Neukum et al., 2004; Vaucher et al., 2009), partial melt is expected to be restricted to limited areas in these two large volcanic provinces. A recent study by Kedar et al. (2021) has analyzed a handful of low frequency Marsquakes, a class of seismic events recorded by InSight that have similar characteristics to quakes produced by volcanic processes on Earth. Although a low-viscosity magma and high effusion rates are required to explain the observed signature of the low frequency events, these physical parameters are still within bounds that have been suggested for Cerberus Fossae. Thus, a magmatic origin of the low frequency events cannot be excluded (Kedar et al., 2021). Analyses of Martian meteorites suggest that the Amazonian Martian mantle contains about 14–72 ppm water that is heterogeneously distributed (McCubbin et al., 2016). These values indicate a drier mantle compared to the Earth, drier even than depleted MORB mantle sources (McCubbin et al., 2016). However, the exact distribution of volatiles in the Martian mantle is somewhat controversial and remains poorly constrained (Udry et al., 2020).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The temperature gradient averaged from the surface to the average depth of the lithosphere versus the average lithospheric thickness for the entire data set of 130 models. Black symbols show cases that produce shadow zones for the TAY13 (Panels a, c, and e) and YOS20 (Panels b, d, and f) compositions, when the source is located at 50 km (Panels a and b), 100 km (Panels c and d), and 150 km (Panels e and f) depth. Similar to Figure 6, the size of the symbols indicates the crustal thickness (i.e., small symbols represent cases with a thin crust, while large symbols indicate cases with a thick crust). Best fit models of Plesa, Padovan, et al. (2018) are shown by symbols filled with orange and cyan colors for the TAY13 and YOS20 compositions, respectively.</figcaption></figure>

Our results show that models with a thin lithosphere (thinner than 400–450 km, on average) produce a strong LVZ that leads to an *S* wave shadow zone, which is incompatible with clear *S* wave phases that are present for the events S0173a and S0235b. For models with an average lithosphere thickness that exceeds 500 km, only a weak LVZ exists. We note, however, that although weak, this zone extends to a greater depth compared to cases where a more pronounced LVZ is present. Our results also show that even though the average lithosphere thickness may exceed 500 km, the lithosphere may be locally thinner below regions
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The percentage of *S*-rays arriving in the 25°-30° epicentral distance range calculated relative to the number of *S*-rays arriving in the same epicentral distance range in a constant velocity model versus the average lithospheric temperature gradient. The left panels (a, c, and e) are for the TAY13 composition and the right panels (b, d, and f) are for the YOS20 composition. The source depth of the seismic events has been set at 50 km (Panels a and b), 100 km (Panels c and d), and 150 km (Panels e and f). The color shows the crustal amount of HPEs relative to the bulk amount of HPEs. Symbols with magenta and blue border represent the best fit models of Plesa, Padovan, et al. (2018). Note that the y-axis scale differs for each row.</figcaption></figure>

of thick crust and thicker below areas covered by a thin crust. Nevertheless, a thick average lithosphere is indicative of a crust highly enriched in HPEs if the bulk amount of HPEs is similar to those estimated by geochemical models (Lodders & Fegley, 1997; G. J. Taylor, 2013; Wänke & Dreibus, 1994; Yoshizaki & McDonough, 2020).
A crust highly enriched in HPEs could significantly affect the distribution of seismic velocities in the lithosphere, where a pattern similar to the crustal thickness pattern may persist up to 400–600 km depth. In addition, local low-velocity regions or global LVZs may be present over a range of depths. Indeed, a shallow global LVZ disappears first in the southern hemisphere and eventually entirely at depths larger than 700 km.

We note that the constraints provided by the absence of a shadow zone between the InSight landing site and Cerberus Fossae exclude thermal evolution models that were previously discarded based on a number of geological, petrological, and geophysical constraints (Plesa, Padovan, et al., 2018; Smrekar et al., 2019). Indeed, previously identified best fit models remain valid when using the shadow-zone constraint discussed in this study. The shadow-zone constraint requires that the crust contains a significant fraction of total HPEs of Mars, as previously suggested by models that reproduce the present-day elastic lithosphere thickness (Plesa, Padovan, et al., 2018; Thiriet et al., 2018). While the absence of shadow zones is less able to constrain the core size compared to geodetic observations, which require a relatively large core (Khan et al., 2018; Plesa, Padovan, et al., 2018; Rivoldini et al., 2011), models employing a small core (<1,700 km) are more prone to produce shadow zones and to defocus the wavefield compared to models using larger cores. This is due to the fact that a smaller core increases the thickness of the mantle. A larger silicate volume that contains a higher amount of HPEs will affect the thermal profiles by decreasing the thickness of the lithosphere with increasing mantle thickness, which, in turn, may lead to shadow zones.

One of the InSight goals is to determine the crustal thickness at the landing site providing thereby a valuable anchor point for global crustal thickness models (Knapmeyer-Endrun et al., 2020; Smrekar et al., 2019). A global crustal thickness model together with the presence or absence of a LVZ will help to put constraints on the distribution of HPEs between the mantle and crust for a given bulk composition model (e.g., G. J. Taylor, 2013; Wänke & Dreibus, 1994; Yoshizaki & McDonough, 2020). A weak LVZ that extends locally to depths of 600 km or deeper would be representative of cases where more than half of the HPEs need to be located in the crust. For compositional models such as Yoshizaki and McDonough (2020) and Lodders and Fegley (1997), which contain a higher amount of HPEs, an equivalent of more than 57% and 67% of the total amount of HPEs, respectively, would need to be located in the crust. Additionally, these scenarios would be characterized by seismic velocities variations that closely follow the crustal thickness variations to depths of 500 km or more. On the other hand, cases that contain less than 20% of the total amount of HPEs in their crust produce a strong LVZ. In these cases, the LVZ is restricted to the uppermost 300–400 km, with only limited low-velocity areas still present at a depth of 400 km. We note, however, that for compositional models with a higher amount of HPEs, such as Yoshizaki and McDonough (2020) and Lodders and Fegley (1997), a pronounced LVZ could be expected for a crustal HPEs content of up to 32% and 47% of the total amount of HPEs, respectively.

Depending on its magnitude and extent, a LVZ at the CMB, as observed in some of the models presented here, could affect the travel times of core-related phases and should be considered in future studies. However, we note that the LVZ at the CMB is not necessarily global, it appears predominantly for TAY13 composition, when the stability field of garnet significantly extends at high pressure (see Figure 2), and the negative seismic velocity gradient is rather small (>-0.00175 s⁻¹). This could lead to a few percent difference between travel times of core-reflected phases for models with and without a LVZ at the CMB. Nevertheless, for such an effect to be observed, seismic events with distinguishable core phases are necessary (see Panning et al., 2017, and references therein for a discussion of the detectability of core phases).

A ray-tracing algorithm, as we have used here, provides a first-order assessment on whether a shadow zone is present on the path between the InSight landing site and the Cerberus Fossae area, which has been proposed as the source location for two of the high-quality events. For these events, epicentral distances were estimated by comparing the observed P–S arrival time differences with a suite of a priori velocity models and a high signal-to-noise ratio allowed for the determination of the back azimuth from polarization (Giardini et al., 2020). Cerberus Fossae is a geologically plausible source region, since it is a young volcanic area and suspected to be tectonically very active (J. Taylor et al., 2013). Our analysis suggests that models containing less than 20% of the total amount of HPEs in the crust have low-velocity zones incompatible with current InSight observations for the S0173a and S0235b events. Additionally, models with a low reference viscosity (i.e., η<sub>ref</sub>~5 × 10²⁰ to 10²⁰ Pa s), a low crustal thermal conductivity, and a smaller core (<1,700 km) are more likely to produce shadow zones or strongly defocus the wavefield. However, we note that the source depth
of the seismic event plays a role in whether a model shows a distinct shadow zone. A shadow zone is more likely to be present for shallow events (source depths of 50–100 km) than for deeper ones (source depths deeper than 150 km). Nevertheless, for cases that show distinct shadow zones for shallow seismic sources, the number of S-rays that reach the relevant epicentral distance is limited, even if the seismic event is located at 150 km depth. A recent study by Brinkman et al. (2021) that investigated the focal mechanism of the events in Cerberus Fossae, finding an extensional regime to be representative for the two high-quality events, also performed a first attempt to constrain the source depth of several Marsquakes. However, the uncertainties for the source depth of the events still remain large (Brinkman et al., 2021). A depth estimate for the seismic events could further help to place constraints on the thermal state and thickness of the lithosphere. In the absence of plate tectonics, Marsquakes are expected to occur at shallower depths compared to subduction zone earthquakes. However, for a thick and cold lithosphere, the seismogenic layer thickness may reach depths of up to 450 km (Plesa, Knapmeyer, et al., 2018).

We note that the calculations presented here are based on averaged velocity profiles computed on a regional scale between InSight and Cerberus Fossae. This approach is justified by the proximity of the Cerberus Fossae events to the InSight landing site that makes lateral variations in crustal and thermal lithosphere thicknesses negligible. For larger epicentral distance, however, the 3D thermal structure needs to be combined with seismic wave propagation simulations to determine to what extent local features such as regional low-velocity zones could affect the arrival of S waves at the InSight location.

## 5. Conclusions

We have used the entire data set of 3D thermal evolution models of Plesa, Padovan, et al. (2018) to investigate the effects of the distribution of HPEs and mantle composition on the variations of seismic velocities in the interior of Mars. The distribution of seismic velocities in the lithosphere follows the crustal thickness pattern to depths of 400 km and deeper if the crust contains a large percentage of the bulk HPEs (>50%). The distribution of HPEs between the mantle and crust and hence the thermal state of the lithosphere strongly affects the depth to which a LVZ may be present and whether this layer is well pronounced or weak. In addition, a LVZ may not be global but rather be present only locally.

Future seismic data from InSight will help us to constrain the distribution of HPEs between mantle and crust for a given bulk amount of HPEs (e.g., Lodders & Fegley, 1997; G. J. Taylor, 2013; Treiman et al., 1986; Wänke & Dreibus, 1994; Yoshizaki & McDonough, 2020). The presence or absence of a low-velocity zone (Giardini et al., 2020) will place constraints on the thermal state of the lithosphere that, in turn, is directly linked to the distribution of HPEs between the mantle and crust. While, so far, surface waves have not been identified in the seismic events observed by InSight, if detected in future seismic data, analysis of surface-wave dispersion could place additional constraints on the seismic velocities of the crust and lithosphere. Together with robust estimates of the crustal thickness at InSight location (Knapmeyer-Endrun et al., 2020) that will constrain crustal thickness models (Wieczorek et al., 2019), the amount of crustal HPEs will help to determine whether the surface abundance of HPEs derived from GRS measurements (Hahn et al., 2011; G. J. Taylor et al., 2006) is representative for the entire crust.

## Data Availability Statement

Data sets containing 3D seismic velocity data are available on Figshare https://doi.org/10.6084/m9.figshare.12625703.v1

## References

Bagheri, A., Khan, A., Al-Attar, D., Crawford, O., & Giardini, D. (2019). Tidal response of mars constrained from laboratory-based viscoelastic dissipation models and geophysical data. *Journal of Geophysical Research: Planets*, **124**, 2703–2727. https://doi.org/10.1029/2019JE006015

Banerdt, W. B., Smrekar, S. E., Banfield, D., Giardini, D., Golombek, M., Johnson, C. L., et al. (2020). Initial results from the InSight mission on Mars. *Nature Geoscience*, **13**, 183–189.

Bissig, F., Khan, A., Van Driel, M., Stähler, S. C., Giardini, D., Panning, M., et al. (2018). On the detectability and use of normal modes for determining interior structure of Mars. *Space Science Reviews*, **214**(8), 114.

## Acknowledgments

The authors thank Hao Hu and an anonymous reviewer for their thoughtful comments that helped to improve an earlier version of the manuscript. The authors also thank the Editor Laurent Montesi for handling the manuscript and for his advice. All the parameters used in the models and their outcomes are listed in the online Supporting Information. A.-C. Plesa gratefully acknowledges the financial support and endorsement from the DLR Management Board Young Research Group Leader Program and the Executive Board Member for Space Research and Technology. E. Bozdağ was supported by NASA Participating Scientist Program grant 80NSSC18K1680. A. Rivoldini acknowledges the support of the Belgian PRODEX program managed by the European Space Agency in collaboration with the Belgian Federal Science Policy Office. S. M. McLennan was supported by NASA Participating Scientist Program grant 80NSSC18K1622. S. Padovan acknowledges support by the DFG within the Research Unit FOR 2440 “Matter Under Planetary Interior Conditions.” M. A. Wieczorek was supported by a grant from the French space agency (CNES). Additional simulations used in this work were performed on the supercomputer ForHLR funded by the Ministry of Science, Research and the Arts Baden-Württemberg and by the Federal Ministry of Education and Research. This is InSight Publication No. 113.
Bozdağ, E., Ruan, Y., Metthez, N., Khan, A., Leng, K., van Driel, M., et al. (2017). Simulations of seismic wave propagation on Mars. *Space Science Reviews*, **211**(1–4), 571–594.

Brinkman, N., Stähler, S. C., Giardini, D., Schmelzbach, C., Khan, A., Jacob, A., et al. (2021). First focal mechanisms of Marsquakes. *Journal of Geophysical Research: Planets*, **126**, e2020JE006546. https://doi.org/10.1029/2020JE006546

Clinton, J. F., Ceylan, S., van Driel, M., Giardini, D., Stähler, S. C., Böse, M., et al. (2021). The Marsquake catalogue from InSight, sols 0-478. *Physics of the Earth and Planetary Interiors*, **310**, 106595. https://doi.org/10.1016/j.pepi.2020.106595

Connolly, J. (2009). The geodynamic equation of state: What and how. *Geochemistry, Geophysics, Geosystems*, **10**, Q1004. https://doi.org/10.1029/2009GC002540

Cottaar, S., Heister, T., Rose, I., & Unterborn, C. (2014). BurnMan: A lower mantle mineral physics toolkit. *Geochemistry, Geophysics, Geosystems*, **15**, 1164–1179. https://doi.org/10.1002/2013GC005122

Cottaar, S., Heister, T., Rose, I., & Unterborn, C. (2016). *BurnMan*. Computational infrastructure for geodynamics. https://doi.org/10.5281/zenodo.546210

Finkelstein, G. J., Dera, P. K., & Duffy, T. S. (2015). Phase transitions in orthopyroxene (En90) to 49 GPa from single-crystal X-ray diffraction. *Physics of the Earth and Planetary Interiors*, **244**, 78–86. https://doi.org/10.1016/j.pepi.2014.10.009

Folkner, W. M., Dehant, V., Le Maistre, S., Yseboodt, M., Rivoldini, A., Van Hoolst, T., et al. (2018). The rotation and interior structure experiment on the InSight mission to Mars. *Space Science Reviews*, **214**(5), 100.

Franck, S., & Kowalle, G. (1994). Seismic velocity models for an internally asymmetric Mars. *Earth, Moon, and Planets*, **65**(3), 277–290.

Giardini, D., Lognonné, P., Banerdt, W. B., Pike, W. T., Christensen, U., Ceylan, S., et al. (2020). The seismicity of Mars. *Nature Geoscience*, **13**(3), 205–212.

Golombek, M., Warner, N., Grant, J., Hauber, E., Ansan, V., Weitz, C., et al. (2020). Geology of the InSight landing site on Mars. *Nature Communications*, **11**(1), 1–11.

Golos, E. M., Fang, H., & van der Hilst, R. D. (2020). Variations in seismic wave speed and VP/VS ratio in the North American lithosphere. *Journal of Geophysical Research: Solid Earth*, **125**, e2020JB020574. https://doi.org/10.1029/2020JB020574

Gudkova, T., & Zharkov, V. (2004). Mars: Interior structure and excitation of free oscillations. *Physics of the Earth and Planetary Interiors*, **142**(1–2), 1–22. https://doi.org/10.1016/j.pepi.2003.10.004

Hahn, B. C., McLennan, S. M., & Klein, E. C. (2011). Martian surface heat production and crustal heat flow from Mars Odyssey Gamma-Ray spectrometry. *Geophysical Research Letters*, **38**, L14203. https://doi.org/10.1029/2011GL047435

Hauber, E., Broż, P., Jagert, F., Jodłowski, P., & Platz, T. (2011). Very recent and wide-spread basaltic volcanism on Mars. *Geophysical Research Letters*, **38**, L10201. https://doi.org/10.1029/2011GL047310

Holland, T. J. B., & Powell, R. (2011). An improved and extended internally consistent thermodynamic dataset for phases of petrological interest, involving a new equation of state for solids. *Journal of Metamorphic Geology*, **29**(3), 333–383.

InSight Marsquake Service. (2020). *Mars seismic catalogue, InSight mission; V3 2020-07-01*. ETHZ, IPGP, JPL, ICL, ISAE-Supaero, MPS, University of Bristol. https://doi.org/10.12686/a8

Johnson, C. L., Mittelholz, A., Langlais, B., Russell, C. T., Ansan, V., Banfield, D., et al. (2020). Crustal and time-varying magnetic fields at the InSight landing site on Mars. *Nature Geoscience*, **13**(3), 199–204.

Kedar, S., Panning, M. P., Smrekar, S. E., Stähler, S. C., King, S. D., Golombek, M. P., et al. (2021). Analyzing low frequency seismic events at Cerberus Fossae as long period volcanic quakes. *Journal of Geophysical Research: Planets*, **126**, e2020JE006518. https://doi.org/10.1029/2020JE006518

Khan, A., Liebske, C., Rozel, A., Rivoldini, A., Nimmo, F., Connolly, J. A. D., et al. (2018). A geophysical perspective on the bulk composition of Mars. *Journal of Geophysical Research: Planets*, **123**, 575–611. https://doi.org/10.1002/2017JE005371

Knapmeyer, M. (2004). TTBox: A MatLab toolbox for the computation of 1D teleseismic travel times. *Seismological Research Letters*, **75**(6), 726–733.

Knapmeyer, M. (2005). Numerical accuracy of travel-time software in comparison with analytic results. *Seismological Research Letters*, **76**(1), 74–81.

Knapmeyer-Endrun, B., Bissig, F., Compaire, N., Joshi, R., Garcia, R., Khan, A., et al. (2020). *Seismic constraints on the crustal structure of Mars from InSight receiver functions*. Paper presented at 51st Lunar and Planetary Science Conference (2326), 1914.

Lodders, K., & Fegley, B. (1997). An oxygen isotope model for the composition of Mars. *Icarus*, **126**(2), 373–394.

Lognonné, P., Banerdt, W. B., Giardini, D., Pike, W., Christensen, U., Laudet, P., et al. (2019). SEIS: Insight's seismic experiment for internal structure of Mars. *Space Science Reviews*, **215**(1), 12.

Lognonné, P., Banerdt, W. B., Pike, W., Giardini, D., Christensen, U., Garcia, R., et al. (2020). Constraints on the shallow elastic and anelastic structure of Mars from InSight seismic data. *Nature Geoscience*, **13**(3), 213–220.

McCubbin, F. M., Boyce, J. W., Srinivasan, P., Santos, A. R., Elardo, S. M., Filiberto, J., et al. (2016). Heterogeneous distribution of H₂O in the Martian interior: Implications for the abundance of H₂O in depleted and enriched mantle sources. *Meteoritics & Planetary Science*, **51**(11), 2036–2060.

Mocquet, A., & Menvielle, M. (2000). Complementarity of seismological and electromagnetic sounding methods for constraining the structure of the Martian mantle. *Planetary and Space Science*, **48**(12–14), 1249–1260.

Mocquet, A., Vacher, P., Grasset, O., & Sotin, C. (1996). Theoretical seismic models of Mars: The importance of the iron content of the mantle. *Planetary and Space Science*, **44**(11), 1251–1268.

Müller, G. (1976). Earth-flattening approximation for body waves derived from geometric ray theory-improvements, corrections and range of applicability. *Journal of Geophysics*, **42**(1), 429–436.

Neukum, G., Jaumann, R., Hoffmann, H., Hauber, E., Head, J. W., Basilevsky, A. T., et al. (2004). Recent and episodic volcanic and glacial activity on Mars revealed by the high resolution stereo camera. *Nature*, **432**, 971–979.

Nimmo, F., & Tanaka, K. (2005). Early crustal evolution of Mars. *Annual Review of Earth and Planetary Sciences*, **33**(1), 133–161. https://doi.org/10.1146/annurev.earth.33.092203.122637

Panning, M. P., Lognonné, P., Banerdt, W. B., Garcia, R., Golombek, M., Kedar, S., et al. (2017). Planned products of the Mars structure service for the InSight mission to Mars. *Space Science Reviews*, **211**(1–4), 611–650.

Plesa, A.-C., Grott, M., Tosi, N., Breuer, D., Spohn, T., & Wieczorek, M. A. (2016). How large are present-day heat flux variations across the surface of Mars? *Journal of Geophysical Research: Planets*, **121**, 2386–2403. https://doi.org/10.1002/2016JE005126

Plesa, A.-C., Knapmeyer, M., Golombek, M. P., Breuer, D., Grott, M., Kawamura, T., et al. (2018). Present-day Mars' seismicity predicted from 3-D thermal evolution models of interior dynamics. *Geophysical Research Letters*, **45**, 2580–2589. https://doi.org/10.1002/2017GL076124

Plesa, A.-C., Padovan, S., Tosi, N., Breuer, D., Grott, M., Wieczorek, M., et al. (2018). The thermal state and interior structure of Mars. *Geophysical Research Letters*, **45**, 12–198. https://doi.org/10.1029/2018GL080728
Rivoldini, A., Van Hoolst, T., Verhoeven, O., Mocquet, A., & Dehant, V. (2011). Geodesy constraints on the interior structure and composition of Mars. *Icarus*, **213**(2), 451–472. https://doi.org/10.1016/j.icarus.2011.03.024

Smrekar, S. E., Lognonné, P., Spohn, T., Banerdt, W. B., Breuer, D., & Christensen, U., et al. (2019). Pre-mission InSights on the Interior of Mars. *Space Science Reviews*, **215**(1), 3.

Sohl, F., & Spohn, T. (1997). The interior structure of Mars: Implications from SNC meteorites. *Journal of Geophysical Research*, **102**(E1), 1613–1635.

Spiga, A., Banfield, D., Teanby, N. A., Forget, F., Lucas, A., Kenda, B., et al. (2018). Atmospheric science with InSight. *Space Science Reviews*, **214**(7), 109.

Spohn, T., Grott, M., Smrekar, S., Knollenberg, J., Hudson, T., Krause, C., et al. (2018). The heat flow and physical properties package (HP') for the InSight mission. *Space Science Reviews*, **214**(5), 96.

Stixrude, L., & Lithgow-Bertelloni, C. (2005). Thermodynamics of mantle minerals—I. Physical properties. *Geophysical Journal International*, **162**(2), 610–632.

Stixrude, L., & Lithgow-Bertelloni, C. (2011). Thermodynamics of mantle minerals—II. Phase equilibria. *Geophysical Journal International*, **184**(3), 1180–1213.

Taylor, G. J. (2013). The bulk composition of Mars. *Geochemistry*, **73**(4), 401–420.

Taylor, G. J., Boynton, W., Brückner, J., Wänke, H., Dreibus, G., Kerry, K., et al. (2006). Bulk composition and early differentiation of Mars. *Journal of Geophysical Research*, **112**, E03S10. https://doi.org/10.1029/2005JE002645

Taylor, J., Teanby, N. A., & Wookey, J. (2013). Estimates of seismic activity in the Cerberus Fossae region of Mars. *Journal of Geophysical Research: Planets*, **118**, 2570–2581. https://doi.org/10.1002/2013JE004469

Taylor, S. R., & McLennan, S. (2009). *Planetary crusts: Their composition, origin and evolution*. Cambridge University Press.

Thiriet, M., Michaut, C., Breuer, D., & Plesa, A.-C. (2018). Hemispheric dichotomy in lithosphere thickness on mars caused by differences in crustal structure and composition. *Journal of Geophysical Research: Planets*, **123**, 823–848. https://doi.org/10.1002/2017JE005431

Treiman, A. H., Drake, M. J., Janssens, M.-J., Wolf, R., & Ebihara, M. (1986). Core formation in the Earth and Shergottite Parent Body (SPB): Chemical evidence from basalts. *Geochimica et Cosmochimica Acta*, **50**(6), 1071–1091. https://doi.org/10.1016/0016-7037(86)90389-3

Udry, A., Howarth, G. H., Herd, C., Day, J. M., Lapen, T. J., & Filiberto, J. (2020). What Martian meteorites reveal about the interior and surface of Mars (Vol. 125). Wiley Online Library.

Vaucher, J., Baratoux, D., Mangold, N., Pinet, P., Kurita, K., & Grégoire, M. (2009). The volcanic history of central Elysium Planitia: Implications for Martian magmatism. *Icarus*, **204**(2), 418–442. https://doi.org/10.1016/j.icarus.2009.06.032

Verhoeven, O., Rivoldini, A., Vacher, P., Mocquet, A., Choblet, G., Menvielle, M., et al. (2005). Interior structure of terrestrial planets: Modeling Mars' mantle and its electromagnetic, geodetic, and seismic properties. *Journal of Geophysical Research*, **110**, E04009. https://doi.org/10.1029/2004JE002271

Wänke, H., & Dreibus, G. (1994). Chemistry and accretion of Mars. *Philosophical Transactions of the Royal Society London*, **A349**, 2134–2137.

Wieczorek, M. A., Beuthe, M., Rivoldini, A., & Van Hoolst, T. (2019). Hydrostatic interfaces in bodies with nonhydrostatic lithospheres. *Journal of Geophysical Research: Planets*, **124**, 1410–1432. https://doi.org/10.1029/2018JE005909

Yoshizaki, T., & McDonough, W. F. (2020). The composition of Mars. *Geochimica et Cosmochimica Acta*, **273**, 137–162.

Zheng, Y., Nimmo, F., & Lay, T. (2015). Seismological implications of a lithospheric low seismic velocity zone in Mars. *Physics of the Earth and Planetary Interiors*, **240**, 132–141.