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

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.


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(Drilleau et al., , 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.

Data and Inversion Methodology
We considered 31 seismic events.The body waves arrival times are detailed in Table S1.1 in Supporting Information 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(Drilleau et al., , 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 constraints 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 3 .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 (Connolly, 2005), using the database of Stixrude and Lithgow-Bertelloni (2011).Following Irving et al. (2018Irving et al. ( , 2023) ) and Samuel et al. (2023), we assume the core to be well mixed, and compute its elastic properties using a thirdorder 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.

10.1029/2023GL105701
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., 2020Drilleau et al., , 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. (2021Drilleau et al. ( , 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).

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-σ 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 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, 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 2 (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 τ FCN , the core radius versus F, and the correlation between τ FCN and F. Gray shaded areas represent 1σ, 2σ, and 3σ uncertainties from Le Maistre et al. ( 2023).The colored diamonds show the three models of each family selected in (d-f) and (g-i).

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 τ 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 τ 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.τ 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 masssheet 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σ (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σ, 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.

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 3 for Families A and B, respectively (see Table S2 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 3 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 3 , 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 3 , 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  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  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   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 graben system (Perrin et al., 2022).

Geophysical Research Letters
10.1029/2023GL105701 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.

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.

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

Figure 3 .
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 TableS7.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 graben system(Perrin et al., 2022).