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

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.

aims at understanding the interior structure and present-day thermal state of Mars and can provide valuable constraints for the planet's thermochemical history (Smrekar et al., 2019).
Since February 2019, the seismometer has continuously recorded Mars' seismic activity 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 , 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 . 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 computations 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 . 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.

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 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 (medi-umE) 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).

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

Table 1 Summary of Crustal Thickness and Crustal Amount of HPEs Used in This Study
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(Cottaar et al., , 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) (2011), only a high-pressure phase was available (akimotoite).
We have also assumed that all iron is present as Fe 2+ . Even with these simplifications, the resulting average density of the crust is 3,150 kg m −3 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.

Ray-Tracing Models
For the 1D ray tracing, we use the TTBOX package, which was validated and verified by Knapmeyer (2004Knapmeyer ( , 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/defocusing 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.

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.

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 PLESA ET AL. 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 PLESA ET AL. 10.1029/2020JE006755 6 of 19 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 be-  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.

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 PLESA ET AL.

10.1029/2020JE006755
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, me-diumE, 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, PLESA ET AL.
10.1029/2020JE006755 9 of 19 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.

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).

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.
PLESA ET AL. 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 PLESA ET AL.

10.1029/2020JE006755
12 of 19 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).
HPEs in the crust, employed a low mantle reference viscosity (i.e., η ref < 5 × 10 20 Pa s), or had a lower thermal conductivity of the crust compared to other cases (2 W m −1 K −1 vs. 3 W m −1 K −1 ). 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 3 /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.

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).
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 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.
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).

10.1029/2020JE006755
15 of 19 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.
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 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 −1 ). 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 . 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., η ref ∼5 × 10 20 to 10 20 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.

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  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.