Structure Along the Martian Dichotomy Constrained by Rayleigh and Love Waves and Their Overtones

Using seismic recordings of event S1222a, we measure dispersion curves of Rayleigh and Love waves, including their first overtones, and invert these for shear velocity (VS) and radial anisotropic structure of the Martian crust. The crustal structure along the topographic dichotomy is characterized by a fairly uniform vertically polarized shear velocity (VSV) of 3.17 km/s between ∼5 and 30 km depth, compatible with the previous study by Kim et al. (2022), https://doi.org/10.1126/science.abq7157. Radial anisotropy as large as 12% (VSH > VSV) is required in the crust between 5 and 40 km depth. At greater depths, we observe a large discontinuity near 63 ± 10 km, below which VSV reaches 4.1 km/s. We interpret this velocity increase as the crust‐mantle boundary along the path. Combined gravimetric modeling suggests that the observed average crustal thickness favors the absence of large‐scale density differences across the topographic dichotomy.

• By jointly analyzing Rayleigh and Love waves, and their overtones in the S1222a record, we obtain the seismic velocity structure of Mars down to 90 km depth • Radial anisotropy up to 12% (V SH > V SV ) is required in the crustal structure along the path to S1222a • Absence of large-scale density differences across the martian dichotomy better explains the average crustal thickness along the propagation path

Supporting Information:
Supporting Information may be found in the online version of this article.
fixing the hypocentral depth and location . However, the limited frequency content and the absence of Love waves in the record prevented crustal structure below 30 km and anisotropy associated with the Martian crust to be constrained. In addition, the paths of the minor-arc Rayleigh waves (R1) were mostly confined to the crust in the northern lowlands in the vicinity of Elysium . A weaker major-arc Rayleigh wave (R2) arrival was also reported, but due to unclear polarization of the data, the interpretation of structural constraints provided by this seismic phase was limited.
The marsquake S1222a, the largest event recorded during the mission to date (M w Ma 4.7; Kawamura et al., 2022), provides an important opportunity to further constrain lateral variations in crustal structure using surface wave analysis. Based on the epicentral location estimated by the Marsquake Service (MQS; InSight Marsquake Service, 2022), the surface waves identified in S1222a travel along the dichotomy on Mars (Figure 1a). Compared to the recordings of the two large impacts, the identified surface waves in S1222a have a broader frequency content and include not only Rayleigh but also Love waves, overtones, and multiple-orbit surface waves. The expected sensitivity of the dispersion measurements extends down to the upper mantle. Fortuitously, the propagation paths for the surface waves are close to the dichotomy boundary, where we expect crustal thickness variation to be greatest (Wieczorek et al., 2022) (Figure 1a).
Here, we report robust group velocity measurements of Rayleigh and Love waves and their first overtones for event S1222a. Using the available frequency content from the surface waves, we invert the group velocities to obtain profiles of S-wave velocity and radial anisotropy down to 90 and 50 km depth, respectively. We compare our results with previously published models derived from group velocity dispersion of fundamental-mode Rayleigh waves based on S1000a and S1094b  and discuss implications for lateral variations in crustal structure across the topographic dichotomy on Mars.

Surface Wave Dispersion Measurements
We use the seismic recording of S1222a and remove glitches from the 20 sample per second UVW channels of the Very Broad Band sensor (Scholz et al., 2020). U, V, and W denote the three non-orthogonal components of the seismic sensor (Lognonné et al., 2019). The deglitched records are rotated to up-down, north-south, east-west components and we confirm that the seismic waveforms are not strongly affected by any known electro-mechanical noise by the sensor or the lander Kim, Davis, et al., 2021b). The P-and S-arrivals are prominent both in the time and spectral domains across the 1-10 s period range and have been assigned a picking uncertainty of ±0.5 s and ±2 s by MQS, respectively. Using the differential travel time and P-wave polarization estimates from these body waves, MQS located this event at a distance of 37 ± 1.6° and a back azimuth range between 98° and 121° from the lander near Cerberus Plains ( Figure 1). About 200 s after the S-arrival, strong dispersive arrivals are evident in vertical-and horizontal-component data and are identified as minor-arc surface waves by the MQS (InSight Marsquake Service, 2022). In contrast to previously reported surface waves in S1094b , the S1222a recording shows both minor-arc Rayleigh (R1) and Love waves (G1) on the vertical and transverse components, respectively. Assuming that propagation occurs along the great circle path (GCP), the apparent time delay between R1 and G1 suggests an anisotropic structure on Mars. Overtones and multiple orbit surface waves have also been detected and cataloged by the MQS. See Kawamura et al. (2022) for more detailed information on event description.
To make group velocity dispersion measurements on the identified surface wave arrivals, we employ a single-station approach using a multi-wavelet transformation as a filter bank (Poppeliers & Preston, 2019). Because the wavelet transform optimizes the trade-off between time and frequency resolution compared to typical narrow-band filtering, this method achieves stable, high-resolution measurements, while providing robust error estimates (Preston et al., 2020). Here, we focus on a 400 s-length window around the surface wave arrivals and use 10 mutually orthogonal wavelets to compute 10 dispersion estimates across the 5-50 s period range of the vertical-and horizontal-component waveforms. Long-period energy beyond 50 s is visible but due to presence of strong atmospheric noise, we focus on periods <50 s that have higher signal-to-noise ratios (SNR) . For each period, we normalize to unity the power of the resulting transform, and pick the maximum envelope amplitude for each of the 10 transforms across different periods (e.g., Figure 2). Picks on vertical and (a) Location of the event S1222a (black symbol). The lander location is denoted by the yellow symbol. The great circle paths for minor-arc (R1, black) and major-arc Rayleigh waves (R2, gray) in S1222a are in solid, while those paths including the back azimuth uncertainty are displayed in dashed line. Location of two large meteoroid impacts (S1094b and S1000a) and the corresponding paths for previously identified surface waves are based on Posiolova et al. (2022) and Kim et al. (2022), respectively. (b) Broadband three-component (ZRT) seismogram of S1222a (light gray) with P and S wave picks. Rayleigh and Love waves are clearly visible on the data bandpass filtered between 10 and 100 s. R1_1 and G1_1 denote the first overtones of R1 and G1, respectively. (c) Vertical and (d) transverse component S-transforms show large amplitude surface wave arrivals with dispersion. Time after origin of the event is converted to group velocity using the equatorial radius of Mars (purple ticks with labels at the top of each spectrogram). radial components are collected for Rayleigh waves while those on the transverse component are used for Love waves. Next, these initial picks are filtered based on the back azimuth and polarization analysis described below.

Back Azimuth and Polarization Analysis
To obtain robust dispersion measurements for inversion, we implement two additional steps that discard those measurements that substantially deviate from the propagation direction and exhibit particle motion inconsistent with that expected for surface waves. For Rayleigh waves, we perform a grid-search to find back azimuth estimates that maximize correlation between the radial-component and the Hilbert transform of the vertical component (i.e., maximize elliptical particle motion in the vertical plane; see Figure 2a). For Love waves, we apply a similar grid-search to the analysis window around G1 and its overtone as we minimize the ratio between the average power of vertical and transverse component data. We find consistent back azimuth estimates between body and surface waves with a small offset of ∼10° possibly indicating complexities associated with lateral varying structure along the wave propagation paths. However, we do not have sufficient sensitivity to resolve such a small back azimuth difference observed in the data.
Next, we conduct frequency-dependent polarization analysis (Park et al., 1987) on the S1222a waveforms to investigate the particle motion of the surface wave arrivals. We employ the S-transform (Stockwell et al., 1996) of the three-component event waveforms and compute a 3 × 3 cross-component covariance matrix at each frequency in 90% overlapping time windows whose duration varies inversely with frequency. The relative sizes of the eigenvalues of this covariance matrix are related to the degree of polarization of the particle motion, while the complex-valued components of the eigenvectors describe the particle motion ellipsoid in each time-frequency window. Our computed polarization attributes (see Tables S2-S1 of Stähler et al. (2021)) are then combined into a metric which highlights signals with elliptically polarized energy in the vertical plane for Rayleigh waves (e.g., Kim et al., 2022;Figure 2b). For Love waves, we examine the phase angle of the particle motion ellipse to ensure our picks have particle motion that is dominantly polarized in the horizontal plane. We discarded picks that show deviations away from the expected propagation direction larger than the measurement uncertainty  or have irregular polarization. The remaining picks resulting from both back azimuth and polarization analyses are considered as our final measurements for inversion ( Figure 2c). This includes four dispersion curves for R1, G1, and their corresponding overtones in the 8-40 s period range. Unlike R1, the direction of propagation and particle motions of the identified R2 and R3 are largely scattered and unclear due to low SNR. (c) Summary of Rayleigh and Love wave dispersion measurements. R1_1 and G1_1 denote the first overtones of R1 and G1, respectively. Measurements from the off-great-circle propagation with low elliptical polarization (less than 0.1 above the global average) are discarded. Uncertainty of 0.1 km/s is assigned to our measurements to account for the MQS epicentral uncertainty of the event. (d) Distribution of the back azimuth estimates for the surface wave arrivals analyzed in this study. The maximum whisker length is specified as 1.0 times the interquartile range. Outliers beyond the whiskers are denoted by circle symbols. Note, the back azimuth estimates for the R2 and R3 arrivals are more scattered compared to those from the minor-arc surface waves due to low signal-to-noise ratios of the data.
Based on our analysis, only a few measurements of R2 and R3 are available at ∼35 s (e.g., Panning et al., 2022) thus we do not use the suggestive R2 and R3 arrivals directly in the inversion.

Inversion of Surface Wave Dispersion Data
We invert the group velocity measurements summarized in Figure 2c using a Markov chain Monte Carlo (McMC) method for sampling the path-averaged S-wave velocity structure with an adaptation of the Metropolis-Hasting algorithm (Hastings & Keith, 1970). We assume a fixed V P /V S ratio of 1.81 estimated for the upper crust beneath InSight using a free surface transform matrix (Kim, Lekić, et al., 2021a). The scaling between V P and density is based on Birch's law (Birch, 1961). We employ a fixed parameterization strategy using b-splines as described in McMC Approach 1 of Kim et al. (2022). We parameterize the crust using eight b-spline functions overlying a mantle halfspace with a constant velocity. The depth to this constant-velocity layer is allowed to vary between 30 and 70 km, the depth range estimated for the average crustal thickness of Mars (Wieczorek et al., 2022). We consider uniform prior distributions for those spline coefficients for V S and radial anisotropy, that is, (V SH −V SV )/ V S in the inversion, but present the anisotropy as = ( SH∕ SV) 2 for easier comparison to other studies. Positive anisotropy corresponds to 1 . Highest-accuracy computation of dispersion in a transversely isotropic medium would require the specification of 5 elastic parameters; here, we assume V PV = V PH , anellipticity or = 1, and use V SV and V SH to compute the dispersion for Rayleigh and Love waves, respectively. This choice dramatically improves computational efficiency, while having negligible impact on accuracy in the period range used (Jiang et al., 2018; see also Beghein et al., 2022 for comparison of different parameterization schemes employed for anisotropic inversion). We compute chi-squared misfit between predicted and observed group velocity dispersion curves assuming a measurement uncertainty of 0.1 km/s based on the MQS event uncertainty . We vary both the assumed a priori distribution and the total number of McMC iterations to ensure that our final model is not biased by these choices. The Rayleigh and Love wave group velocity kernels calculated using the mean posterior velocity model yield the sensitivity of available group velocity measurements is weak at depths shallower than 5 km, similar to the case with S1000a and S1094b . However, the longest-period R1 and overtones in S1222a extend the previously reported sensitivity (5-30 km) to 90 km depth (50 km depth for G1) (Figures 3a and 3b). . Kernels are computed based on the average model from the inversion. R1_1 and G1_1 denote the first overtones of R1 and G1, respectively. Mean predicted dispersion curves are denoted by gray lines. (C) Posterior distribution of V SV and (D) radial anisotropy structure inverted from the group velocity dispersion curves of S1222a. Posterior distribution and prediction are based on the best-fitting 10,000 models after two million iterations. Depths where sensitivity is inadequate (<40% in cumulative kernel strength) are muted. Note our V SV is constrained by a combination of both Rayleigh and Love waves while the radial anisotropy is primarily constrained by Love waves. Posterior distributions from the isotropic inversion of S1000a (light blue) and S1094b R1 (magenta) are denoted by horizontal lines at each depth  10.1029/2022GL101666 6 of 10

Results and Discussion
Our V SV profile is characterized by a positive velocity gradient of 0.015 km/s per km with an average velocity of 3.45 km/s between 5 and 60 km depth (Figure 3c). Like models derived from Rayleigh wave group velocity measurements from the two large impacts (S1000a & S1094b), which have the average V SV of 3.2 km/s between 5 and 30 km depth , the S1222a models also show little depth-variation in V SV with a slightly slower average V SV of 3.17 km/s in the overlapping sensitivity depth ranges. Below 30 km depth, the posterior distribution of S1222a V SV is shown to be compatible with the broader distribution of S1000a V SV (cyan, Figure 3c), approaching 3.8 km/s at 60 km. We observe a large discontinuity in the V SV profile at 63 ± 10 km depth with a velocity jump of ∼0.4 km/s (representing a total impedance contrast of ∼20%). This velocity jump accounts for the steep increase in the group velocity of R1 and its overtone seen near 25 and 12 s, respectively. Similar velocity increases at periods with sensitivity near the crust-mantle boundary are observed on Earth in both continental and oceanic settings (e.g., , 1952. Below the discontinuity, the average velocity of 4.1 km/s is consistent with that inferred for the upper mantle for Mars from body wave analyses (Duran et al., 2022;Khan et al., 2021). Therefore, we interpret the 63 km deep interface to represent the crust-mantle boundary along the R1 path. Its depth is well within global crustal thickness estimates on Mars (Wieczorek et al., 2022). The abruptness of the velocity jump across the discontinuity is in part due to how the mantle property is being parameterized by a constant value. Moreover, because the surface wave sensitivity functions are fundamentally broad over a wide depth range (e.g., Figures 3a and 3b) and the R1 path traverses near the dichotomy, we are unable to constrain whether the martian crust-mantle boundary is sharp or gradational ( Figure S4 in Supporting Information S1).
Identification of both R1 and G1 in S1222a allows us to determine radial anisotropy of shear wave speeds within the crust. We find that a model where radial anisotropy steadily decreases from ∼ 1.3 (equivalent to 12% for ( SH − SV) ∕ S ) at 5 km depth with a gradient of −0.01 per km depth is required to fit the dispersion measurements of the Love waves (Figure 3d). Our resolution test on the anisotropic inversion shows that the observed anisotropy in the model is resolvable within the sensitivity depth range of the surface wave data ( Figure S4 in Supporting Information S1). Indeed, no simple isotropic crust can explain the group velocities of G1 or its first overtone if we preclude the presence of anisotropy in our inversions ( Figure S5 in Supporting Information S1). While our posterior distribution of prefers values less than one (V SV > V SH ) in the lower crust (below 40 km) we do not believe this to be robust because we lose sensitivity below ∼50 km depth due to the absence of long-period Love wave dispersion measurements; this is reflected in the large model uncertainties between 50 and 60 km depth. Furthermore, we observe little change in the Chi-squared misfit when values for > 1 are allowed during the inversion ( Figures S6 and S9 in Supporting Information S1). Variation of data uncertainty does not strongly affect our inversion results either (Figures S7-S9 in Supporting Information S1).
To quantify the implication of our inverted models, including the observed crustal thickness in a global context, we generate Rayleigh wave phase/group velocity maps of Mars by linearly extrapolating our S-wave velocity profile based on the crustal thickness model constrained by gravity data (Wieczorek, 2021). Following the modeling steps described in Wieczorek et al. (2022), each crustal thickness model is produced by fixing the crustal thickness to 45 km at the location of the lander (Knapmeyer-Endrun et al., 2021) assuming a density contrast over the boundary of topographic dichotomy as mapped by Andrews-Hanna et al. (2008). We compute global crustal thickness models for spherical harmonics up to degree 120. Then, we linearly scale the velocity profiles in Figure 3c and from Kim et al. (2022) with the relative depth variations to the crust-mantle boundary. Motivated by the crustal models of Mars discussed in Wieczorek et al. (2022), two end-member dichotomy models are tested: (a) type I-a model with a uniform density ranging from 2,550 to 3,050 kg/m 3 (e.g., Baratoux et al., 2014;Wieczorek et al., 2019) and (b) type II-a model with a density contrast of 100-500 kg/m 3 across the dichotomy that may indicate crust that originated exogenically (e.g., Andrews-Hanna et al., 2008;Nimmo et al., 2008) (Figures 4a and 4b).
Due to a known trade-off between crustal thickness and density, the larger the assumed density contrast, the smaller the average crustal thickness variation across the dichotomy becomes. Hence, these models represent two limiting cases for describing the crustal dichotomy structure on Mars. We compute predicted travel times of the R1 as well as R2 and R3 arrivals by kinematic ray tracing of the surface waves through phase velocity maps at different periods (e.g., Equations 16.185 and 16.186 in Dahlen & Tromp, 1998). We find that the resulting travel times deviate by less than 1% between the GCP and the ray theoretical path for R1-R3 arrivals ( Figure S10 in Supporting Information S1), justifying our use of the great circle approximation in the inversions. Therefore, Figure 4. Conceptual models of the martian crust (a) with (type II) and (b) without a density contrast (type I) and the corresponding global group velocity maps at 40 s using the inverted velocity profile of surface waves in S1222a (c, d), respectively. Regions shaded indicate the off-propagation paths based on the back azimuth analysis. Inset below (c, d) shows crustal thickness profiles along the S1222a R1 for each type of models for all tested density values in Supporting Information S1 (Figures S11-S16). Horizontal solid and dashed lines indicate the average crustal thickness and its uncertainty observed in Figure 3c. (e, f) Average group velocities along the propagation paths for all of the models using the velocity profile in S1222a (K1222aR1G1; Figures S11-16 in Supporting Information S1) and (g, h) S1094b (K1094bR1; Figures S17-S22 in Supporting Information S1). Group velocity at each period increases as density contrast increases across the dichotomy boundary. Note the uncertainties of the R2 and R3 arrivals in the MQS catalog V12 (gray symbol) are substantially larger than those provided by this study (yellow symbol) due to the absence of back azimuth and polarization constraints. lateral thickness variations are unlikely to substantially affect the travel times of the surface waves particularly for long-periods >25 s. However this assumption may not be optimal for higher-orbit surface waves beyond R3 ( Figure S10 in Supporting Information S1) and a 3D sensitivity kernel should be taken into account rather than a simple range of GCPs drawn from the back azimuth uncertainties. More realistic 3D wavefield simulations through candidate crustal thickness models are beyond the scope of this study, but such an effort should provide further information on 3D surface wave propagation (e.g., Bozdag et al., 2017).
In the group velocity maps shown in Figures 4c and 4d, we gray out regions outside the vicinity of the GCP based on the back azimuth uncertainty of the event location. The R2 and R3 arrivals in our S1222a data travel with an average speed of 2.88 km/s traversing much greater distances along the potential velocity contrast across the dichotomy. This speed is 0.13 km/s faster than the R1 at 35 s period (Figure 2c). We find the spread in group velocities due to GCP propagation is substantially smaller than the range predicted by different crustal thickness models (Figures 4e-4h). For example, the type II model with a 300 kg/m 3 density contrast at 35 s implies an average crustal thickness of 43 km which yields ∼7% larger average group velocity for R2 than the uniform density model. Because the crust in type I model is thicker on average, the corresponding velocity profile is stretched downward to larger depths and the apparent speed at which the surface waves travel at a given depth is expected to be slower for such a type of model. The group velocity predictions for R2 and R3 arrivals in S1222a are strongly dependent on the choice of velocity profiles used in the modeling (Figures 4e-4h). While type I crustal model shows the best-fit between the group velocity predictions and measurements of R2 and R3 arrivals based on the previous velocity profile of S1094b, we were not able to explain the data with the new V SV profile of S1222a ( Figure 3c) and the predictions are largely under-estimated. If we assume the scaling approach in the global extrapolation is optimal, our analysis suggests that the low velocity structure beneath the lander (Knapmeyer-Endrun et al., 2021) which is also evident below the S1222a R1 path (Figure 3c) may not be prevalent along the equatorial dichotomy. Crucially, the average Moho depth of 63 km shown in our inversion is only compatible with the absence of a density contrast across the dichotomy (inset, Figures 4c and 4d), independent of the velocity profiles used in the extrapolation procedure (Figures S11-S22 in Supporting Information S1).
Regardless of the choice of models, the Hellas impact basin is expected to be a significant outlier, with velocities close to that of the mantle across different periods (Figure 4). Once corrected for the fraction of the path traversing Hellas, it has been suggested that the crustal wave speed at 5-30 km depth is similar between the northern lowlands and the southern highlands . Here, we provide another independent constraint indicating that there is no large-scale dichotomy in average crustal density. We also constrain the crustal velocity structure at those depths to be largely similar (difference less than 5%) with a caveat that the non-linearity of the surface wave sensitivity to depth is difficult to be implemented with such sparse data collected on Mars. At greater depths, the propagation path samples both with crust and mantle, and periods larger than 50 s using the higher-order multiple orbit surface waves would have to be further analyzed and reviewed (InSight Marsquake Service, 2022; Kawamura et al., 2022).