# Constraints on Lateral Variations of Martian Crustal Thickness From Seismological and Gravity Field Measurements

Mélanie Drilleau, Henri Samuel, Raphaël F Garcia, Attilio Rivoldini, Clément

Perrin, Mark Wieczorek, Philippe Lognonné, William B Banerdt

► To cite this version:

Mélanie Drilleau, Henri Samuel, Raphaël F Garcia, Attilio Rivoldini, Clément Perrin, et al. Constraints on Lateral Variations of Martian Crustal Thickness From Seismological and Gravity Field Measurements. Geophysical Research Letters, 2024, 51 (4), pp.e2023GL105701. 10.1029/2023gl105701. hal-04467414

HAL Id: hal-04467414

https://hal.science/hal-04467414v1

Submitted on 20 Feb 2024

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
# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2023GL105701

**Key Points:**

* We infer lateral variation of Martian crustal thickness using InSight seismic data combined with topography and gravity field modeling

* Marsquakes are relocated and maps of the seismicity recorded by InSight are proposed

* Two families of models are found, leading to different interpretations of Mars' thermochemical evolution, with hot and cold scenarios

**Supporting Information:**

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

**Correspondence to:**

M. Drilleau
melanie.drilleau@isae-supaero.fr

**Citation:**

Drilleau, M., Samuel, H., Garcia, R. F., Rivoldini, A., Perrin, C., Wieczorek, M., et al. (2024). Constraints on lateral variations of Martian crustal thickness from seismological and gravity field measurements. *Geophysical Research Letters*, *51*, e2023GL105701. https://doi.org/10.1029/2023GL105701

Received 27 JULY 2023

Accepted 31 JAN 2024

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

## Constraints on Lateral Variations of Martian Crustal Thickness From Seismological and Gravity Field Measurements

Mélanie Drilleau¹ ID, Henri Samuel² ID, Raphaël F. Garcia¹ ID, Attilio Rivoldini³ ID, Clément Perrin⁴ ID, Mark Wieczorek² ID, Philippe Lognonné² ID, and William B. Banerdt⁵ ID

¹ISAE-SUPAERO, Université de Toulouse, Toulouse, France, ²Université Paris Cité, Institut de physique du globe de Paris, CNRS, Paris, France, ³Royal Observatory of Belgium, Brussels, Belgium, ⁴Nantes Université, Université d'Angers, Le Mans Université, CNRS, UMR 6112, Laboratoire de Planétologie et Géosciences, UAR 3281, Observatoire des Sciences de l'Univers de Nantes Atlantique, Nantes, France, ⁵Jet Propulsion Laboratory, California Institute of Technology, Pittsburgh, PA, USA

**Abstract** Using body wave arrival times from 31 seismic events recorded on Mars by the InSight mission, combined with topography and gravity field modeling, we constrained lateral variations of crustal thickness through a Bayesian inversion approach. The parameterization of the seismic structure relies on quantities that influence the thermochemical evolution of Mars, enabling the seismic velocities and densities in the different planetary envelopes to be consistently linked through common physical assumptions. Compared to a 1D structure, models with lateral variations of crustal thickness show two possible interpretations of the thermal evolution of Mars, with either a hot or cold scenario at the present-day. We found the hot scenario to be more compatible with InSight's radiotracking data and the tidal Love number. We relocated the marsquakes and derived maps of seismicity recorded by InSight, which is mostly located along or North of the boundary between the Northern lowlands and the Southern highlands.

**Plain Language Summary**

Thanks to the seismometer of the InSight mission, which recorded ground vibration measurements emanating from marsquakes and meteorite impacts during almost 4 years, the 1D interior structure of the crust, mantle, and core have been revealed. These models are mainly based on the assumption that the crustal thickness is similar everywhere on Mars. However, both the InSight lander and the most of the seismic events are located between the Northern lowlands and the Southern highlands where the crustal thickness varies widely, which can bias the interpretation of a 1D crustal model. In this study, combining the InSight seismic data with other independent geophysical measurements (topography and gravimetry data), we investigated to what extent the interior structure models are modified if lateral variations of crustal thickness are considered. Our results show that two different interpretations of Mars' thermal history can be considered, with either a hot or a cold scenario. We assessed the compatibility of our results with independent observations of oscillation of Mars' rotational axis, and found that the hot scenario is most likely.

## 1. Introduction

During almost 4 years of operation on Mars, the seismometer of the InSight mission (Banerdt et al., 2020; Lognonné et al., 2019), deployed in the Northern lowlands close to both the crustal dichotomy boundary between the Southern highlands and the Elysium volcanic provinces to the North, has recorded more than 1,300 seismic events (InSight Marsquake Service, 2023). Among these events several marsquakes show energy predominantly below 1 Hz, for which the identification of body wave arrivals was possible. Using these marsquakes, and several meteorite impacts (Garcia et al., 2022; Posiolova et al., 2022), 1D average models of the crust (Knapmeyer-Endrun et al., 2021; Li et al., 2022; Lognonné et al., 2020), mantle (Drilleau et al., 2022; Durán et al., 2022; Khan et al., 2021) and core (Irving et al., 2023; Samuel et al., 2023; Stähler et al., 2021) have been deduced under the assumption of a spherically symmetric planet. However, seismic wavefield disturbances from a 3D shallow structure would undeniably complicate interpretations based on 1D radial models. While the majority of the marsquakes occurred in the vicinity of Cerberus Fossae fault system (e.g., Clinton et al., 2021; Drilleau et al., 2022; Durán et al., 2022; Stähler et al., 2022; Zenhäusern et al., 2022), the location of the InSight lander and the marsquake sources close to the crustal dichotomy between the Southern and Northern hemispheres imply that significant lateral variations of relief along the crust-mantle interface and the surface topography could potentially
affect the seismic wave arrival times, especially for multiples, which are seismic phases that are reflected once or
several times at the planet's surface (i.e., PP, PPP, SS, SSS).

In contrast to investigations based purely on seismic data, which are sensitive to the crustal structure beneath the seismic source and the receiver and beneath the bounce points in the case of multiples, gravity and topography data allow one to address the question of lateral crustal variations at the global scale (Wieczorek et al., 2022). Assuming an average value for the crustal thickness of the planet together with densities for the crust and uppermost mantle, this latter approach has the main advantage of being able to estimate crustal thickness everywhere, but suffers from the non-uniqueness of the interpretation of potential fields (Wieczorek et al., 2022). Using the seismological constraints presented in Knapmeyer-Endrun et al. (2021) for a three-layer model beneath the InSight lander as an anchoring point, in combination with gravity and topography data, Wieczorek et al. (2022) estimated an average thickness of the Martian crust lying between 30 and 72 km. Recently, using the same crustal thickness range beneath the InSight lander (Knapmeyer-Endrun et al., 2021), and the distribution of the average crustal thickness along the Rayleigh waves great circle path from the largest marsquake recorded on Mars (Kawamura et al., 2023), Kim et al. (2023) refined the average thickness of the Martian crust to 42–56 km.

In this study, we present a new approach where we use the most recent database of body wave arrival times
combined with gravity field modeling and geodetic data to constrain 3D crustal thickness models of Mars, and 1D
models of the mantle and core. Our results are compared to a case considering a 1D crustal model. Our inversion
scheme fully integrates the thermal history of Mars directly into the forward problem (Drilleau et al., 2021, 2022),
allowing to constrain geodynamical parameters and performing marsquakes' relocation. The strength of our
geodynamical parameterization is that seismic velocities and densities are consistently linked using common
physical assumption, enabling a natural integration of the different geophysical data in the inversion scheme.

## 2. Data and Inversion Methodology

We considered 31 seismic events. The body waves arrival times are detailed in Table S1.1 in Supporting In-
formation S1. We used the database of Drilleau et al. (2022) consisting of 17 low-frequency marsquakes, which
we expanded with additional marsquakes, using the method detailed in Drilleau et al. (2022). In addition to the
direct P- and S-phases, we included phases reflecting below the surface (PP, PPP, SS, SSS) and at the core-mantle
boundary (ScS), and depth phases (pP, sP, sS). Among the seismic events considered, two meteorite impacts
(S1000a, S1094b, Kim et al., 2022; Posiolova et al., 2022) and the largest marsquake recorded on Mars (S1222a,
Kawamura et al., 2023) are included. We also considered core-transiting seismic phase (SKS) recently measured
in Irving et al. (2023) from two farside seismic events detected on Mars (S1000a and S0976a).

To model the 1D interior structure of Mars, we relied on a geodynamical parameterization of the forward problem
in terms of quantities that influence the thermochemical evolution of the planet, accounting for 4.5 Gyr of
planetary evolution (Drilleau et al., 2021, 2022; Samuel et al., 2023). This approach accounts for the heat and
chemical exchanges between several planetary envelopes (a) a convecting liquid iron-rich core; (b) an adiabatic
convecting silicate mantle; and (c) a time-evolving stagnant lithospheric lid that includes a crust enriched in Heat-
Producing Elements (HPE) with respect to the underlying mantle. The thermochemical evolution is computed
following the approach described in detail in Text S1 in Samuel et al. (2019). The list of inverted parameters and
the prior bounds considered, are given in Table S2.1 in Supporting Information S1. While the crustal thickness is
an inversion output (considered here as the global mean crustal thickness of Mars), its complex seismic structure
is decoupled from temperature effects and is directly inverted for using three crustal layers. Based on the con-
straints provided by receiver function estimations and gravity and topography considerations (Knapmeyer-
Endrun et al., 2021; Wieczorek et al., 2022), our maximum permissible average crustal density is set to 3,130 kg/
m³. Following Drilleau et al. (2022), Huang et al. (2022), Irving et al. (2023) we considered a bulk silicate mantle
composition in major and HPE corresponding to that of EH45 (Sanloup et al., 1999). Models based on the bulk
silicate composition of Yoshizaki and McDonough (2020) tested in Irving et al. (2023) showed more difficulty to
fit the SKS data. In the mantle, the computation of density and seismic velocities from thermal profiles are
performed using a thermodynamic model via the Perple_X Gibbs free energy minimization software (Con-
nolly, 2005), using the database of Stixrude and Lithgow-Bertelloni (2011). Following Irving et al. (2018, 2023)
and Samuel et al. (2023), we assume the core to be well mixed, and compute its elastic properties using a third-
order isentropic Birch-Murnaghan equation of state (eos). The eos parameters are estimated from the data together
with all the other model parameters and as such, incorporate the information about the composition of the core.
With this approach, no prior assumptions about the core composition are required and the eos implicitly relates density and acoustic wave velocity in a thermodynamic consistent manner.

For a given 1D interior structure model, a 3D crust of uniform density is generated using the ctplanet software package (Wieczorek, 2021), which uses the observed gravity field and topography as constraints. This approach enables one to determine the relief along the crust-mantle interface, if a density contrast across the crust-mantle interface, a mean crustal density, and a global mean crustal thickness are specified. The values of these quantities are provided by our geodynamical parameterization, which outputs the radial density profile of Mars.

Following previous works (e.g., Drilleau et al., 2020, 2022; Irving et al., 2023; Samuel et al., 2023), given that the origin time of the seismic events remains unknown, the cost function (see Text S.3 in Supporting Information S1) is defined using differential times relative to the P- and S-waves phase arrivals. Because the degree-two tidal Love number ($k_2$) is sensitive to the core radius and mantle rigidity, $k_2$ is computed for each model and included into the cost function, and is compared to the estimate provided in Konopliv et al. (2020). Note that in computing $k_2$ we consider only the elastic, real part of the Love number, which is by far the dominant contribution (Stähler et al., 2021).

The eight main steps of the forward modeling can be summarized as follows. (a) The geodynamical parameters and the marsquake locations are drawn from the prior (Table S2.1 in Supporting Information S1). In the case of 3D models, the back azimuths are randomly sampled within the range specified in Table S1.1 in Supporting Information S1. (b) Following the approach described in Drilleau et al. (2021, 2022) and Samuel et al. (2019), the present-day thermo-chemical state of the planet is used to compute the density, and P- and S- waves radial seismic velocity profiles of Mars assuming compositional, mineralogical, and thermodynamic models for the different planetary envelopes. (c) If the crustal thickness of this 1D model—that we consider to be the average crustal thickness of Mars—does not lie between 30 and 72 km (Wieczorek et al., 2022), the model is rejected. (d) The relief along the crust-mantle interface is determined using the ctplanet software package (Wieczorek, 2021). (e) If the crustal thickness at the InSight landing site is located outside the range estimated in Knapmeyer-Endrun et al. (2021) (39 ± 8 km), the model is rejected. If the model is retained, (f) the body waves arrival times are computed using the TauP toolkit (Crotwell et al., 1999), and a travel time correction is applied for each phase following the method described by Nolet (2008) (Chapter 13.3). The travel time corrections are propagation delays induced by lateral variations of crust thickness relative to the homogeneous crust used in the TauP computations. (g) The $k_2$ Love number is computed following the formalism of Alterman et al. (1959); Takeuchi and Saito (1972). (h) Finally, the cost function is determined. In the case of purely 1D models, the steps (d) and (e) are skipped, and no travel time corrections are applied in step (f).

Due to the ill-posed nature of the problem, the inverse problem is solved considering a Bayesian approach based on a Monte Carlo Markov chain method (e.g., Mosegaard & Tarantola, 1995; Tarantola, 2005) (Text S.3 in Supporting Information S1). To infer the posterior distribution of the parameters, we use the Metropolis algorithm (Metropolis et al., 1953).

### 3. Constraints on the Seismic Structure and Geodynamical Parameters

We considered two sets of inversions: A first set consisting of radial 1D models (referred to as “1D”), and a second set including lateral variations of crustal thickness (referred to as “3D”). The output velocity models are shown in Figure 1. Figure 2 shows the a posteriori distributions of several present-day quantities and initial quantities from the geodynamical forward problem. The mean values of a selection of geodynamical parameters, and their $1-\sigma$ standard deviations, are summarized in Table S2.2 in Supporting Information S1.

The results of the 1D set (Figure 1a) are relatively close to those presented in recent works relying on a similar inversion approach (Drilleau et al., 2022; Irving et al., 2023; Samuel et al., 2023). Indeed, the additional recorded events considered in this study reinforce the previous findings, without necessarily refining them, because (a) the differential time uncertainties of the new seismic events are similar to those used in the previous data set, and (b) half of the new quakes are located near Cerberus Fossae, and their associated body waves follow almost identical paths.

Concerning the models with a 3D crust, an important finding is that the results show a bimodal distribution, with two distinct families retrieved (the a posteriori distributions of $V_P$ and $V_S$ are shown in Text S.4 in Supporting Information S1). Three quarters of the seismic velocity models (Figure 1b) accepted by the Bayesian algorithm
<figure><img src="image_3.png" /><figcaption>Figure 1.</figcaption></figure>
are relatively similar to the 1D models output (Figure 1a), and the associated geodynamical parameters show relatively close values (Table S2.2 in Supporting Information S1). In the following, this subset is referred to as “Family A (3D).” The remaining quarter of the 3D models shows a clearly different seismic structure (Figure 1c). We refer to this second subset as “Family B (3D).”

Figures 1d–1f, display the temperature, $V_p$ and $V_s$ in the mantle for three representative models from these families. Compared to the 1D case and to Family A, the Family B case is characterized by a present-day mantle that is on average 300 K colder (Figures 1d, 2c1–2c3, and 2d1–2d3) and that has an about 50% thinner thermal lithosphere at present-day (Figure 1d; Table S2.2 in Supporting Information S1). This results from the fact that for Family B, the initial thermal state is hotter (Figures 2f1–2f3), which reduces the lithospheric growth compared to the 1D case and Family A. The thermal gradient in the lithosphere (corresponding to the region between the base of the crust and the top of the uppermost thermal boundary layer) for Family B is generally 15%–20% smaller than that for Family A and the 1D case (Table S2.2 in Supporting Information S1; Figures 2b1–2b3). Despite the thinner lithosphere, this is mainly due to the fact that the mantle is colder and has a smaller effective activation energy, $E^*$ (Table S2.2 in Supporting Information S1). Indeed, the lithospheric gradient $dT_l/dr$ shows a dependence with the uppermost mantle temperature $T_m$ and mantle rheology as ~$T_m(1 - 20 T_m/E^*)$ (see Text S.2 in Supporting Information S1), and the combination of smaller $T_m$ and $E^*$ for Family B results in smaller values of the lithospheric gradient (1.9 K/km vs. 2.0 K/km for Family A and the 1D case, see Table S2.2 in Supporting Information S1; Figure 2b1–b3). This leads to distinct temperature profiles that result in different mineralogical phase proportions as a function of depth (Figures 1g–1i). Due to its exothermic nature, the olivine-to-wadsleyite phase transition is located ~40–70 km shallower for the cold temperature profile of Family B, compared to the 1D case and to Family A. The selected model for Family B also shows a clear seismic discontinuity at ~750 km depth, due to the phase change from orthopyroxene to high-pressure clinopyroxene. For the selected models of the 1D case and Family A, this phase change occurs almost at the same depth as the olivine-to-wadsleyite phase transition, therefore their associated seismic discontinuities appear combined near 1,030–1,050 km depth. Using observed triplicated P and S waves from five teleseismic events, Huang et al. (2022) constrained the depth of the post olivine transition to be 1,006 ± 40 km. Most of our models have their olivine-to-wadsleyite phase transition occurring within this uncertainty range.

Other noticeable differences appear for Family B, compared to the 1D case and Family A. Family B shows a mantle initially 90 K hotter and with a reference mantle viscosity 1–2 orders of magnitude smaller (Table S2.2 in Supporting Information S1), leading to a faster mantle cooling compared to the other models. In addition, Family B models have half the activation energy, but a slightly larger effective activation volume. The crustal enrichment, defined as the ratio of heat production in the crust to the total heat production in the bulk silicate envelope crust plus mantle, is 20% smaller than that for the 1D and Family A cases (Figures 2e1–2e3), indicating that the crust is less enriched in HPE with respect to the primitive mantle. For all models, the surface heat flux lies within similar ranges: 20–23 mW m⁻² (Figures 2a1–2a3). This results from the fact that the same bulk content in HPE is used for all models, which eventually controls the present-day heat loss at the surface (Drilleau et al., 2021). Quantities associated with the core are also different, with a core radius 100 km smaller and thus 6% denser than the other sets of models due to constraints on the mass of Mars (Table S2.2 in Supporting Information S1). However, $V_p$ in the core are relatively similar between the different sets of models, because of the two SKS phases (Irving et al., 2023) we use in our data set, which constrain the velocities in the core.

The gap in the 3D output models distribution leading to the two distinct families A and B is due to combinations of geodynamical parameters (mainly the mantle rheology, the crustal enrichment, and to a lesser extent the initial thermal state) that result in present-day crustal thicknesses larger than the maximum bound of 72 km, and/or extreme HPE depletion of the mantle (see Text S.2.3 in Supporting Information S1).

**Figure 1.** Inversion results. 1D models are shown in red, and 3D models are shown in blue and green for Families A and B, respectively. In (a–c) are displayed the output $V_p$ and $V_s$ models in color. The corresponding pdfs are represented in Figure S4.2 in Supporting Information S1. The gray areas indicate the prior bounds. (d) Shows the temperature profile in the mantle for three representative selected models, and (e, f) show the corresponding $V_p$ and $V_s$ selected models. (g–i) display the phase proportion of the minerals for the selected models in (d–f). Abbreviations for mantle minerals are as follows: olivine (ol), wadsleyite (wad), ringwoodite (ring), ferropericlase (fp), clinopyroxene (cpx), orthopyroxene (opx), high-pressure clinopyroxene (C2/c), garnet (gt), feldspar (fs). (j–l) show, for a selection of models, the core radius versus $\tau_{FCN}$, the core radius versus F, and the correlation between $\tau_{FCN}$ and F. Gray shaded areas represent 1$\sigma$, 2$\sigma$, and 3$\sigma$ uncertainties from Le Maistre et al. (2023). The colored diamonds show the three models of each family selected in (d–f) and (g–i).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Posterior distributions of models considering the 1D case (left, in red), and 3D cases for Families A (middle, in blue) and B (right, in green). Present-day quantities are shown for (a) the average surface heat flow, (b) the average temperature gradient in the lithosphere (excluding the crust), (c) the mantle potential temperature, and (d) the temperature profile. The crustal enrichment factor is shown in (e) and the initial uppermost mantle temperature is shown in (f).</figcaption></figure>
4. Datafit and Compatibility With Mars Nutations

The fit to the seismic data is displayed in Figure S3.1 in Supporting Information S1. The three families of models are able to reproduce the observed differential arrival times equally well. We thus attribute the emergence of Family B to the multiple trade-offs between the velocity model, the location of the seismic events, and the travel time corrections due to the lateral crustal thickness variations. 33% of Family B shows a $k_2$ value smaller than the 2σ uncertainty in Konopliv et al. (2020), because of their stiffer mantles due to their colder temperature profiles (Figure S3.2 in Supporting Information S1). In Text S.5 in Supporting Information S1, we present additional inversion sets using the body wave differential times and $k_2$ separately.

Additionally, we tested the compatibility of our models (using a representative subset) with the nutations of Mars recently measured by the Rotation and Interior Structure Experiment (RISE) on InSight (Le Maistre et al., 2023) (Figures 1j–1l). Nutations can be characterized by the Free Core Nutation (FCN) period $\tau_{FCN}$, a rotational normal mode related to the misalignment between the core and mantle rotation axes, and the core amplification factor $F$, which is related to the Moment of Inertia (MoI) of the core (e.g., Dehant & Mathews, 2015). Both $\tau_{FCN}$ and $F$ depend on the MoI of the core, and on the deformation of the mantle resulting from the deformation of the core due to tidal forcing and rotational rate variation. $\tau_{FCN}$ is also proportional to the dynamic shape of the core, which mainly depends on its rotation rate, on the density jump between the core and mantle, and to a lesser extent on mass anomalies within the mantle and at the surface (Le Maistre et al., 2023).

The geometric and dynamic shape of Mars can be explained with two mass-sheet anomalies, one placed deep within the lithosphere and the other resulting from the surface topography (Le Maistre et al., 2023). Both mass-sheet anomalies significantly affect the shape of the core. By varying the depth of the internal load, the effect on the core shape can be modulated. Le Maistre et al. (2023) showed that for the case of the two mass-sheet assumption, a depth of more than 500 km is required to explain the observed FCN period. Since the convective part of the mantle behaves as a fluid over long time scales, the maximal depth of the internal load is limited by the thickness of the lithosphere.

We use representative subsets of 1D, Family A, and B models to assess how RISE data could improve our knowledge about the interior structure of Mars. The subset of 1D models are selected among models that agree with the MoI estimates of Konopliv et al. (2020) at 2σ.

All models except Family B models with a core radius smaller than about 1,750 km agree with $F$ at $2\sigma$ (Figure 1k). For a given core radius, $F$ is larger for 1D models since the latter have a denser core (see Table S2.2 in Supporting Information S1) and a larger core MoI compared to Family A and B models. A smaller fraction of all models (1D: 40%, Family A: 90%, Family B: 21%) agree simultaneously with $F$ and the FCN period (Figure 1l) at $2\sigma$, even if the depth of the internal load of all models is more than about 500 km. All these models have an internal load placed at the bottom of their lithosphere if their crustal density lies on the lower end of the range obtained in the joint inversion, or up to 300 km higher if the crust density is on the high side. This requirement is fulfilled for about 90% of the Family A models, but for both 1D and Family B, a large fraction of models have a lithosphere thickness-crust density combination that results in a FCN period above that measured by RISE.

Among the models considered here, Family B models have the smallest density jump at the CMB (i.e., those corresponding to coldest mantle, highest mantle density, and lowest core density). For this reason, the effect of rotation on the shape of the core and FCN period is the largest (Le Maistre et al., 2023). The internal load affects the FCN by a less than 10 days, an effect that is too small to bring most of Family B models in line with RISE data, but large enough for a large subset of Family A models to agree with RISE. Unlike Family A and B, for 1D models, the rotational flattening alone results in an FCN period in reasonably a good agreement with RISE. Consequently, the internal load has to compensate for the effect of the surface load on the core shape. Although 1D models have a thick lithosphere, the latter is too thin for 60% of models to agree with RISE because of their low crust density.

5. Crustal Thickness and Location of Seismic Events

In the absence of constraint on the MoI in our inversions, the crustal density of 1D models is unconstrained (Figure S6.1d in Supporting Information S1). 3D models provide similar results for the crustal density, with mean fitting values of 3,014 and 2,992 kg m⁻³ for Families A and B, respectively (see Table S2.2, Figures S6.1e, and S6.1f in Supporting Information S1). Compared to Wieczorek et al. (2022), our obtained crustal
densities lie in the upper range of their estimations, and the implied density contrasts between the crust and mantle lie in their lower range (Figures S6.1h and S6.1i in Supporting Information S1). The low density contrasts that we find (447 and 466 kg m⁻³ for Families A and B) have the effect of increasing the amplitude of the relief along the crust-mantle interface in the global crustal thickness models when compared to smaller values. Our crustal densities are comparable to those expected from the average upper crustal composition of Taylor and McLennan (2009). They are, however, somewhat lower than what would be expected based on orbital gamma-ray spectroscopy measurements of surface materials that are basaltic in nature with high FeO contents (3,200–3,450 kg m⁻³, see Baratoux et al., 2014). Our densities are also in the upper range of those expected from in situ analyses of felsic rocks on the surface (from 2,680 to 3,020 kg m⁻³, see Foley et al., 2003; Brückner et al., 2003; Sautter et al., 2015). Consistent with the interpretation of Wieczorek et al. (2022), our bulk crustal density suggests the presence of substantial amounts of materials in the crust that are more felsic than the surface materials.

For the three families of models, the retrieved average crustal thickness are in the upper range of Wieczorek et al. (2022), with 68.7, 64.2, and 63.3 km for 1D, Family A and Family B models, respectively (Table S2.2 in Supporting Information S1). The crustal thickness is largely influenced by the bulk HPE content (Drilleau et al., 2021; Knapmeyer-Endrun et al., 2021), the latter being associated with the planet composition (in our case the EH45). HPE favor crustal production, resulting in crustal thickness values larger than 50 km (Figures S6.1a–S6.1c in Supporting Information S1). However, to satisfy receiver function constraints on the crustal thickness beneath InSight, models associated with a thick crust (>72 km) are rejected. This can be seen by comparing the prior and posterior distributions for the crustal thickness (Figures S6.1a–S6.1c in Supporting Information S1). An appreciable fraction of models is rejected because of this criteria. Though this lies in the range of 30–72 km of Wieczorek et al. (2022), our thickness is somewhat larger than the estimation of the 44–58 km interquartile range from Kim et al. (2023) that was deduced from multi-orbiting surface waves along great circle paths associated with the S1222a marsquake. Nevertheless, we note that the average crustal thickness range along this great circle path for our Family A and Family B models show an interquartile range of 60–67 and 60–68 km, respectively, overlapping with the upper ~25% of models in the posterior distribution of Kim et al. (2023). Because the InSight landing site is located at low elevations in the Northern hemisphere of Mars, the crustal thickness at this location is expected to be thinner than the global average crustal thickness. Our models predict a total thickness of about 44.9 and 44.6 km (Families A and B) at the InSight landing site (Table S2.2 in Supporting Information S1). These values are consistent with the uppermost values obtained by receiver function analyses, namely 39 ± 8 km (Knapmeyer-Endrun et al., 2021) and 37–47 km (Durán et al., 2022).

The global crustal thickness model associated with the best cost function for Family A is shown in map view in Figure 3. The mean locations of the seismic events for the three different models are shown (see also Table S7.1 and Figure S7.1 in Supporting Information S1). For the 31 seismic events used in this analysis, the mean deviation compared to the locations retrieved with 1D models are respectively for Families A and B, 0.5 and 0.9° for epicentral distance, and 2.8 and 4.8° for back azimuth. It corresponds to mean deviations of ~32 km and ~56 km for Families A and B, respectively (Table S7.2 in Supporting Information S1). The locations associated with Family A are closest to the 1D case given that their seismic profiles are relatively similar (Figure 1). The seismicity distribution is more spread for Family B, in particular at Cerberus Fossae (Figure 3c). We redid the same inversions by removing the marsquakes located at Cerberus Fossae, where the crustal thickness is relatively homogeneous. Similar velocity distributions were obtained.

The predicted crustal thicknesses at the 31 epicenter locations are displayed in Figure S6.2 in Supporting Information S1. For the events located near Cerberus Fossae, the crustal thickness lies between 40 and 60 km. When the locations are closer to the crustal dichotomy (e.g., S0325a, S1222a) or in the vicinity of major geological structures, such as Utopia Planitia (S0185a and S0861a), Olympus Mons (S1102a) or Valles Marineris (S0976a), the crustal thickness distributions are clearly multimodal. The differential travel times corrections for the 3D models (Figures S8.1 and S8.2 in Supporting Information S1) are directly linked to the crustal thickness values at the epicenter locations. If the absolute travel time correction values are smaller than 1.5 s for the S-P differential time, they can reach 5 and 8 s for PP-P and SS-S, and up to 8 and 13 s for PPP-P and SSS-S, respectively.

Most of the events described in this study have been interpreted as originating from possible tectonic activities in Drilleau et al. (2022). Excluding both impact events (S1000a and S1094b), 12 new events provide additional
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Best-fitting global map of Martian crustal thickness and event locations. For each event, red, blue and green dots are the inferred locations from the 1D, Family A (3D) and Family B (3D) models, respectively (see Table S7.1 in Supporting Information S1). For better readability, white and black lines connect the same event. Known impact locations detected from satellite images (Posiolova et al., 2022) are indicated by black stars. Close-up views in the Cerberus Fossae region are shown in (b, c) and present respectively the 1D-3D (Family A) and 1D-3D (Family B) solutions. The background maps represent the crustal thickness (best model of Family A from this study) an overlying shaded-relief map based on laser altimetry (Smith et al., 2001). Red and black lines are global compilation of normal and reverse faults, respectively (Knapmeyer et al., 2006). The thick red lines highlight the Cerberus Fossae grahen system (Perrin et al., 2022).</figcaption></figure>
information to the distribution of seismicity on Mars (Figure 3). Even when considering potential large uncertainties (especially in back azimuth determination, (see Table S7.1 and Figure S7.1 in Supporting Information S1), six of these events (S1012d, S1015f, S1022a, S1048d, S1133c, and S1157a) are located in the vicinity of large fossae located East of InSight. This confirms that Cerberus Fossae (and possibly Grjótá Valles and Elysium Fossae) are major tectonic structures (e.g., Perrin et al., 2022) that generated most of the seismic activity in the vicinity of the InSight lander (e.g., Giardini et al., 2020; Jacob et al., 2022; Stähler et al., 2022; Zenhäusern et al., 2022). Three new events (S1039b, S1222a and S1197a) are located in the Southern highlands or on the dichotomy boundary. The latter is characterized by large North to South variations in elevation and crustal thickness that could give rise to tectonic stresses in the lithosphere, potentially accounting for these events (e.g., Watters, 2003a,b). Finally, three distant seismic events (S0976a, S1102a and S1153a) are located close to the Tharsis Rise. The large uncertainties in back azimuth encompass major structures West of Tharsis, mostly extensional, from large fossae in the South (Sirenum, Memnonia, Mangala), to the Northern flanks of Olympus Mons (Lycus Sulci) in the North, through Gordii Dorsum. East of Tharsis, the numerous faults surrounding Valles Marineris might also trigger seismic events, like S0976a.

## 6. Conclusions

In this study, we improved our probabilistic inversion method that was initially developed to constrain the 1D interior structure of Mars (Drilleau et al., 2022; Irving et al., 2023; Samuel et al., 2023; Stähler et al., 2021) by accounting for lateral variations in Martian crustal thickness inferred from gravity data. When the effect of a 3D crust is considered on body wave travel times, our results showed that two families of models are possible, leading to different interpretations of Mars' evolution: (a) a hot case scenario (Family A), representing the majority of models (~3/4) with seismic velocity profiles and geodynamical parameters almost identical to the 1D case, and (b) a cold case scenario (Family B), characterized by a present-day mantle that is about 300 K colder. Although both families provide a similar fit to the seismic data, a smaller fraction of Family B models agree with RISE data (Le Maistre et al., 2023) and $k_2$. We found that the crustal thickness of Mars is ~66 km on average, and that the thickness beneath the InSight lander is about 45 km, both of which are on the upper end of the previous estimates of Knapmeyer-Endrun et al. (2021); Wieczorek et al. (2022). Our predictions of the crustal density are compatible with estimates of the average upper crustal composition from Martian soils (Taylor & McLennan, 2009), but are lower than what is expected for typical basaltic materials on Mars (e.g., Baratoux et al., 2014) and somewhat higher than what is expected for felsic surface rocks analyzed in situ (e.g., Sautter et al., 2015). We updated the seismicity map of Mars of Drilleau et al. (2022) by adding 12 new events. Even if the presence of heterogeneities in the crust and mantle (e.g., Samuel et al., 2023; Smrekar et al., 2019), and the shape of the core (Le Maistre et al., 2023), are inevitably sources of error on our results, this study is a first step to go beyond the assumption of a spherically symmetrical planet. Further inversion sets considering the recently discovered basal mantle layer above Mars' core (Samuel et al., 2023) will investigate the potentially substantial influence of deep mantle stratification on our findings.

## Data Availability Statement

All raw waveform data is available through IPGP data center (InSight Mars SEIS Data Service, 2019b), IRIS-DMC (InSight Marsquake Service, 2023), and NASA-PDS (InSight Mars SEIS Data Service, 2019a). The hundred best models of each family are available in Drilleau (2024).

## References

* Alterman, Z., Jarosh, H., & Perekis, C. L. (1959). Oscillations of the Earth. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 252(1268), 80-95.

* 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*, 1-14(3), 183-189. https://doi.org/10.1038/s41561-020-0544-y

* Baratoux, D., Samuel, H., Michaut, C., Toplis, M. J., Monnereau, M., Wieczorek, M., et al. (2014). Petrological constraints on the density of the Martian crust. *Journal of Geophysical Research: Planets*, 119(7), 1707-1727. https://doi.org/10.1002/2014JE004642

* Brückner, J., Dreibus, G., Rieder, R., & Wänke, H. (2003). Refined data of alpha proton x-ray spectrometer analyses of soils and rocks at the Mars Pathfinder site: Implications for surface chemistry. *Journal of Geophysical Research*, 108(E12), 8094. https://doi.org/10.1029/2003JE002060

* 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. A. D. (2005). Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation. *Earth and Planetary Science Letters*, 236(1-2), 524-541. https://doi.org/10.1016/j.epsl.2005.04.033

## Acknowledgments

We acknowledge NASA, CNES, partner agencies and institutions (UKSA, SSO, DLR, JPL, IPGP-CNRS, ETHZ, ICL, and MPS-MPG), and the operators of JPL, SISMOC, MSDS, IRIS-DMC, and PDS for providing SEED SEIS data. M.D. was granted access to the GENCI HPC resources of IDRIS under allocation AD010413017R1. Numerical computations were partly performed on the S-CAPAD/DANTE platform, IPGP, France. Co-authors affiliated to French laboratories thank the French Space Agency CNES and ANR fund (ANR-19-CE31-0008-08). A.R. was supported by the Belgian PRODEX program managed by the European Space Agency in collaboration with the Belgian Federal Science Policy Office, contract number PEA4000140326. This is InSight Contribution Number 329. The authors acknowledge Angela Marusiak and an anonymous reviewer for their constructive comments that helped to improve this paper.
Crotwell, H. P., Owens, T. J., & Ritsema, J. (1999). The TauP toolkit: Flexible seismic travel-time and ray-path utilities. *Seismological Research Letters*, 70(2), 154–160. https://doi.org/10.1785/gssrl.70.2.154

Dehant, V., & Mathews, P. M. (2015). *Precession, nutation and wobble of the Earth*. Cambridge University Press.

Drilleau, M. (2024). Constraints on lateral variations of Martian crustal thickness from seismological and gravity field measurements [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.10517921

Drilleau, M., Beucler, E., Lognonné, P., Panning, M., Knapmeyer-Endrun, B., Banerdt, B. W., et al. (2020). MSS/1: Single-station and single-event marsquake inversion. *Earth and Space Science*, 7(12). https://doi.org/10.1029/2020EA001118

Drilleau, M., Samuel, H., Garcia, R. F., Rivoldini, A., Perrin, C., Michaut, C., et al. (2022). Marsquake locations and 1-D seismic models for Mars from InSight data. *Journal of Geophysical Research: Planets*, 127(9), e2021JE007067. https://doi.org/10.1029/2021JE007067

Drilleau, M., Samuel, H., Rivoldini, A., Panning, M., & Lognonné, P. (2021). Bayesian inversion of the Martian structure using geodynamic constraints. *Geophysical Journal International*, 226(3), 1615–1644. https://doi.org/10.1093/gji/ggab105

Durán, C., Khan, A., Ceylan, S., Zenhäusern, G., Stähler, S., Clinton, J., & Giardini, D. (2022). Seismology on Mars: An analysis of direct, reflected, and converted seismic body waves with implications for interior structure. *Physics of the Earth and Planetary Interiors*, 325, 106851. https://doi.org/10.1016/j.pepi.2022.106851

Foley, C. N., Economou, T., & Clayton, R. N. (2003). Final chemical results from the Mars Pathfinder alpha proton X-ray spectrometer. *Journal of Geophysical Research*, 108(E12), 8096. https://doi.org/10.1029/2002JE002019

Garcia, R. F., Daubar, I. J., Beucler, É., Posiolova, L. V., Collins, G. S., Lognonné, P., et al. (2022). Newly formed craters on Mars located using seismic and acoustic wave data from InSight. *Nature Geoscience*, 15(10), 774–780. https://doi.org/10.1038/s41561-022-01014-0

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. https://doi.org/10.1038/s41561-020-0539-8

Huang, Q., Schmerr, N. C., King, S. D., Kim, D., Rivoldini, A., Plesa, A.-C., et al. (2022). Seismic detection of a deep mantle discontinuity within Mars by InSight. *Proceedings of the National Academy of Sciences of the United States of America*, 119(42). https://doi.org/10.1073/pnas.2204474119

InSight Marsquake Service. (2023). Mars seismic catalogue, InSight mission; v13 2023-01-01 [Dataset]. ETHZ, IPGP, JPL, ICL, University of Bristol. https://doi.org/10.12686/a19

InSight Mars SEIS Data Service. (2019a). InSight SEIS data bundle. PDS geosciences (GEO) node [Dataset]. InSight Mars SEIS Data Service. https://doi.org/10.17189/1517570

InSight Mars SEIS Data Service. (2019b). SEIS raw data, InSight mission [Dataset]. IPGP, JPL, CNES, ETHZ, ICL, MPS, ISAE-Supaero, LPG, MFSC. https://doi.org/10.18715/SEIS.INSIGHT.XB_2016

Irving, J. C. E., Cottaar, S., & Lekić, V. (2018). Seismically determined elastic parameters for Earth's outer core. *Science Advances*, 4(6), eaar2538. https://doi.org/10.1126/sciadv.aar2538

Irving, J. C. E., Lekić, V., Duran, C., Drilleau, M., Kim, D., Rivoldini, A., et al. (2023). First observations of core-transiting seismic phases on Mars. *Proceedings of the National Academy of Sciences of the United States of America*, 120(18), e2217090120. https://doi.org/10.1073/pnas.221709012

Jacob, A., Plasman, M., Perrin, C., Fuji, N., Lognonné, P., Xu, Z., et al. (2022). Seismic sources of insight marsquakes and seismotectonic context of Elysium Planitia, Mars. *Tectonophysics*, 837, 229434. https://doi.org/10.1016/j.tecto.2022.0229434

Kawamura, T., Clinton, J. F., Zenhäusern, G., Ceylan, S., Horleston, A. C., Dahmen, N. L., et al. (2023). S1222a—The largest marsquake detected by InSight. *Geophysical Research Letters*, 50(5), e2022GL015143. https://doi.org/10.1029/2022GL015143

Khan, A., Ceylan, S., van Driel, M., Giardini, D., Lognonné, P., Samuel, H., et al. (2021). Upper mantle structure of Mars from InSight seismic data. *Science*, 373(6553), 434–438. https://doi.org/10.1126/science.abf2966

Kim, D., Banerdt, W. B., Ceylan, S., Giardini, D., Lekić, V., Lognonné, P., et al. (2022). Surface waves and crustal structure on Mars. *Science*, 378(6618), 417–421. https://doi.org/10.1126/science.abq7157

Kim, D., Duran, C., Giardini, D., Plesa, A.-C., Stähler, S. C., Boehm, C., et al. (2023). Global crustal thickness revealed by surface waves orbiting Mars. *Geophysical Research Letters*, 50(12), e2023GL013482. https://doi.org/10.1029/2023GL013482

Knapmeyer, M., Oberst, J., Hauber, E., Wählisch, M., Deuchler, C., & Wagner, R. (2006). Working models for spatial distribution and level of Mars' seismicity. *J. Geophys. Res. Planets*, 111(11), 1–23. https://doi.org/10.1029/2006JE002708

Knapmeyer-Endrun, B., Panning, M. P., Bissig, F., Joshi, R., Khan, A., Kim, D., et al. (2021). Thickness and structure of the Martian crust from InSight seismic data. *Science*, 373(6553), 438–443. https://doi.org/10.1126/science.abf8966

Konopliv, A. S., Park, R. S., Rivoldini, A., Baland, R. M., Le Maistre, S., Van Hoolst, T., et al. (2020). Detection of the chandler wobble of Mars from orbiting spacecraft. *Geophysical Research Letters*, 47(21), 1–9. https://doi.org/10.1029/2020GL090568

Le Maistre, S., Rivoldini, A., Caldiero, A., Yseboodt, M., Baland, R.-M., Beuthe, M., et al. (2023). Spin state and deep interior structure of Mars from InSight radio tracking. *Nature*, 619(7971), 733–737. https://doi.org/10.1038/s41586-023-06150-0

Li, J., Beghein, C., Wooley, J., Davis, P., Lognonné, P., Schimmel, M., et al. (2022). Evidence for crustal seismic anisotropy at the InSight lander site. *Earth and Planetary Science Letters*, 593, 117654. https://doi.org/10.1016/j.epsl.2022.117654

Lognonné, P., Banerdt, W. B., Pike, W. T., Giardini, D., Christensen, U., Garcia, R. F., et al. (2020). Constraints on the shallow elastic and anelastic structure of Mars from InSight seismic data. *Nature Geoscience*, 13(3), 213–220. https://doi.org/10.1038/s41561-020-0536-y

Lognonné, P., Banerdt, W. B., Giardini, D., Pike, W. T., Christensen, U., Laudet, P., et al. (2019). Seis: Insight's seismic experiment for internal structure of Mars. *Space Science Reviews*, 215(1), 12. https://doi.org/10.1007/s11214-018-0574-6

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. (1953). Equation of state calculations by fast computing machines. *Journal of Chemical Physics*, 21(6), 1087–1091. https://doi.org/10.1063/1.1699114

Mosegaard, K., & Tarantola, A. (1995). Monte-Carlo sampling of solutions to inverse problems. *Journal of Geophysical Research*, 100(B7), 12431–12447. https://doi.org/10.1029/94jb03097

Nolet, G. (2008). *A breviation of seismic tomography*. Cambridge University Press.

Perrin, C., Jacob, A., Lucas, A., Myhill, R., Hauber, E., Batov, A., et al. (2022). Geometry and segmentation of Cerberus fossae: Mars: Implementations for marsquake properties. *Journal of Geophysical Research: Planets*, 127(1), e2021JE007118. https://doi.org/10.1029/2021JE007118

Posiolova, L. V., Lognonné, P., Banerdt, W. B., Clinton, J., Collins, G. S., Kawamura, T., et al. (2022). Largest recent impact craters on Mars: Orbital imaging and surface seismic co-investigation. *Science*, 378(6618), 412–417. https://doi.org/10.1126/science.abq7704

Samuel, H., Drilleau, M., Rivoldini, A., Xu, Z., Huang, Q., Garcia, R. F., et al. (2023). Geophysical evidence for an enriched molten silicate layer above Mars's core. *Nature*, 622(7984), 712–717. https://doi.org/10.1038/s41586-023-06601-8

Samuel, H., Lognonné, P. H., Panning, M., & Lainey, V. (2019). The rheology and thermal history of Mars revealed by the orbital evolution of Phobos. *Nature*, 569(7757), 523–527. https://doi.org/10.1038/s41586-019-1202-7
Sanloup, C., Jambon, A., & Gillet, P. (1999). A simple chondritic model of Mars. *Physics of the Earth and Planetary Interiors*, 112(1-2), 43-54.
https://doi.org/10.1016/s0031-9201(98)00175-7

Sautter, V., Toplis, M. J., Wiens, R. C., Cousin, A., Fabre, C., Gasnault, O., et al. (2015). In situ evidence for continental crust on early Mars. *Nature Geoscience*, 8, 605-609. https://doi.org/10.1038/NGE02474

Smith, D. E., Zuber, M. T., Frey, H. V., Garvin, J. B., Head, J. W., Muhleman, D. O., et al. (2001). Mars orbiter laser altimeter: Experiment summary after the first year of global mapping of Mars. *Journal of Geophysical Research*, 106(E10), 23689-23722. https://doi.org/10.1029/2000je001364

Smrekar, S. E., Lognonné, P., Spohn, T., Banerdt, B., Breuer, D., Christensen, U., et al. (2019). Pre-mission insights on the interior of Mars. *Space Science Reviews*, 215(1), 3. https://doi.org/10.1007/s11214-018-0563-9

Stähler, S. C., Khan, A., Banerdt, W. B., Lognonné, P., Giardini, D., Ceylan, S., et al. (2021). Seismic detection of the Martian core. *Science*, 373(6553), 443-448. https://doi.org/10.1126/science.abi7730

Stähler, S. C., Mittelholz, A., Perrin, C., Kawamura, T., Kim, D., Knapmeyer, M., et al. (2022). Tectonics of Cerberus Fossae unveiled by marsquakes. *Nature Astronomy*, 6(12), 1376-1386. https://doi.org/10.1038/s41550-022-01803-y

Stixrude, L., & Lithgow-Bertelloni, C. (2011). Thermodynamics of mantle minerals - II. Phase equilibria. *Geophysical Journal International*, 184(3), 1180-1213. https://doi.org/10.1111/j.1365-246X.2010.04890.x

Takeuchi, H., & Saito, M. (1972). Seismic surfaces waves. In B. A. Bolt (Ed.), *Methods in computational physics* (Vol. 1, pp. 217-294). Academic Press.

Tarantola, A. (2005). *Inverse problem theory—And methods for model parameter estimation*. SIAM.

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

Watters, T. R. (2003a). Lithospheric flexure and the origin of the dichotomy boundary on Mars. *Geology*, 31(3), 271-274. https://doi.org/10.1130/0091-7613(2003)031(0271:LFATOO)2.0.CO;2

Watters, T. R. (2003b). Thrust faults along the dichotomy boundary in the Eastern hemisphere of Mars. *Journal of Geophysical Research*, 108(E6), 5054. https://doi.org/10.1029/2002JE001934

Wieczorek, M. A. (2021). *Ctplanet, version 0.2.1*. Zenodo. https://doi.org/10.5281/zenodo.4439426

Wieczorek, M. A., Broquet, A., McLennan, S. M., Rivoldini, A., Golombek, M., Antonangeli, D., et al. (2022). InSight constraints on the global character of the Martian crust. *Journal of Geophysical Research: Planets*, 127(5), e2022JE007298. https://doi.org/10.1029/2022JE007298

Yoshizaki, T., & McDonough, W. F. (2020). The composition of Mars. *Geochimica et Cosmochimica Acta*, 273, 137-162. https://doi.org/10.1016/j.gca.2020.01.011

Zenhäusern, G., Stähler, S. C., Clinton, J. F., Giardini, D., Ceylan, S., & Garcia, R. F. (2022). Low-frequency marsquakes and where to find them: Back azimuth determination using a polarization analysis approach. *Bulletin of the Seismological Society of America*, 112(4), 1787-1805.
https://doi.org/10.1785/0120220019

## References From the Supporting Information

Breuer, D., & Spohn, T. (2003). Early plate tectonics versus single-plate tectonics on Mars: Evidence from magnetic field history and crust evolution. *Journal of Geophysical Research*, 108(E7), 5072. https://doi.org/10.1029/2002JE001999

Grott, M., & Breuer, D. (2008). The evolution of the Martian elastic lithosphere and implications for crustal and mantle rheology. *Icarus*, 193(2), 503-515. https://doi.org/10.1016/j.icarus.2007.08.015

Hauck, S. A., & Phillips, R. J. (2002). Thermal and crustal evolution of Mars. *Journal of Geophysical Research*, 107(E7), 6-1-6-19. https://doi.org/10.1029/2001JE001801

Morschhauser, A., Grott, M., & Breuer, D. (2011). Crustal recycling, mantle dehydration, and the thermal evolution of Mars. *Icarus*, 212(2), 541-558. https://doi.org/10.1016/j.icarus.2010.12.028

Thiriet, M., Breuer, D., Michaut, C., & Plesa, A.-C. (2018). Scaling laws of convection for cooling planets in a stagnant lid regime. *Physics of the Earth and Planetary Interiors*, 286, 138-153. https://doi.org/10.1016/j.pepi.2018.11.003

Tosi, N., Godolt, M., Stracke, B., Ruedas, T., Grenfell, J. L., Höning, D., et al. (2017). The habitability of a stagnant-lid Earth. *Astronomy and Astrophysics*, 605, A71. https://doi.org/10.1051/0004-6361/201730728