Focal Mechanism Determination of Event S1222a and Implications for Tectonics Near the Dichotomy Boundary in Southern Elysium Planitia, Mars

On 4 May 2022 the InSight seismometer SEIS‐VBB recorded the largest marsquake ever observed, S1222a, with an initial magnitude estimate of MWMa ${M}_{W}^{\mathrm{M}\mathrm{a}}$ 4.6. Understanding the depth and source properties of this event has important implications for the nature of tectonic activity on Mars. Located ∼37° to the southeast of InSight, S1222a is one of the few non‐impact marsquakes that exhibits prominent surface waves. We use waveform modeling of body waves (P and S) and surface waves (Rayleigh and Love) to constrain the focal mechanism, assuming a double‐couple source, and find that S1222a likely resulted from reverse faulting in the crust (source depth near 22 km). We estimate the scalar moment to be 2.5 × 1015–3.5 × 1015 Nm (magnitude MW 4.2–4.3). Our results suggest active compressional tectonics near the dichotomy boundary on Mars, likely due to thermal contraction from planetary cooling.


10.1029/2023JE007793
2 of 13 due to unknown three-dimensional (3D) structure; and (d) the moment tensor solution requires knowledge of the source location, which may have large uncertainty.However, since different seismic phases have different radiation patterns and sample different portions of the focal sphere, single station moment tensor inversions are possible for sufficiently high-quality events with known locations, provided that sufficiently accurate structural models are available.
After more than two Martian years on the surface, the InSight (Interior Exploration using Seismic Investigations, Geodesy, and Heat Transport, Banerdt et al., 2020) lander and associated seismometer (Lognonné et al., 2019) have helped transform our understanding of the interior and dynamics of Mars.InSight landed in the Elysium Planitia of Mars on 26 November 2018 and has recorded over 1,000 marsquakes (Ceylan et al., 2022) with distinct spectral characteristics (e.g., Giardini et al., 2020), revealing seismo-tectonically active zones (e.g., Stähler et al., 2022).The majority of the observed marsquakes are low-magnitude high-frequency (HF) events with energy around and above 2.4 Hz, which have been interpreted as having crustal propagation paths and have helped determine the physical properties of the shallow subsurface (Karakostas et al., 2021;Lognonné et al., 2020;Menina et al., 2021).Although larger low-frequency (LF) and broadband (BB) marsquakes with mantle traversing paths are less commonly observed, these events have been key for constraining the deep interior structure of Mars, including upper mantle seismic wave speed (Drilleau et al., 2022;Durán, Khan, Ceylan, Zenhäusern, et al., 2022;Khan et al., 2021), the radius and composition of the core (Irving et al., 2023;Khan et al., 2022;Stähler et al., 2021), structure near the core-mantle boundary (Durán, Khan, Ceylan, Charalambous, et al., 2022), and mineral phase transitions in the deep mantle (Huang et al., 2022).Constraints on the crustal structure nearby (Kim, Lekić, et al., 2021;Knapmeyer-Endrun et al., 2021;Li, Beghein, Davis, et al., 2022) and away (Kim, Banerdt, et al., 2022;Kim, Stähler, et al., 2022;Li, Beghein, Lognonné, et al., 2022;Li, Beghein, McLennan et al., 2022) from the InSight lander have also been possible from analysis of BB events.Accurate models of crustal structure, including global variation in crustal thickness (Kim et al., 2023), may be particularly useful for modeling wave propagation and source properties of marsquakes.For an overview of the results of SEIS/InSight, see Lognonné et al. (2023).
Prior to the arrival of InSight on Mars, reverse faulting due to planetary contraction driven by secular cooling was expected to be one of the principal drivers of seismicity (Solomon et al., 1991).Thus far, however, seismic evidence of ongoing tectonic activity has come from the Cerberus Fossae extensional graben system which has been a prominent source of both HF and LF marsquakes (Stähler et al., 2022).Moment tensor estimates based on detailed waveform fitting of body waves from Cerberus Fossae events S0173a and S0235b implied extension along steeply dipping normal faults with roughly E-W strike orientations (Brinkman et al., 2021).Alternatively, based on body wave polarities and relative amplitudes, Sita & van der Lee (2022) found that S0173a can best fit as a marsquake doublet that starts as a thrust followed by oblique normal faulting.They also found that S0235b likely represents vertical dip-slip motion with a small normal faulting component.Recently, Jacob et al. (2022) estimated moment tensor solutions from nine high-quality LF and BB marsquakes located in Cerberus Fossae and surrounding regions using P and S waveform fits, secondary phase amplitudes, and upper bounds on surface wave amplitude.They found that two of the events (S0325a and S0784a) located south of Cerberus Fossae near the Martian crustal dichotomy boundary are consistent with thrust faulting source mechanisms.Additionally, they located the hypocenters of all nine tectonic events to moderate depths in the crust (15-36 km).
On 4 May 2022, InSight recorded the seismic event S1222a, the largest marsquake observed during the mission, with an initial estimated moment magnitude of   Ma   4.6 (InSight Mars SEIS Data Service, 2019).Despite occurring in the early afternoon during a period of high background noise due to strong winds, InSight's very broadband seismometer, SEIS-VBB (Lognonné et al., 2019), recorded some of the clearest seismic waveforms of any marsquake observed during the InSight mission (Figure 1).S1222a was designated a quality A event by the Mars Quake Service (MQS) due to its high signal-to-noise ratio and the identification of both body and surface wave phases with clear polarization (Kawamura et al., 2022).Several prominent instrument glitches are present during the event (Figure 1d).Although glitches have been shown to be capable of influencing seismic analysis and interpretation (Kim, Davis, et al., 2021), no glitches are apparent near the phases of interest in this study.Despite the high SNR of the event, P and S waves are relatively emergent in character compared to other high-quality marsquakes that have yielded focal mechanisms (e.g., Brinkman et al., 2021), which poses a challenge for determining the polarities of seismic phases.Importantly, both Love and Rayleigh waves were identified in the waveforms of S1222a, which could provide complementary constraints on the moment tensor that have not been available to other marsquake moment tensor inversions.

10.1029/2023JE007793
3 of 13 In this paper, we use seismic waveforms from both body waves (P and S) and, for the first time, surface waves (Rayleigh and Love) to estimate the best-fitting magnitude, depth, and focal mechanism of S1222a.We take advantage of available models of the crustal structure along the S1222a minor-arc paths, which were obtained by modeling the group velocity dispersion of Rayleigh and Love minor-arc phases (Beghein et al., 2022) and their overtones (Kim, Stähler, et al., 2022).

Moment Tensor Inversion
For S1222a, the MQS reported an origin time of UTC 2022-05-04T23:23:07 (±4.8 s) and an epicentral distance of 37 ± 1.6°, based on the S-P differential travel times (InSight Marsquake Service, 2023) and a set of 250 interior models compatible with P, PP, S, SS, and SCS travel times (Stähler et al., 2021).Based on the P-wave polarization, the back azimuth of the event was reported to be 101° with uncertainty bounds of 96°-112°.Analysis of minor and major-arc surface waves and overtone data yields a larger back azimuth uncertainty range, with a peak likelihood near 125° (Kim, Stähler, et al., 2022).Using narrow-band estimates of high-quality data, Panning et al. (2023) estimated the most likely back azimuth to be 140°.Figure 1a shows the likely source region of S1222a based on the uncertainties in distance and back azimuth.
We estimate the best-fitting moment tensor solution of S1222a by inverting the waveforms of P, S, Rayleigh, and Love wave phases recorded by SEIS-VBB.Even though moment tensor estimation is a linear inverse problem for a fixed velocity model, the determination of source depth is not.Optimal phase alignment between synthetic and observed waveforms is achieved by alignment based on cross-correlation, as commonly done in regional moment tensor inversions on Earth (e.g., Dreger et al., 2021).This alignment introduces an additional nonlinearity into the moment tensor inversion.Therefore, we choose to use a grid search approach.The synthetic waveforms are computed using Instaseis (van Driel et al., 2015), which allows rapid retrieval of pre-computed Green's functions based on AxiSEM waveform simulations (Nissen-Meyer et al., 2014).The velocity model used to compute synthetic waveforms is a modified version of the KKS21_GP model (one of the 250 models

10.1029/2023JE007793
4 of 13 used by MQS; Stähler et al., 2021), in which the upper 80 km is replaced with a radially anisotropic model based on fitting both Rayleigh and Love wave dispersion measurements of S1222a (Beghein et al., 2022) (Figure S1 in Supporting Information S1).The resolved positive radial anisotropy in the crust (V SH > V SV ) in the model of Beghein et al. (2022) is representative of the path average structure, although negative radial anisotropy has been observed near the InSight landing site (Li, Beghein, Wookey, et al., 2022).Isotropic models with strong crustal layering have also been shown to fit the observed Rayleigh and Love dispersion data of S1222a (Xu et al., 2023).The shear attenuation quality factor Q μ is assumed to be 600 in the crust and mantle (e.g., Giardini et al., 2020).
The event location and origin time are essential parameters for performing a moment tensor inversion.Here, we fix the origin time to the time reported by MQS and locate the event using the reported back azimuth of 101°.We assume an epicentral distance of 35.9° based on the body wave travel times predicted from our velocity model.While uncertainties in the location of the event will introduce uncertainties in the estimations of source properties, we find that our analysis is insensitive to small shifts in source location.
Prior to waveform fitting, raw data are instrument-corrected and converted to displacement, then rotated to the radial (R), transverse (T), and vertical (Z) components.Body wave data are filtered between 3-12 s and surface wave data are filtered between 14-36 s using a fourth-order nonzero-phase Butterworth filter.The number of cycles for these waves during their propagation from origin is therefore relatively low and about 40, 70, and 30 for P, S, and surface waves (see Table 1), justifying a model with low attenuation.However, scattering effects could still reduce the amplitudes.For each seismic phase, we manually select the windows in which to fit seismic waveforms.The window durations are 35 s for the P-wave and either 25 s or 40 s for the S-wave.The arrival of S on the Z component is unclear, so we omit it from the inversion.Minor arc Rayleigh waveforms are fitted on both Z and R and the Love waveform is fitted on T. Both Rayleigh and Love waves are fitted in 250 s long windows.
To account for small travel time shifts between observed and predicted data, synthetic waveforms are aligned by cross-correlation prior to calculating the misfit.The best-fitting double-couple solution for each possible depth is found by searching over a grid spanning the range of possible values of strike, dip, and rake, with 120 grid points in each dimension.M W is varied between 4.0 and 5.0 in increments of 0.05.This search is performed for source depths between 0-100 km, with a depth increment of 2 km.For every step in the grid search the misfit χ is calculated using Equation 1 following alignment of observed and synthetic waveforms that maximizes their cross-correlation: where   obs  and   syn  are the observed and synthetic displacement waveforms in the ith window, respectively, N is the number of windows, T is the window length, and w i is a weighting factor.We validate our approach by performing a single-station moment tensor inversion using terrestrial data (Figure S2 in Supporting Information S1).
Although the choice of weights is subjective, it is an important parameter in waveform-based moment tensor inversions.Additionally, limiting the maximum allowable cross-correlation shift for optimal phase alignment may provide additional constraints on the solution.Here, we take two different inversion strategies to explore the importance of these choices.In strategy 1, we apply weights to body and surface waves that are inversely proportional to the L2 norm of the observed data vector to ensure that each phase contributes approximately equally to the total misfit.The windows and weights used for the inversion are shown in Table 1.In this case, we use an S-wave window with a 25 s duration.Additionally, we allow large travel time shifts, Δt, between data and synthetics of ±10 s for P and S waves and ±20 s for Rayleigh waves.To honor the resolved radial anisotropy in the velocity model, we require the travel time shifts applied to Rayleigh and Love waves to be within ±5 s of each other.We refer to this value as the relative travel time shift, Δt R .In strategy 1, the onset time of the seismic phases is not required to be fit due to the arbitrary cross-correlation alignment that is allowed.This means that there is no requirement that the direct body waves of the synthetics are aligned with the observed direct waves, and may instead be aligned on later arriving phases such as depth phases.In strategy 2, we limit the permissible travel time shift of P and S waves to be within ±2 s of the travel time picks.Additionally, to emphasize the importance of the direct P-wave arrival, we amplify the weight in the first half of the P-wave window by a factor of 20 compared to the rest of the window.An extended S-wave window with a 40 s duration is used so that the importance of fitting depth phases is emphasized.The maximum allowable R1 shift is ±20 s.

Strategy 1
Figure 2 shows the results of the double-couple grid search over source depths of 0-100 km, using inversion strategy 1. Inversions were performed using both body waves and surface waves together (Figure 2a), body waves only (Figure 2b), and surface waves only (Figure 2c).From here on, we refer to inversions of both body and surface waves as BW + SW inversions, and inversions of only body waves or only surface waves as BW or SW inversions, respectively.The lowest misfits of BW + SW inversions are found for sources in the depth range of ∼20-28 km.In this depth range, the predominant focal mechanism of the best fitting BW + SW solutions is a thrust fault with a roughly NW-SE striking fault plane.In general, the misfit varies strongly with depth for both BW + SW and BW inversions, likely owing to the influence of depth phases.For SW inversions, the misfit varies more smoothly with depth, with the best-fitting solutions found for a source near 28 km.Here, the best fitting solution is an E-W striking normal fault.The discrepancy between solutions of BW + SW inversions and individual BW or SW inversions emphasizes the importance of considering all phases together.A second minimum in the misfit function near 4 km depth suggests that a shallow source cannot be ruled out from SW inversions alone.Events deeper than ∼50 km are considered unlikely because the misfit generally increases with increasing depth beyond this range for all inversion types.Below, we outline two scenarios for focal mechanisms at different depths in the crust based on the total misfit and detailed assessment of important aspects of the waveform fits such as body wave first motions.
The first scenario we consider is a surface source (red beachball in Figure 2a), which is shown in Figure 3a.Although a source depth at 0 km is not in the region of lowest misfit, we include it because it is representative of the fits achieved for a shallow source.The best-fitting focal mechanism for the BW + SW grid search is a M W 4.6 (scalar moment, M 0 = 1.0 × 10 16 Nm), predominantly reverse faulting event, with nodal planes striking either approximately E-W or NW-SE.This solution is consistent with an inversion based on analysis of direct wave amplitudes alone (Figure S3 in Supporting Information S1).The P-wave is fitted as the large pulse with an onset near 10 s in the P-wave window.Although the synthetics provide a good match to the amplitude of the first large positive arrival, the first arriving energy at the P-wave onset time (t = 0) is not fitted.This is because the signal near the P-wave onset time has a small amplitude, and thus will have a minor influence on the overall misfit in the P-wave window.Predicted S-waveforms align well with the observed signal and show general agreement in amplitude.The predicted Rayleigh waves are in excellent agreement with the observed amplitudes, although the early portion of the signal is not well fitted.Synthetic Love waves are in good agreement with both amplitude and phase of the observations.The absolute travel time shift of Rayleigh is ∆ t = 19.7 s and the relative shift between Rayleigh and Love waves is ∆ t R = −5.0s.Absolute travel time shifts of body or surface wave phases could reflect uncertainties in source location or errors in the velocity model between their different paths through the interior.Nevertheless, because our crustal velocity model already includes radial anisotropy that explains well the observed group velocities of minor arc Love and Rayleigh waves, focal mechanisms that imply small relative shifts between Rayleigh and Love should be preferred.
Figure 3b shows the second scenario we consider, which is for a source at 20 km depth (green beachball in Figure 2).The magnitude of the best-fitting solution is M W 4.3 (M 0 = 3.5 × 10 15 Nm) and the source represents reverse faulting along an NW-SE oriented fault plane.This solution provides a substantial improvement to the P-wave fit compared to the shallow source, with the amplitude and waveform of the largest pulse matched well on both the vertical and radial components.In this scenario, the direct P-wave is modeled as a small negative pulse, followed by a larger pP phase that arrives ∼6 s later.This is similar in character to the observed P-wave signal, suggesting that the largest pulse in the P-wave window is the depth phase (Figure S4 in Supporting Information S1).Although the small P-wave pulse strongly resembles the observed signal near the P-wave onset, the phase arrives slightly later than is observed.The S-waveform fits also show good agreement with observed amplitudes and appear to satisfy the first motions.The synthetic Rayleigh waves closely match the observed phase and amplitude near the signal onset, and the travel time shift is ∆ t = 11.6 s.The Love wave fit is comparable to the previous scenario.

Strategy 2
Inversion results using strategy 2 are shown in Figure 4.In this case, sources shallower than ∼16 km depth are strongly disfavored by the inversion because the direct P-wave of the synthetic is always aligned closely with the observed P-wave onset and the larger secondary pulse is not able to be well fitted.Again, the best fitting solutions for BW + SW inversions are found for a NW-SE striking thrust fault, although in this case there is a clear preference for a source near 22 km depth.For both BW + SW and BW inversions, oblique thrusting or predominantly strike-slip mechanisms are preferred at depths between ∼30-50 km.Although solutions near 40 km provide a low level of total misfit, we do not prefer these solutions since they are disfavored by surface wave observations.The waveform fits for the best-fitting thrust fault solution at 22 km depth (green beachball in Figure 4a) are shown in Figure 5.In this case, both the low amplitude signal near the P-wave onset and the early portion of the larger subsequent arrival are well fitted.Additionally, signals in the extended S-wave window are well fitted as depth phases (sS) on both the radial and transverse components.Although we prefer to exclude the vertical component S-wave in the inversion, we find similar results if it is included (Figure S5 in Supporting Information S1).
Modeled surface waves are similar to those shown in Figure 3b (using a thrust fault solution at 20 km depth) and provide excellent fits to both Rayleigh and Love.The scalar moment of this solution is 2.5 × 10 15 Nm (M W 4.2).

Discussion
Our identification of body wave depth phases and the inclusion of surface wave constraints resulted in an optimal source depth near 22 km for S1222a; additional information from the frequency content and signal duration may 10.1029/2023JE007793 9 of 13 also help reduce uncertainty in depth determination.Giardini et al. (2020) proposed that LF events likely originate in the uppermost mantle, which could explain the lack of high frequencies that would be attenuated away in the mantle.Since then, other observed marsquakes have challenged this paradigm.For example, high frequency energy (>5 Hz) has been observed from teleseismically detected large impacts (Kim, Banerdt, et al., 2022;Posiolova et al., 2022) with mantle traversing body waves.
The signal duration may provide an additional depth discriminant because shallow events are likely to produce longer coda due to extensive scattering in the near-surface (e.g., Karakostas et al., 2021;van Driel et al., 2021).S1222a has one of the longest durations of any recorded marsquake, lasting more than 8 hr for the multiple orbit surface waves, which orbit up to three times (Kim et al., 2023), and about an hour in the body wave coda (see Figure 1), suggesting a shallow origin.However, without more precise knowledge of the attenuation and scattering along the path, this remains fairly qualitative speculation.For Cerberus Fossae events, the corner frequency may be indicative of source depth because the lower corner frequencies of LF events compared to HF events are thought to result from their origin in deeper, warmer, structurally weaker zones (Stähler et al., 2022).Similarly, shallow moonquakes exhibit higher corner frequencies than deep moonquakes, likely indicating a high-stress drop (Oberst, 1987).Events outside of Cerberus Fossae generally show higher corner frequencies, making the relationship between source depth and corner frequency less clear.The high corner frequency observed for S1222a (∼4 Hz, see Kawamura et al., 2022) could result from a combination of high stress drop and a cold, weakly attenuating lithosphere.
If the source location of S1222a were known precisely, our source mechanism inferences could help inform our understanding of the relationship between geologic setting and crustal deformation regime.However, due to the relatively large uncertainty of the back azimuth, the uncertainty bounds of the source region span a broad range and include areas of both the northern lowlands and southern highlands (Figure 6).According to the mapped fault  catalog of Knapmeyer et al. (2006), the likely source region of S1222a contains both compressional faults and extensional faults.Reverse faults in the source region are clustered in the northern lowlands with predominantly N-S strike orientations, aligning with the wrinkle ridges of Avernus Dorsa.The formation of this wrinkle ridge system has been dated to the late-Noachian to early Hesperian period (Tanaka et al., 2014), which was likely a time of rapid global cooling and contraction (Watters, 1993).The strike orientation of the best-fitting thrust solutions of S1222a tends to be NW-SE, which is in possible agreement with the crustal structures expected to be associated with Avernus Dorsa.For example, the focal mechanism shown in Figure 5 has a strike orientation rotated 38° from the N-S direction and appears to be roughly aligned with mapped fault segments (see beachball in Figure 6).Mechanical modeling of wrinkle ridge formation has shown that reverse faults can extend deep in to the crust and upper mantle, with a wide range of possible dip angles (e.g., Andrews-Hanna, 2020).Thus, it is plausible that S1222a resulted from reverse faulting on a moderately dipping blind thrust associated with the Avernus Dorsa wrinkle ridge system, that is currently active under a regime of compressional tectonics due to planetary thermal contraction.
Further south, near the dichotomy boundary, normal faults with a NW-SE strike direction are common within the source uncertainty bounds.Additionally, the Al Qahira Vallis, an approximately N-S oriented graben system, is located near the southern end of the source uncertainty region.Although the location of S1222a reported by MQS (Kawamura et al., 2022) is in the northern lowlands, there are several lines of evidence that are suggestive of a source in the southern highlands.First, the back azimuth estimate based on minor and major-arc surface waves (e.g., Kim, Stähler, et al., 2022) differs from that based on the direct P-wave polarization and suggests a back azimuth near 125°, which would place the source in the southern highlands near the Al Qahira Vallis.Second, the average crustal thickness inferred from surface waves along the S1222a path is ∼60 km which is inconsistent with crustal propagation primarily in the northern lowlands, where thinner crust is expected (Kim et al., 2023).Third, scattered wave analysis by Menina et al. (2023) suggests a thick crust (e.g., 60 km) in the source region of S1222a, consistent with a source in the southern highlands.There are no mapped reverse faults in the southern highlands within the source uncertainty bounds that align with the orientation of the preferred thrust fault solution, so if S1222a was located in the southern highlands, it may not have taken place on a fault with a clear surface expression.However, the orientation of several mapped normal faults south of the dichotomy boundary is roughly consistent with the orientation of best-fitting thrust fault mechanisms, although the fault slip of S1222a is opposite of what is expected for an extensional feature.It is possible that existing dichotomy boundary faults that were once extensional have been reactivated in a reverse regime under thermal contraction from planetary cooling (Jacob et al., 2022).
Although a reverse fault is preferred in our inversion, it is not possible to entirely rule out a normal fault mechanism.We find that if a stricter limit of ±5 s is placed on the maximum allowable shift of R1, a normal fault with a NW-SE strike orientation is preferred (Figure S6 in Supporting Information S1).However, given the uncertainty on both the origin time of S1222a and the velocity model used in the inversion, there is no strong justification for requiring such small travel time shifts, especially given that imposing the stricter requirement comes at the expense of reducing the quality of the body wave fits.
The structural model used to generate Green's functions is a key component of moment tensor inversions.Ideally, 3D models that account for structural variations near the source and receiver would be used so that reflected and converted body wave phases can be accurately modeled.However, in this study we are limited to one-dimensional (1D) structural models, as are commonly used for moment tensor inversions on Earth (e.g., Dreger & Helmberger, 1993).The structural model used here is based on fitting fundamental mode dispersion of Rayleigh and Love waves (Beghein et al., 2022), although our results are similar if we use the model of Kim, Stähler, et al. (2022), which fits both fundamental mode and overtone data (Figure S7 in Supporting Information S1).Although using models constructed from surface wave dispersion enables accurate fitting of Rayleigh and Love waveforms, the model does not incorporate constraints on crustal layering below InSight.
A structural model based on joint inversion of surface wave dispersion and receiver functions could potentially improve moment tensor estimations because waveform predictions would include converted and reflected body wave phases from discontinuities below InSight.Near-source and along-path structure, however, would remain unaccounted for and would require the calculation of Green's functions from models with lateral variations in velocity and crustal thickness (Kim et al., 2023).
In this study, we find moment tensor magnitudes for preferred focal mechanisms range between M W = 4.2-4.3(M 0 = 2.5 × 10 15 -3.5 × 10 15 Nm), which is lower than the initial MQS estimate (   Ma  = 4.6).The magnitude misfit between this inversion and the MQS magnitude based on Böse et al. ( 2021) is similar to the discrepancy found by Brinkman et al. (2021) for Cerberus Fossae marsquakes.Since the MQS model was calibrated using all possible focal mechanisms, it is possible that events with specifically high body wave radiation in the direction of the station are overestimated in magnitude, which can be resolved by a waveform-based source inversion, as presented in this paper.
Over the course of the InSight mission, the majority of significant marsquakes have been related to extensional tectonics in the Cerberus Fossae system, and compressional faulting due to planetary contraction did not appear associated with the observed seismic activity.Stähler et al. (2022) estimate that the Cerberus Fossae events account for an annual seismic moment release of 1.4-5.6 × 10 15 Nm/yr, or over half of the seismic moment release in the InSight hemisphere.The estimated scalar moment of S1222a (2.5 × 10 15 -3.5 × 10 15 Nm) suggests that this event alone accounts for a large fraction of the annual seismic budget of Mars, estimated to be between 1.55 × 10 15 and 1.97 × 10 18 Nm/yr (Knapmeyer et al., 2023).If S1222a did result from compressional faulting, it provides evidence that planetary thermal contraction is an important source of seismicity on Mars, as was expected prior to the arrival of InSight.

Conclusions
Based on waveform fits from P, S, Rayleigh, and Love wave phases, we estimated the source properties of S1222a, the largest marsquake recorded during the InSight mission.Our approach, which included minimizing the L2 misfit between broadband observations and synthetic seismograms via a grid search, allowed us to estimate the best-fitting magnitude, depth, and focal mechanism.We find that S1222a likely resulted from reverse faulting along a NW-SE oriented fault plane at moderate depth in the crust.Potential depth phases suggest a source depth near 22 km.The presence of active thrust faulting on Mars indicates that thermal contraction of the martian crust may be an important source of seismicity.The estimated moment magnitude of S1222a is M W = 4.2-4.3(M 0 between 2.5 × 10 15 -3.5 × 10 15 Nm), lower than initially determined by MQS.The source uncertainty region of S1222a includes mapped compressional faults in the northern lowlands associated with the Avernus Dorsa wrinkle ridge system, suggesting that S1222a may be associated with these tectonic structures.However, a source in the southern highlands, where predominantly extensional faults have been mapped, cannot be ruled out.Future work may be able to provide tighter constraints on the source parameters of S1222a (and other marsquakes) if improved structural models of Mars become available or if additional seismic phases are clearly identified.Our moment tensor inversion was based on a modified version of the MTUQ code package (https://github.com/uafgeotools/mtuq), and we thank the developers for making this tool openly available.Seismic data from the S1222a event, and other marsquakes, are available for download through the IRIS Data Management Center.This paper is InSight contribution number 290.The authors acknowledge the NASA, the CNES, their partner agencies and Institutions (UKSA, SSO, DLR, JPL, IPGP-CNRS, ETHZ, IC, and MPS-MPG) and the flight operations team at JPL, SISMOC, MSDS, IRIS-DMC, and PDS for providing the SEED SEIS data.D.K. acknowledges support from the ETH+ funding scheme (ETH+02 19-1: "Planet Mars").N.S. and V.L. were supported by NASA Grant 80NSSC18K1628.C.B. and J.L. were funded by NASA InSight PSP Grant 80NS-SC18K1679.P.L acknowledges support from CNES and ANR (ANR-19-CE31-0008-08 MAGIS, ANR-18-IDEX-0001).JCEI acknowledges support from UKSA Grant ST/ W002515/1.We thank Göran Ekström and Julien Thurin for their thorough and constructive reviews that helped improve this manuscript.

Figure 1 .
Figure 1.(a) Map showing the location of InSight (green triangle) and the Mars Quake Service (MQS) reported location of S1222a (blue star) (Kawamura et al., 2022).The dashed blue lines depict the range of uncertainty in epicentral distance (37.0° ± 1.6°) and back azimuth (96°-140°).Thin white and black lines indicate mapped compressional and extensional faults, respectively, from the database of Knapmeyer et al. (2006).The bold line separating the northern lowlands and southern highlands marks the approximate dichotomy boundary.(b) Raw recordings of S1222a from the SEIS-VBB instrument.The dashed vertical lines show the MQS arrival times of P and S. (c) Same as (b) but the data are instrument corrected and rotated to vertical (Z), radial (R), and transverse (T) components using a back azimuth of 101°.All components are bandpass filtered between 0.1 and 0.6 Hz.(d) Scalogram of BHZ data.Vertical green lines on the lower time axis mark clearly apparent instrument glitches.

Figure 2 .
Figure 2. Inversion misfit for a range of source depths using strategy 1. Panels (a, b, and c), show results for BW + SW, BW, and SW inversions, respectively.In each panel, the top figure shows the normalized misfit for the full grid search depth range, and the bottom panel shows the best-fitting moment tensor solutions for source depths between 0 and 50 km.The red and green focal mechanisms in (a) represent the two scenarios shown in Figure 3.The horizontal dashed line marks the boundary between the crust and mantle according to our velocity model.

Figure 3 .
Figure 3. Best-fitting moment tensors and the corresponding waveform fits shown for sources at the surface (a), and at 20 km depth (b), using inversion strategy 1. Observed waveforms are shown in black and synthetic waveforms of BW + SW inversions are shown in red.Beachball solutions shown in the top left of each panel represent the best-fitting moment tensor for BW + SW inversions.The label to the upper left of each waveform indicates the phase (P,S,R1,L1) and component of ground motion (Z,R,T).In the P-wave and S-wave windows, t = 0 corresponds to phase onset times.

Figure 4 .
Figure 4. Inversion misfit and best-fitting focal mechanisms for a range of source depths using inversion strategy 2. (a) Shows results for BW + SW inversions and (b) shows results for BW inversions.SW inversions are not shown because they are the same using strategies 1 and 2.

Figure 5 .
Figure5.Preferred best-fitting moment tensor and the corresponding waveform fits shown for a source at 22-km depth, using inversion strategy 2.

Figure 6 .
Figure 6.Map of the source region of S1222a.The beachball represents the best fitting solution for a source at 22 km depth found using inversion strategy 2 and is shown at the location reported by Kawamura et al. (2022).The blue dashed lines illustrate a range of back azimuths between 96°-140°, and an epicentral distance of 35.4°-38.6°.White and black lines represent compressional and extensional faults, respectively, from the database of Knapmeyer et al. (2006).
Note.Diff start time is the differential time between middle of the window and quake origin time.Number of Q cycles are computed with central periods of 7.5 and 25 for body and surface waves.