High‐ and Low‐Frequency Waveform Analysis the Marsquake of Sol 1222: Focal Mechanism, Centroid Moment Tensor Inversion and Source Time Function

The seismometer onboard InSight NASA Mars mission discovered a seismically active planet. We focused on the strongest event named S1222a (4 May 2022, Mw ∼ 4.7), which was recorded by the Very Broad Band sensors and associated channel ELYSE and is located 37.2° away from InSight. We use two different methods based on a point source approach for an elastic, horizontally layered medium to retrieve source parameters of S1222a. In the first case, the seismic moment tensor inversion of high‐frequency seismogram data is calculated using a matrix method for the direct waves. The process includes the generation of records in displacement using the frequency‐wavenumber integration technique. A method of inversion of the moment tensor of direct P‐ and S‐waves, less sensitive to path effects than reflected and transformed waves, is presented, which significantly increases the accuracy and reliability of the method. In the second case, tensors were calculated using common low‐frequency full‐waveform inversion and the tests to verify the plausibility of this solution obtained from the single station calculation were performed and the uncertainty estimations for inversions can be useful in future research.


Introduction
Seismic data recorded by the SEIS experiment (Lognonné et al., 2019(Lognonné et al., , 2020) ) onboard the InSight mission (Banerdt et al., 2020) have shown that Mars is seismically active with more than 1,300 events detected and cataloged by the InSight Marsquake Service (MQS) (Ceylan et al., 2022;Clinton et al., 2020;Giardini et al., 2020).The global seismic event rate sets Mars as moderately active, between the weak lunar activity and the terrestrial intraplate seismicity (Banerdt et al., 2020) and provided enough seismic records to perform the first structure inversions of the planet (see Lognonné et al., 2023 for a review of the results).On 4 May 2022, NASA's InSight Mars lander detected the largest "extraterrestrial" quake ever observed on another planet: an estimated magnitude M w 4.7 temblor that occurred on the 1,222nd Martian day, or Sol, of the mission (Kawamura et al., 2022).We can state that the double couple is very likely a valid approximation due to the small magnitude of the quake and the fact that only relatively long-period (LP) waves are used, while a 1D model is an assumption we know to be far from reality, especially due to the large scattering we do see on seismic waves (Karakostas et al., 2021;Lognonné et al., 2020;Menina et al., 2021).
The basis for quantifying the sources of seismic events is the seismic moment tensor (MT).The inversion of full waveforms for MT parameters is a widely used approach applicable to all scales: from small to large earthquakes.The accuracy and reliability of MT inversions depend on whether two major assumptions hold.First, it is assumed that the point source approximation is valid and second, that the effect of Mars's structure on seismic waves is modeled correctly.If either of these assumptions does not hold, the resulting MT may contain a large non-doublecouple component, even if the source mechanism is a double-couple.To comply with the point source approximation in the MT inversion presented here, seismic waves are used only with wavelengths longer than the fault plane dimensions.Also, the analytical expressions are drawn out relating the components of MT to the components of displacements in the immediate vicinity of the source.The second assumption provides challenges because Green's functions used in MT inversions depend on the structure between the source and receiver, which is not precisely known.Obviously and until more seismic stations constrain Mars, these two limitations will remain for SEIS data analysis.
The first study of the focal mechanisms of three well-recorded events on Mars (S0173a, S0183a, S0235b) is presented by Brinkman et al. (2021).They have shown a method that is adapted to the case of a single, multicomponent receiver and based on fitting waveforms of P and S waves against synthetic seismograms computed for the initial crustal velocity model.Later, Jacob et al. (2022) developed a method of seismic moment tensor solution and applied it to nine tectonic Marsquakes.For this purpose, body waves P and S waveforms were inverted, and amplitudes of secondary phases (PP, SS, PPP and SSS).More information about Marsquakes detected during the InSight mission can be found in for example, Giardini et al. (2020), Clinton et al. (2020), Stähler et al. (2021), Ceylan et al. (2022).
Therefore, the aim of this study is to determine the focal mechanism of high-and low-frequency waveform analysis of the Marsquake of Sol 1222.To achieve this goal, the authors consider two different methods.One of them uses only the direct waves which are less sensitive to path effects than reflected and transformed waves (Malytskyy, 2010(Malytskyy, , 2016;;Malytskyy & Mikesell, 2021).Second method (low-frequency full-waveform MT inversion) is based on ISOLA software (Sokos & Zahradník, 2008).Here is an opportunity to show the advantages and disadvantages of each of them.

Theory: Waveform Inversion
Two different approaches were chosen to address the problem of the unavoidable inaccuracy of seismic waves modeling: (i) focusing only on direct waves and (ii) using low-frequency inversion.The choice of these methods is associated with high-frequency and low-frequency analysis of waveforms to determine the focal mechanism of the Marsquake of Sol 1222.

Theory-High Frequency Analysis (i)
In the first case we propose to invert only the direct waves instead of the full field.An advantage of inverting only the direct P-and S-waves is that, compared to reflected and converted waves, they are less sensitive to the structural model used in the inversion.For example, waveforms of converted and reflected waves depend strongly on velocity contrasts below the source and receiver, and thus imprecise knowledge of subsurface structure will lead to inaccurate modeling.Waveforms of direct phases are less sensitive to subsurface layering, and scattering and may carry a less distorted imprint of the source.The advantage of choosing a matrix method for calculating synthetic seismograms is its ability to analytically isolate direct waves from the full wave field.In the earlier version of our method, as well as in most other MT inversions, waveforms at several seismic stations are simultaneously inverted (Malytskyy & D ' Amico, 2015; Malytskyy, 2016).Although much more information on the source should be contained in the waveforms from several stations, we show in our study that all the components of seismic moment tensor contribute to the waveforms at only one station and, at least theoretically, can be retrieved from them.
We use a point-source approximation, assuming the location and origin time proposed by Kawamura et al. (2022).Based on forward modeling, a numerical technique is developed for the inversion of observed waveforms for the components of seismic moment tensor M(t), obtained by generalized inversion (Malytskyy & Mikesell, 2021).
Our method enables us to obtain the solution of the focal mechanism using records at only one station.Moment tensor solution by waveform inversion consists of two parts: forward modeling (mathematic modeling of seismic Earth and Space Science 10.1029/2023EA003272 waves in layered media by matrix method) and inverse modeling (spectra of moment tensor components are calculated using a solution of generalized inversion).Thus, based on the Thomson and Haskell matrix method (Haskell, 1951;Thomson, 1950), we develop a new analytical approach to calculate synthetic seismograms on the upper surface of layered inhomogeneous media and present an approach for determining the seismic MT and source time function (STF) from observed waveforms.We use the point source model.The source is located inside a layer and is represented by seismic moment tensor (Malytskyy, 2016).

Forward Modeling
We consider the equation of motion in an elastic layered medium.The vector of displacement is represented through vector and scalar potentials (Malytskyy, 2016).The solution of the equation of motion for scalar and vector potentials is represented in the form of FBM (Fourier-Bessel-Mellin) integrals.As a result, we obtain components of displacements u (0)  z (t,r,ϕ) , u (0) r (t,r,ϕ) and u (0) ϕ (t,r,ϕ) on the free surface in the matrix form Malytskyy (2010) and Malytskyy and Mikesell (2021): where The source is located at r = 0 and represented by six independent seismic moment tensor components: φ and k are the station's azimuth and horizontal wave number respectively (Figure 1).Functions g i = g zi ,g ri ) T and g φi contain parameters of structure model and propagation effects between the source and the receiver.
Equation 1 is then expressed for the direct P-and S-waves at one station in the spectral domain (Malytskyy, 2016): where vector T contains the components of displacement of direct Pand S-waves, vector T consists of the six independent components of MT, and matrix K contains parameters of the structure model.Mathematical expressions for the parameters of the matrix K can be found in the article by Malytskyy and Mikesell (2021).

Inverse Modeling
Following Aki and Richards (2002), vector M can be obtained from Equation 2 in the matrix form (a generalized inversion): where the wave denotes complex conjugation and transposition, and 1 inversion.
The seismic moment tensor M(ω) is calculated using Equation 3and can be transformed into M(t) by applying the inverse Fourier transform: where STF(t) is the source time function and M ij is the seismic moment tensor.
In the article Vavryčuk and Kuhn ( 2012) is shown that the focal mechanism depends on the frequency range of the studied seismic waves.We obtained the results of seismic tensor solution (Equation 3) and time-independent focal mechanism using records at one station for only direct P-and S-waves.Thus, the inverse problem consists of determining six independent components of the moment tensor M under the conditions that the source location and structure model are known.

Low-Frequency Full-Waveform Moment Tensor Inversion (ii)
In this case the standard MT inversion was performed.This method (low-frequency full-waveform MT inversion) is based on ISOLA software (Sokos & Zahradník, 2008), and it is commonly used for local seismic networks.The method closely follows the description outlined in previous work (e.g., Křížová et al., 2016).We consider a seismic source at a point (Aki & Richards, 2002) where u stands for displacement which is changeable in time and could be represented as convolution * of moment tensor M and Green's tensor G inscribed using Cartesian coordinates p, q.Then MT could be represented by six elementary moment tensors According to AXITRA code (Bouchon, 1981;Coutant, 1989) The final MT could be written as Earth and Space Science 10.1029/2023EA003272 The M 1 to M 5 tensors in Equation 7 represent five DC (double-couple) focal mechanisms, M 6 is a purely isotropic MT so only the parameter a 6 carries information about the isotropic part of MT (a 6 = Tr(M)/3).
When we convolve elementary MT from Equation 7 with Green's function then we get elementary seismogram E i each for every elementary MT.We get overestimated linear inverse problem for displacements which could be solved by least-squares method (Tarantola, 2005) The comparison between real u and synthetic s seismograms during moment tensor inversion is also standard.We are trying to maximize the waveform fit between real (observed) and synthetic seismograms, here is represented by global variance reduction We calculate full MT and deviatoric MT.The deviatoric inversion is just a simplification when the parameter a 6 is set to zero.The deviatoric MT is an approximation and can be used only when a negligible isotropic component of MT is assumed.
The main uncertainty comes from insufficient knowledge of Green's function (which is connected with inaccurate crustal model and errors in event location) and 3D structure.To reduce the influence of input parameters uncertainties we try to use as lowest frequencies as possible.During the inversion the MT is searched together with centroid position and origin time.This inverse problem is linear in MT and is solved using the least square method but is non-linear in time and space, so a grid search is used.The low-frequency full-waveform MT inversion is used to minimize the errors of MT components in accordance with Vavryčuk (2007), where the authors mentioned that complete MT can't be determined only from the P and S wave amplitudes when we have data from a single three-component station.Limitations in determining the focal mechanism due to the use of only one station (Šílený et al., 2022) are taken into account for the low-frequency full-waveform inversion.To show how the obtained solutions are affected by using only one station we calculate synthetic seismograms where we added three more fictional stations (forward simulation).Synthetics correspond to the solutions given in Maguire et al. (2023).We made full-and deviatoric-MT inversion using one, two, three or four stations.The tests were made with noise free data and additional uncertainties estimation with random noise were also performed.

Inversion Results for S1222a
In this section, we present the inversion results for the S1222a event on Mars (2022-05-04, P-arrival 23:27:45, 3:54 LMST, M w 4.7, back azimuth 109°; Kawamura et al., 2022) which is located on Aeolis Southeast at 37.2°d istance from InSight.(station name XB.ELYSE).We used the S1222a waveform data (InSight Mars SEIS Data Service, 2019a, 2019b) provided by the VBB sensors (Lognonné et al., 2019(Lognonné et al., , 2020)).The raw data have only been corrected by the instrument transfer function and rotated.As it was mentioned above, the two different approaches were applied in this study.

Results-High Frequency Analysis (i)
Like a validation tests we first present the focal mechanism of the S0235b event on Mars (26 July 2019), located 25°from the epicenter (Brinkman et al., 2021).We compare the two methods: in the first, we propose to invert only direct waves (Malytskyy, 2016), and in the second we consider direct inversion for the full moment tensor (Brinkman et al., 2021).We tried three different source depths: 17, 32, and 56 km.The TAYAK velocity model was used (Jacob et al., 2022).The focal mechanisms for the source depth of 32 km shown in Figure 2 look very similar to each other.The focal mechanisms at depths of 17 and 56 km were also very similar.
(i) We use the direct P-and S-waves.Only the beginning of the P and S waveforms are needed and there were free of glitches (Lognonné et al., 2020;Scholz et al., 2020).Modification of the TAYAK velocity model (Lognonné et al., 2020) used in the inversions of waveforms and for calculation of seismic tensor is listed in Table 1.The velocity model includes the upper crustal model based on Lognonné et al. (2020).The original displacement seismograms are filtered correspondingly between 0.2 and 9.0 Hz (see Figure 3a).The durations of direct P-and S-waves at the station are estimated visually from the records and delays of the reflectionconversion phases at the respective epicentral distance and source depth are considered.For the event S1222a, we estimate them to vary between 1.1 and 4.5 s for the P-wave (Figure 3b), and 1.4 s and 3 s for the S-wave (Figure 3c) at the station XB.ELYSE.The highest frequency, on the other hand, is controlled by the assumption of the point source and corresponds to a wavelength larger than the linear dimensions of the fault (often less than 1 km in small earthquakes).We know the time when the direct wave (P and S) arrives.The duration of the direct wave segment in the record corresponds to the duration of the process in the source.
After the end of the direct wave, the record should be zero, unless other phases have arrived.Since we do not know the duration of the source, we cut off the segment of the direct wave in places where the amplitude of the direct wave is zero.
The component of the moment tensor resulting from the inversion of the direct P-and S-waves forms at only the station ELYSE (at epicentral distance of 37.2°) using Equation 3 and the corresponding focal mechanism are shown in Figure 4.
Figure 5 shows the fit of synthetic waveforms to the observed ones for direct waves with Qmu = 600 in the crust and mantle.Synthetics are calculated using data of moment tensor inversion for the direct P-and S-waves at a source depth of 22 km (see Figure 4b).They were computed in displacement for E, N, and Z components using the modified TAYAK velocity model.We also used Instaseis (van Driel et al., 2015) and AxiSEM (Nissen-Meyer et al., 2014).Prior to waveform fitting, we obtained the vertical Z, transverse T, and radial R components of synthetic and observed waveforms.Direct Pand S-waves were filtered between 3 and 12 s.We manually select the windows to fit seismic waveforms: the vertical Z and radial R components (Pz, Pr) for the direct P waves, and radial R and transverse T components (S R,   S T ) for the direct S waves.The alignment was done based on cross correlation between the vertical components (for the P wave) and the transverse component (for the S wave).Estimated scalar moment of S1222a is M 0 = 3 × 10 15 Nm.

Results-Low-Frequency Full-Waveform Moment Tensor Inversion (ii)
(ii) In this case we obtained the origin time and centroid depth using a grid search during MT inversion.As an initial condition the origin time 23:23:07 and epicenter 171°E; 5°S (Panning et al., 2022) was chosen and the depth was 22 km (same as in the previous case = High frequency analysis (i)).Although we also performed calculations in the 1D layered velocity model mentioned in  several frequency ranges, the solutions for the interval 0.028-0.072Hz were chosen to illustrate outcome from the full-waveform MT inversion.This range was used in Maguire et al. (2023).
For the S1222a event, for which a faulting origin is most likely.Data from the single station does not have enough information to constrain an isotropic signature in the data and we only use a 6 from MT Equation 8 like parameter which stabilized inversion and with no physical meaning except contribution to the final value of scalar seismic moment M 0 .The isotropic component for the full MT solution is most likely associated with differences between the selected model and the real averaged structure along the source-station path and might include furthermore the signature of scattering and other 3D seismic structures not modeled.Despite of this fact the tests with synthetic data showed that for single station MT inversion is suitable to calculate full MT instead of commonly used deviatoric MT inversion.It is clearly visible in Figures 6a and 6b, where we at first in forward simulation calculated input data (for the depth 6 km and focal mechanism with strike = 292°, dip = 36°, rake = 61°) and then obtained solution for low-frequency MT inversion.In Figure 6a for full MT inversion you can see how the focal mechanisms changed less for different depths than in case of deviatoric MT inversion (Figure 6b).Because input data were without noise it is obvious that the best solution with highest correlation is in both cases right and it is in the depth of 6 km.This is the reason why we prefer to show the results for the full MT inversion and consider them more realistic.As has been written above, calculation of MT could be simplified by setting parameter a6 = 0 in Equations 9-11 only in case when negligible isotropic component is assumed.In such cases solutions for full MT and deviatoric MT inversions should be considered as almost the same if the accurate crustal model and huge data set are available.Although deviatoric MT inversion used to be more stable and it is widely used for tectonic events on Earth, in this special case when the single station was available and due to glitches in seismograms we assumed the deviatoric MT inversion is not suitable for this data set.
The focal mechanism for full MT as a result of inversion is shown in Figure 6c at the depth of 20 km.(strike, dip, rake: 63°, 82°, 5°and    see the stability of the focal mechanism close to the best solution and differ less than 5°in each parameter strike, dip and rake.Figure 6c is only cutout from the more extensive grid search in depths 5-30 km and in many origin time shifts. Even though comparing the results of different methods does not guarantee that the identical solutions are correct, we try to calculate inversion with fixed mechanisms presented in Maguire et al. (2023).He used "cut and paste" method which is less affected by inaccuracies in velocity model because the P and S waves are during calculation separated from the surface waves, but the main idea of the calculation is similar like in this study.In Table 3 there are focal mechanisms for four solutions ordered from the best to the worst and the found depth for global and local variance reduction (correlation) maximum for these prescribed mechanisms is also presented here.
In the uncertainty tests we prescribed due to real noise variance in data samples 5.0e 15 in all cases and we make calculations for the 100 possible solutions across the depths used in inversions (for more details see Sokos and Zahradník (2013)).Here we show more interesting results in Figure 7 (results for our best solution and best solution in Figure 7a).In Figure 7b is obvious that we can easily obtained solution with opposite rake angle and reverse mechanism.The Figure 7c shows that deviatoric solution is more scattered than full MT solution (Figure 7b).Synthetic tests with more stations showed that solution for the single station ELYSE (Figure 7d) is more ambiguous than results obtained from the two or more stations (in this case is azimuth and distance: 129°+ 1,195 km; 52°+ 1,202 km; 254°+ 1,261 km; 289°+ 2,190 km).Just for illustration (Figure 8) we demonstrate that in the case of large distances between the seismic source and the receiver and noisy data we can obtain only partial agreement for real and synthetic seismograms.It is shown on the example of surface waves, which are crucial for our calculations.

Discussion and Conclusions
In our study, we first explored the possibility of retrieving source parameters of the S1222a event on Mars (2022-05-04, M w 4.7).Addressing the problem, we've chosen to invert only the direct waves instead of the full field.An advantage of the direct P-and S-waves consists in their much lesser distortion, if compared to reflected and converted waves, by inaccurate modeling of velocity contrasts.So, the advantage of the matrix method of calculating the wave field that we have developed is its ability to mathematically isolate direct waves from the full field.
We presented a method to obtain the moment tensor solution of the direct waveforms at only one station.The method that we used in this study is based on an inversion approach described in Malytskyy (2016).The inversion scheme consists of two parts: forward modeling (propagation of seismic waves in vertically inhomogeneous media is considered and a version of the matrix method for calculation of synthetic seismograms on the upper surface is developed); inverse modeling (spectra of moment tensor components are calculated using a solution of generalized inversion).So, as another option modification of ISOLA software was used and uncertainty tests were performed, although it is primarily intended for local and regional distances.
Recently, Maguire et al. (2023) used waveform fitting of body waves and surface waves to estimate the moment tensor solution of S1222a and found that either an E-W to SE-NW striking thrust fault or normal fault can explain the data.For the same structural model and frequency range used in this study and in Maguire et al. (2023) for low-frequency MT inversion waveform (within the estimated error) similar solutions were found.For example, the full-moment tensor solution we find here based on long-period waveform inversion (Figure 6c with error estimation in Figure 7a) closely resembles the reverse solution of Maguire et al. (2023).Additionally, both studies find a similar optimal source depth approximately 20 km.Brinkman et al. (2021) used waveform inversion of polarization of filtered body waves and found a NW-SE striking thrust fault solution best explains the data.This is similar to the direct P-and S-wave inversion solution we find here (Figure 4), although we find the best fitting fault planes strike NE-SW.Despite the range of possible moment tensor solutions of S1222a, it is encouraging that independent studies based on different methodologies, and using different structural models, point to similar solutions.However, further work should be done to understand the sources of uncertainty in single station  moment tensor inversions of marsquakes, which may help us understand the discrepancies between solutions and provide more robust constraints.
The differences among solutions obtained in this article are due to different approaches and frequency ranges.Method (i) does not use surface waves in contrast with method (ii), but absence of this part of the waveform in the first mentioned case is not considered important for us to compare the two proposed methods for determining the focal mechanism.In each case the best results were grid searched in several depths to obtain best fit between real and synthetic data.

Figure 2 .
Figure 2. (a) Focal mechanisms of the S0235b event obtained by inversion of only direct waves(Malytskyy, 2016) and (b) by direct inversion for the full moment tensor(Brinkman et al., 2021) for a source depth of 32 km.

Figure 3 .
Figure 3. (a) The waveforms (converted to displacements) of the event S1222a.(b and c) The durations of direct P-and Swaves at the station.Records were filtered in the frequency range between 0.2 and 9 Hz and are shown for N, E and Z components in green, blue and red lines respectively.

Figure 5 .
Figure 5. Waveform fits of body waves at 22 km depth using the modified TAYAK velocity model.Direct P-and S-waves are filtered between 3 and 12 s.The synthetic direct waves are shown in red and the observed ones in black lines respectively.

Figure 6 .
Figure 6.Solutions for the frequency range 0.028-0.072Hz.(a) The variability of focal mechanisms for different depths are shown for synthetic test for full MT inversion.(b) The variability of focal mechanisms for different depths are shown for synthetic test for deviatoric MT inversion.(c) The full MT solution-Results for grid search beneath epicenter.

Figure 7 .
Figure 7. Variation of nodal lines and P, T axes for 100 results of uncertainty tests.(a) our best low-frequency full-waveform MT solution (b) solution "C" from Table 3 (c) solution "C" from Table3when we consider deviatoric MT (d) solution for synthetic data for mechanism "D" for single station ELYSE (e) solution for synthetic data for mechanism "D" for four stations.

Figure 8 .
Figure 8.The normalized real (black) and synthetic (red) seismograms for north-south and east-west component showing a match for the part of surface waves in the frequency range 0.028-0.072Hz.

Table 1
Modification of the TAYAK Velocity Model Used in the Inversion of Waveforms and forCalculation of Seismic Tensor (This Model Differs From  Jacob et al., 2022)

Table 1
Maguire et al. (2023)to present the results only for the slightly modified model Kim et al. (2023) (Table2) used in the paperMaguire et al. (2023).After testing

Table 3
The Order of Solutions With Prescribed Focal Mechanism of Low-Frequency Full-Waveform Inversion for Frequency Range 0.028-0.072Hz MALYTSKYY ET AL.