Numerical modelling of impact seismic sources using the stress glut theory

SUMMARY

Meteorite impacts have proved to be a significant source of seismic signal on the Moon, and have now been recorded on Mars by InSight seismometers. Understanding how impacts produce seismic signal is key to the interpretation of this unique data, and to improve their identification in continuous seismic records. Here, we use the seismic Representation Theorem, and particularly the stress glut theory, to model the seismic motion resulting from impact cratering. The source is described by equivalent forces, some resulting from the impactor momentum transfer, and others from the stress glut, which represents the mechanical effect of plasticity and non linear processes in the source region. We condense these equivalent forces into a point-source with a time-varying single force and nine-component moment tensor. This analytical representation bridges the gap between the complex dynamics of crater formation, and the linear point-source representation classically used in seismology. Using the multiphysics modelling software HOSS, we develop a method to compute the stress glut of an impact, and the associated point-source from hypervelocity impact simulations. For a vertical and an oblique impact at 1000 m s⁻¹, we show that the moment tensor presents a significant deviatoric component. Hence, the source is not an ideal isotropic explosion contrary to previous assumptions, and draws closer to a double couple for the oblique impact. The contribution of the point force to the seismic signal appears negligible. We verify this model by comparing two signals: (1) HOSS is coupled to SPECFEM3D to propagate the near-source signal elastically to remote seismic stations; (2) the point-source model derived from the stress-glut theory is used to generate displacements at the same distance. The comparison shows that the point-source model is accurately simulating the low-frequency impact seismic waveform, and its seismic moment is in trend with Lunar and Martian impact data. High-frequencies discrepancies exist, which are partly related to finite-source effects, but might be further explained by the difference in mathematical framework between classical seismology and HOSS’ numerical modelling.

1 INTRODUCTION

As space exploration and instrumentation progress, seismic investigation of extra-terrestrial bodies has become possible and two space missions have been able to detect seismic waves generated by quakes and surface impacts. The first mission was the Apollo Lunar Surface Experiments (ALSEP), which deployed a network of four seismometers from 1969 to 1977 (Latham et al. 1969) on the Moon. The second was the InSight mission (Banerdt et al. 2020), which operated the Seismic Experiment for Interior Structure (SEIS), composed of a Short-Period (SP) and a Very Broad-Band (VBB) seismometer, from November 2018 to December 2022 on Mars (Lognonné et al. 2019, 2020; Giardini et al. 2020).

In agreement with pre-landing estimates (McGarr et al. 1969), ALSEP detected seismic signals from more than 1750 meteoroid impacts (Latham et al. 1970; Nakamura et al. 1982; Oberst & Nakamura 1987, 1991). Similarly on Mars, impacts were thought to be a potential source of detectable signals prior to the landing of InSight (Davis 1993; Teanby & Wookey 2011; Lognonné & Johnson 2015; Teanby 2015; Daubar et al. 2018). Despite the lack of detections during the nominal mission from 2018 to 2020 (Daubar et al. 2020; Miljkovic et al. 2021; Fernando et al. 2022), several impacts were recorded during the extended mission. The first ones, located at close distances from InSight (less than 300 km), generated six very-high frequency (VF) seismic events identified by the Marsquake Service (Clinton et al. 2018; Ceylan et al. 2022) and were notably accompanied by late low-frequency, dispersed wave trains. These dispersed arrivals, predicted in pre-landing studies (Garcia et al. 2017), correspond to acoustic waves travelling through the Martian night-time acoustic waveguide (Xu et al. 2022). Using their seismic arrival times and backazimuth (Garcia et al. 2022; Daubar et al. 2023), the location of these events was determined and associated with fresh craters imaged around InSight by the Mars Reconnaissance Orbiter (MRO) cameras (Malin et al. 2007). Two other very large craters, 150 ± 10 and 130 ± 12 m in diameter, at 3460 and 7455 km distance from InSight, were later detected by the MRO imager (Posiolova et al. 2022), and associated with seismic events recorded by SEIS on 24 December and on 18 September 2021. These two significant events produced broad-band seismic motion with a moment larger than 4.

Impact seismic observations on both the Moon and Mars allow for the first time a comparative study, and as such raise a number of questions. For example, can a single approach relate the long period amplitude of the seismic signals to the crater size measured on both celestial bodies? Can we predict seismic directivity effects for impacts with relatively oblique velocity such as the 18 September 2021 event (Posiolova et al. 2022)? In order to address these questions, we must integrate differences between the Martian or Lunar subsurface, and deconvolve the combined effects of wave propagation and wave generation on the observed signal in the time and frequency domain. Achieving this goal requires a model of the impact seismic source, accounting for the dynamic process during which the energy and momentum of the impactor is transferred into elastic energy within the target body. Such a model must quantify the magnitude of the impact seismic event, as well as the various timescales of the source.

Existing models of the impact seismic source rely on small-scale laboratory experiments, scaling laws developed using far-field impact seismic signals on the Moon, or analytical models based on explosion seismology. In preparation for the Apollo mission, McGarr et al. (1969) proposed to describe the impact seismic source by measuring the ratio between the seismic energy E_s within the target and the kinetic energy E_k of the impactor. This energy ratio quantifies the amount of impactor energy that is converted into seismic waves and is called the seismic efficiency, k_s. It has been used in numerous impact related studies: the seismic efficiency of impacts can indeed be estimated using recordings of missiles, surface explosions and lunar impacts, or via hydrocode modelling (Latham et al. 1970; Shishkin 2007; Güldemeister & Wünnemann 2017; Rajšić et al. 2021a).

More recently, in the frame of InSight pre-launch activities, scaling laws between the impact crater diameter and the seismic moment of the source, M₀, were proposed (Teanby & Wookey 2011; Teanby 2015; Wójcicka et al. 2020; see the review of Daubar et al. 2018). One advantage of the seismic moment over the seismic efficiency is that the former scales linearly with the peak displacement or the peak velocity of a seismogram, thus facilitating detectability analyses. It derives from a more complex mathematical object known as the moment tensor, which contains a representation of each force couples exerted on the target in every direction.

In the past, various models have been proposed to evaluate the seismic moment, either from modelling or data. Wójcicka et al. (2020) tested different methods to compute the seismic moment from hydrocode impact simulations. One method uses an analytical expression of the seismic moment of an explosion, obtained from the reduced displacement potential of compressional waves by Müller (1973):

where K and G represent the bulk modulus and shear modulus of the target material, and <D> represents the residual radial seismic displacement on a sphere of surface S surrounding the impact region. Another method by Walker (2003) uses the seismic impulse to derive a measure of the radial (in the cylindrical sense) component of the moment tensor:

with v_r the radial velocity at a radial distance r from the impact and t the time since the impact.

In both methods, the moment is computed directly from seismic amplitudes. Lognonné et al. (2009) on the other hand, estimated M₀ from a scaling of artificial impact seismic data recorded by the Apollo seismometers. More precisely, in Lognonné et al. (2009) and Gudkova et al. (2011), the seismic amplitude of lunar impacts, and hence the seismic moment, is found to be proportional to the vertical component of the seismic impulse, which is amplified by ejecta.

All the aforementioned models show limits in describing the impact seismic source. The seismic efficiency proposed by McGarr et al. (1969) varies by several orders of magnitude across different studies: Daubar et al. (2018) report values of k_s ranging from 10⁻⁶ to 10⁻¹ depending on the calculation method. The alternative approaches with the seismic moment or seismic impulse reported above are not entirely satisfactory either. Often, an analogy between impacts and shallow explosive sources is made. Thus, the impact seismic source is considered isotropic, ignoring the momentum transfer and the impactor directivity evidenced by the ejecta and surface expression of recent impacts on Mars (Posiolova et al. 2022). Although authors have also dedicated efforts in measuring the seismic impulse of impacts (McGarr et al. 1969; Walker 2003; Gudkova et al. 2015), or in developing more complex moment tensor (Lognonné et al. 1994; Gudkova et al. 2015), no study to date was able to reconcile the respective contribution of momentum exchange and explosive cratering on the impact seismic source. Furthermore, both the moment tensor and the seismic impulse are point-source models, which overlook finite-size effects at the source. These limitations justify a more physics-based approach to the development of impact source models.

One of the most exhaustive representation of seismic sources was established by Backus & Mulcahy (Backus & Mulcahy 1976a, b), initially to describe earthquakes. They introduced the notion of ‘stress glut’, which allows to quantify the plastic processes happening in the source region and translate them into a field of equivalent forces. Lognonné et al. (1994) and Gudkova et al. (2015) were the first studies to use the stress glut theory, with simplified hypotheses on source dynamics, in order to analytically derive a moment tensor representation for an impact. In this work, we propose to revisit this approach, by combining a rigorous application of Backus & Mulcahy’s theory with state-of-the art modelling of impact dynamics. We obtain a semi-analytical description of the impact seismic source.

Today, several validated and benchmarked hydrocodes have proven successful in resolving the complex dynamics of impact cratering (Pierazzo et al. 2008; Wünnemann et al. 2011; Collins et al. 2012; Güldemeister & Wünnemann 2017; Wójcicka et al. 2020; Rajšić et al. 2021a; Caldwell et al. 2021). Here, we introduce a new numerical method to compute the stress glut and seismic source of an impact based on the HOSS hydrocode (Knight et al. 2020). Our method is tested by comparing seismic motion recorded in the impact simulation with the seismic motion modelled using our semi-analytical seismic source.

The combination of a new analytical representation of the impact source with a physics-based numerical model allows for a detailed understanding of impact wave generation processes, with the ability to investigate their sensitivity to a wide range of parameters such as impactor velocity and mass, and target seismic velocity and strength. Moreover, this approach enables direct comparison of the impact source with other types of seismic sources, which may in the future allow better identification of impact-generated seismic signals on the solar system bodies.

This article consists of three sections, in which we present (i) an application of the seismic representation theorem to the impact problem, (ii) a numerical method to calculate the stress glut and associated seismic source terms, results and interpretation for a vertical and oblique impact simulations and (iii) tests conducted to assess the ability of this source model in representing impact seismic signals, both numerically and using comparisons with Lunar and Martian data.

2 THE REPRESENTATION THEOREM APPLIED TO IMPACTS

The goal of this section is to develop an analytical representation of the impact seismic source. In what follows, vectors are written as lower-case bold letters, while upper-case bolds letters are used to name tensors of order two or more. Operations on tensor and vector components follow Einstein summation conventions.

2.1 Representation of seismic sources

2.1.1 Elastic and non-linear equations of motion

The dynamics of seismic waves is developed in the framework of continuum and linearized mechanics. For a planetary body of volume V and surface Σ, seismic displacements are expressed by:

These three lines represent the elasto-dynamics equation, the boundary conditions and the constitutive relation for an elastic medium of density ρ⁰ and stiffness tensor |${\boldsymbol{C}}$|⁠. These equations relate displacement |${\boldsymbol{u}}({\boldsymbol{x}},t)$| and strains |${\boldsymbol{\varepsilon }}$| within the material, to linear-elastic stresses |${\boldsymbol{\Psi }}$|⁠, surface forces |${\boldsymbol{f^\Sigma }}$| and volume forces |${\boldsymbol{f^V}}$|⁠. They are thoroughly detailed in the literature (Aki & Richards 2002; Dahlen & Tromp 1998). However, they are only valid for a linear-elastic material which exhibits little changes in shape and no change in mass during its deformation. Such conditions are not fulfilled when dealing with large deformations and stresses caused by an hypervelocity impact. We now relax several assumptions made in eq. (3):

• The linearization of the left-hand side of the equation of motion is abandoned in order to account for the variations in the true local density ρ and the advection of momentum caused by the impact.

• The material of V is no longer considered ideally linear-elastic, that is, stresses within it are no longer dictated by Hooke’s law.

• Loss of mass and momentum through surface Σ is introduced to account for the reactive effect of impact ejecta.

The loss of mass through surface Σ is the major difference with respect to previous works, such as Lognonné et al. (1994), where the impact source (in this specific study, associated to the Shoemaker Levy 9 impact) was embedded in the volume. We develop the formalism in Appendix  A by making use of the Reynolds Transport Theorem, following Irschik & Holl (2004). The equation of motion in volume V remains the same as in fluid mechanism (e.g. Landau & Lifshitz 1987) and can therefore be either written as:

or:

with |$\text{d}/\text{d}t$| the material derivative, and where |${\boldsymbol{v}}$| is the velocity, |${\boldsymbol{h}}^V$| is the volumetric density of forces and |${\boldsymbol{S}}$| the non-linear stress tensor. The stress |${\boldsymbol{S}}$| is different from the ideal elastic stress |${\boldsymbol{\Psi }}$| of eq. (3) and does not a priori follow Hooke’s law of elasticity.

In distinction to the equation of motion, the continuity of stress at the mass-less boundary of Σ must integrate the transfer of momentum through this surface by ejecta, which leads to:

with |${\boldsymbol{v^\Sigma }}$| the velocity of the surface bounding the non-ejected material.

2.1.2 Equivalent forces of the impact and seismic wavefield representation

An accurate representation of the seismic source of an impact requires non-linear phenomena of Section 2.1.1 (eqs 4 and 6) to be accommodated into the elastic system of eq. (3). To do so, following the method of Backus & Mulcahy (1976a), non-linear effects are introduced in the form of equivalent volume and surface forces |${\boldsymbol{\gamma }}^V$| and |${\boldsymbol{\gamma }}^\Sigma$|⁠. The updated system is the following:

The equivalent forces are obtained by equating system (eq. 7) with the true boundary conditions in eq. (6) and the true equation of motion in eq. (4):

Eq. (8) introduces the tensor |${\boldsymbol{\Pi }}$|⁠, which is the difference between ideal, elastic stresses and true, non-linear stresses associated to the current deformation state. It is named the stress glut and was first discussed in the work of Backus & Mulcahy (1976a), who demonstrated its key role in explaining indigenous seismic sources such as earthquakes. The stress glut measures the deviation of true stress from the stress predicted by Hooke’s law, and is consequently related to the amount of plastic processes taking place in a seismic source region. The ability of a local thermoplastic stress change to generate elastic motion was originally proven in Eshelby’s famous inclusion problem (Eshelby 1957). As noted by several authors (Backus & Mulcahy 1976a; Madariaga 2015), the stress glut is related, in a dynamic sense, to the stress-free strain proposed in Eshelby’s static approach. We also note here the similarity of eq. (8) with the expressions of equivalent forces proposed by Takei & Kumazawa (1994) and Lognonné et al. (1994), who generalized Backus & Mulcahy’s source formulation by including a non linear mass advection term.

We will now definitively ignore gravity and will therefore assume that |$h_i^V$| and |$f_i^V$| are both null. Thanks to the seismic Representation Theorem, equivalent forces of eq. (8) can now be used to build the response of the media to the seismic source, as further detailed in Appendix  B. When the surface Σ is a free surface, the expression of the n^th component of displacements at time t and coordinates |${\boldsymbol{x}}$| in V is given by:

Hence, seismic motion derives from a field of equivalent forces spread over coordinates |${\boldsymbol{\xi }}$| in space and exerted at time τ. G_in is the Green’s function of the target medium, which propagates information elastically from the sources at |${\boldsymbol{\xi }}$| and τ to the receiver at |${\boldsymbol{x}}$| and t. The expression of the Green’s function depends on the elastic equations of motion: if gravity and rotation of a planet are accounted for, the Green’s function must be defined using the operators of gravito-elasticity and of Coriolis forces. We refer the reader to Appendix  D or for example Dahlen & Tromp (1998) for more details on the linearized gravito-elastic equations of motion. Despite this change in equations, the Representation Theorem of eq. (9) remains in essence the same (Dahlen & Tromp 1998, section 5.3).

Note that most studies making use of the seismic Representation Theorem assume that the mass enclosed by surface Σ is constant. Again, the formation of ejecta during an impact constitutes a mass and momentum loss, which in turn adds equivalent terms in the Representation Theorem of eq. (9). In this study, we do not account for the effects of variable mass and volume on the source, as we will find them to be negligible for the simulated impacts. However, the reader can find an exact version of Betti’s relation for a variable-mass system in Appendix  C, as proposed by Minster (1974) and reported by Archambeau & Scales (1989).

2.2 Point-source of an impact

2.2.1 The point-source approximation.

Given the equivalent forces and expression of displacements developed above, it is possible to further simplify the representation of the source. Following the point source approximation, we assume that the source is exerted at a point |${\boldsymbol{\xi ^{*}}}$| in space. This approximation is valid if the source is sufficiently small compared to the receiver’s distance and other typical spatial variations (⁠|$|{\boldsymbol{\xi }}-{\boldsymbol{\xi ^{*}}}| \lt \lt |{\boldsymbol{x}}|$|⁠). Assuming that the Green’s function varies smoothly and that its derivative is defined in the inner side of the free surface, its Taylor’s expansion of order 2 is conducted around |${\boldsymbol{\xi ^{*}}}$|⁠, before being reinserted in eq. (9). Details on this technique can easily be found in the literature (see e.g. Stump & Johnson 1977; Jost & Herrmann 1989; Julian et al. 1998b). Through this multipole expansion, displacements now become:

Or, expressed as a convolution (Stump 1985; Julian et al. 1998a):

Two key source parameters are revealed:

The above developments mean that, once integrated over the volume V^t of the source region, the force field |${\boldsymbol{\gamma }}^V$| and tractions |${\boldsymbol{\gamma }}^\Sigma$| reduce to the following approximations:

• A vector force of the form |$\longrightarrow$|⁠, also called a monopole. It corresponds to the term F_i(τ) in eq. (12). This term encompasses all momentum changes caused by the seismic source.

• Force-less couples of the form |$\longleftarrow {{\cdot }} \longrightarrow$|⁠, also called dipoles. They are contained in the moment tensor terms M_ij(τ) in eq. (12).

If an exact non-linear solution of the wavefield and stress field is modelled in the whole source volume, for example numerically, then the equivalent forces can be evaluated. We develop below the expression of the monopole and dipole in terms of these wavefields.

2.2.2 Expression of the monopole of the impact source

The expressions of equivalent volume and surface forces from eqs (8) and (7) are inserted into the definition of F_i(τ) (eq. 12). The Gauss–Green–Ostrogradsky Theorem allows us to simplify the expression of the net vector force exerted on Earth by the impact:

or:

This vector force is merely the variation of momentum in the instantaneous target volume following the impact, in the approximation of constant density. Note here that the mass leaving volume V will lead, by momentum conservation, to a variation of the target momentum. Therefore, even though F_i(τ) is an integral of momentum restricted to the target volume V, it includes by definition the impulse resulting from ejecta motion.

2.2.3 Expression of the dipole of the impact source

An expression of moment tensor components is obtained by replacing equivalent volume and surfaces forces in eq. (12) using their definition in eqs (8) and (7):

Integration by part is used to further simplify the second and third integrals. We get:

which yield after rearranging the stress terms:

eq. (12) can also be re-arranged differently if we recall from eq. (7) that |$\gamma ^V_i = \rho ^0\frac{\partial v_i}{\partial t} - \partial _j \Psi _{ij}$|⁠. This time, using an integration by part on ∂_jΨ_ij, surface integrals can be eliminated and we obtain:

This expression, which only requires volume integrals, is better suited for an evaluation with a numerical calculation. Like relation (14), it only requires the computation of velocity, elastic and non-linear stress fields within finite volume elements.

2.2.4 Remarks on these results

We have now gone from a source composed of distributed force fields, to a point source approximation composed of a monopole |${\boldsymbol{F}}$| and dipole |${\boldsymbol{M}}$|⁠. Further description of the source requires an evaluation of physical fields such as density, velocity, true stresses and ideal stresses within the material, which will be made possible with numerical modelling. Meanwhile, the expression of |${\boldsymbol{F}}$| and |${\boldsymbol{M}}$| already provide insights into the mechanism at play at the impact seismic source:

The vector force|${\boldsymbol{F}}$| (eq. 14) takes a rather elegant form, which is simply the integral of the linearized momentum within the instantaneous target volume. It hints at the fact that part of the impact seismic signal is due to momentum exchange between the impactor and ejecta. It also means that |${\boldsymbol{F}}(t)$| should be an impulse-shaped function, which decays with time as the impacted volume slowly relaxes to a new equilibrium. As such, it is not able to describe the permanent deformation of the surface, that is the crater, as this would require the application of a constant equivalent force on the surface. Such permanent deformations are linked to the creation of plasticity, and will therefore be accommodated by the stress glut present in the dipole term of the point source.

The typical duration of the monopole impulse is strongly dependant on the dynamics of the impactor and ejecta. The deceleration of the impactor occurs mainly in the first stage of crater formation, called the contact and compression stage. Its duration is typically the time needed for the impactor to burry itself, that is r_i/(v_isin (θ)), with r_i, v_i the radius and velocity of the impactor, and θ its incidence angle (Melosh 1989; Melosh & Ivanov 1999). This contact and compression time should thus be one of the dominant time scales of the monopole source.

The moment tensor: The expression of M_ij(t) in eq. (18) contains two integral terms. The first term is the first moment of momentum, which is similar in dimensions to the angular momentum (⁠|${\boldsymbol{L}}={\boldsymbol{r}}\times {\boldsymbol{p}}$| for a point in space with momentum |${\boldsymbol{p}}$|⁠). The second term is the ideal elastic stress, that is the stress computed from the non-linear strain field using an ideal elastic stiffness tensor. Therefore, M_ij(t) originates both from momentum exchange and from the non-linear elastic behaviour of material close to the source.

From eq. (18), and recalling that Ψ_ij = Π_ij + S_ij, we see that the time-evolution of M_ij(t) depends on the time evolution of v_i(t), S_ij(t) and Π_ij(t). We can infer some properties about the first integral of eq. (18) from the behaviour of v_i(t) in the far field. As the impactor momentum dissipates into the target volume, regions of non-zero velocity will concentrate in the seismic P and S waves. In the far-field, these waves decay proportionally to |$1/|{\boldsymbol{\xi }}_{P,S}-{\boldsymbol{\xi ^{*}}}|$|⁠, where |${\boldsymbol{\xi }}_{P,S}$| is the position of the P or S pulse with respect to the point-source (Aki & Richards 2002). Therefore, the first term of M_ij(t) should reach a constant value in the far field. Note also that this term is analogous to the expression of the radial seismic moment proposed by Walker (2003) (eq. 2).

The second term of M_ij(t) should also converge to a constant value, illustrating the presence of a residual plastic deformation at the source location. Similarly to Backus & Mulcahy (1976a), we can write that |$\lim _{t\rightarrow +\infty } \Psi _{ij}({\boldsymbol{\xi }},t) = \lim _{t\rightarrow +\infty } \Pi _{ij}({\boldsymbol{\xi }},t)= C_{ijkl}\, \varepsilon ^{P}({\boldsymbol{\xi }})$|⁠, with |$\varepsilon ^{P}({\boldsymbol{\xi }})$| being the residual plastic strain field of the source. This non-zero value of M_ij after the seismic event mean that couples of forces are being permanently exerted within V in order to maintain its new permanently deformed shape. We point out to the reader that there exists an interesting relationship between the integral of the stress glut and the expression of the seismic moment of an explosion due to Müller (1973) presented in section 1 (eq. 1). In a discussion in 2005, Richards & Kim pointed out the debate surrounding two different expressions of the moment of an explosion. The expression of Müller (1973), |$M_{0,1} = (K+\frac{4G}{3})\, \delta V$|⁠, with δV = S < D >, K and G the bulk and shear modulus of the medium, competed with another expression: M_{0, 2} = KΔV. The difference between the volume changes δV and ΔV involved in the two expressions was yet unclear. Richards & Kim (2005) resolved this debate by stating that these two moment definitions are in fact equivalent, but that the meaning of each volume is different. Precisely, the volume ΔV in M_{0, 2} corresponds to the permanent volume change experienced by the strongly loaded material of the source. Again, we find a parallel between the second integral of M_ij(t) (eq. 18) and Müller’s equivalent M_{0, 2}. Indeed, if we consider a purely isotropic, compressive source and assume the plastic strain to be related to the volume change by |$\varepsilon ^P_v = \frac{\Delta V}{V_{source}}$|⁠, the second integral of M_ij(t) becomes:

Thanks to Richards & Kim (2005), we can thus conclude that the last term of the definition of the seismic moment in eq. (18) is equivalent to the expression of Müller (1973) used in Wójcicka et al. (2020).

Finally, we indicate a possible further simplification of the Moment Tensor expression. We recall from an integration by part that |$\int _V S_{ij} \;\mathrm{d}V=\int _\Sigma S_{ij} (\xi _j-\xi ^{*}_j) \;\mathrm{d}\Sigma -\int _V \partial _jS_{ij} (\xi _j-\xi ^{*}_j) \;\mathrm{d}V$|⁠. In the far-field, where density changes and Reynolds inertial effects can be neglected, we have: |$\partial _j S_{ij}(t)(\xi _j-\xi ^{*}_j)\sim \rho _0 \frac{\partial v_i}{\partial t}(\xi _j-\xi ^{*}_j)$| from eq. (4). This means that in the far field, the ‘angular momentum’ integral and the ‘true stress’ integral of eq. (18) should amount to:

and that the moment should be reduced to:

This approximation and the respective amplitudes of the two remaining source terms will be further discussed in Section 3.5.

As a final note, we point out that, while this second integral of eq. (18) yields a symmetric tensor by definition of the Cauchy stress tensor, this is not necessarily true for the first integral. In particular, if the impact problem is not cylindrically symmetric, such as during an oblique impact, the first integral conveys the change of angular momentum imparted to the celestial body by the bolide.

2.3 Towards a more detailed source model

All the above development aims at obtaining a simple model of the impact seismic source, in the form of a point of origin for a force and six couples of forces. While the point-source approximation has proved to be successful in a large variety of seismic studies (inversion of sources, estimations of signal amplitude and magnitude, etc.), it does not provide a full description of the seismic source physics. One of its biggest limitation is that it overlooks the effects of the source kinematics and of its finite dimensions on the observed seismic signal. In particular, the geometry of a source has an influence on the seismic signal. Indeed, in eq. (9), the source energy radiates from multiple points |${\boldsymbol{\xi }}$| in space, like an antenna. Hence, as the shock front progresses, the delayed radiation of two distinct points causes interferences in the generated seismic signal and alters its frequency content. This phenomenon has been at the focus of many earthquake studies, for example by Savage (1966), Haskell (1969) or Madariaga (1976) for elliptical, rectangular or circular fault ruptures. The geometry of an impact source is likely to be very different from an essentially 2-D fault surface.

The development of analytical models for an extended impact source and its kinematics is beyond the scope of this study. However, we will show in Section 4.2.2 that we can reproduce an extended source numerically using the stress glut and velocity fields described above. The effects of the finite source dimension on the signal spectrum can then be assessed.

3 NUMERICAL COMPUTATION OF THE SEISMIC SOURCE

The above analytical developments quickly reach a limit in their ability to describe the source. Further understanding of the impact seismic source requires us to compute several fields (velocity, stress and displacements) evolving in a strongly non-linear regime. Numerical modelling is adequate to our development as it provides all the needed fields in a discretized space. In this section, we present a numerical method to quantify the stress glut field associated to impacts, and to estimate the source terms developed in Section 2.2. We show an application of this method to two impact scenarios with different incidence angles and we analyse their source mechanism.

3.1 Numerical model and methods

3.1.1 The HOSS software

We develop algorithms and numerical methods to compute the stress glut field using the Hybrid Optimization Software Suite (HOSS). HOSS is based on the Finite-Discrete Element Method (FDEM; Munjiza 2004; Knight et al. 2020), which combines the Finite Element method for the description of continuum, with the Discrete Element (DE) Method to simulate fractures, fragmentation and interparticle interactions. HOSS supports large-scale parallelization and is equipped with a variety of tunable material models for an accurate description of metals and geomaterials (Lei et al. 2014). As such, it has been applied to simulate various impact problems (Froment et al. 2020; Caldwell et al. 2021), and to model dynamic fracture processes at play during earthquake ruptures and their effects on the frequency content of near-field seismic radiation (Okubo et al. 2019).

3.1.2 Material models

HOSS’ simulation domain is a user-defined mesh of tetrahedral elements. In the framework of FDEM, each tetrahedron constitutes a Finite Element (FE), but the bounding surfaces of each FE are also discretized by a set of DEs. The interaction of adjacent elements is calculated through these DEs. Hence, neighbouring elements can be bounded and form a continuous volume, or they can behave independently and interact via contact and friction. Consequently, HOSS material models are of two kinds, in order to define the response of both FEs and DEs to deformation.

The response of FEs is typically separated into a volumetric and a deviatoric material model. The volumetric model specifies the evolution of pressure as a function of the volumetric strain of an element, while the deviatoric model defines the six remaining deviatoric components of stress. In this work, we develop a stress glut calculation method based on the volumetric and deviatoric models used by Froment et al. (2020), described below. Indeed, as shown by the authors, these models successfully reproduce the response of a porous and granular material, similar to the regolith present on the surface of Mars and the Moon, to impacts in laboratory conditions. Each response models handles elasticity and plasticity in its own way, which leads to two separate methods of computation of the stress glut. FEs may also be subject to an additional damping stress, which is added in order to reproduce the visco-elastic behaviour of realistic geomaterials. This stress term is purely inelastic, and can be simply added to the measure of the stress glut. In this study, the contribution of damping stress on the stress glut was determined to be negligible, and it will hence be ignored in the following sections.

In contrast with FEs, the response of DEs depends on the current conformation of simulation elements. For two initially bounded, adjacent elements, it takes the form of a strain softening curve: elements with opposite motion gradually dissipate energy into a fracturing process, until their maximal fracture aperture is reached (see e.g. Rougier et al. 2014). Beyond this point, the two elements become independent, and will interact together via a contact detection and friction algorithm associated to a penalty that prevents interlocking. In the following approach, for simplicity, we do not consider the effect of fracturing on the impact seismic source. To ensure that the fracture phenomenon does not affect the estimate of the source, we define the material strain softening curve as having a small fracture energy compared to the energies of other dynamic processes. On the other hand, friction between particles is an important mechanism at play in the response of geomaterials and constitutes an additional force in our simulations. Although this interelement friction has not been included in the equations of Section 2.2 and in the following development, we discuss its possible effects on the seismic source in Section 5.

3.2 Computation of the stress-glut

We present numerical methods adopted to compute the stress glut. Numerical tests of these algorithms for a single simulation element are found in the Supplementary Information (Section SI-2).

3.2.1 The pressure glut

The volumetric deformation of the material is represented using the SocCrush model, based on an algorithm by Schatz (1974). In this model, the deformation of an element is separated into three domains. At low pressure, the material follows Hooke’s law: the pressure P and volumetric strain ε_v are related by |$P=K_{el} \, \varepsilon _v$|⁠, with K_el being the bulk modulus of the elastic material. In this regime, the modelled pressure and the ideal, elastic pressure are equal: by definition, the volumetric term of the stress-glut (or ‘pressure glut’) is thus zero (see Fig. 1). Above a certain limit stress, the pressure–strain relationship departs from this ideal linear behaviour and the material starts accumulating plasticity. Froment et al. (2020) as well as laboratory studies (Yamamuro et al. 1996; Luo et al. 2011; Housen et al. 2018) showed that an exponential curve provides a good description of the pore-crushing and compaction phenomenon that occur in a regolith-like material. In this regime, we have:

Once pressure in the material starts decreasing, it follows a new linear-elastic unloading path back to P^True = 0 (see Fig. 1). Each modelled element that has entered the pore-crush regime retains a final volumetric strain ε^P. Following eq. (22), this means that in the last stages of deformation, its pressure-glut converges to |$\lim _{t\rightarrow \infty }\Pi ^V(t)=K_{el}\, \varepsilon ^P$|⁠, as mentioned in the remarks of section 2.2.4.

3.2.2 The deviatoric stress glut

We now consider deviatoric stresses using the algorithm described in Lei et al. (2020). The deviatoric stresses of an element are related mathematically to seven different modes of deformation: one representing pure volumetric strain, which is handled by the method of section 3.2.1 above, and six modes of pure shear deformation. The six shear modes refer to six angles of deformation, written φ_{1 − 6}, which are measured within a reference material element. At each new simulation time step t_{N + 1}, the algorithm predicts the angles of the deformed material element, φ_{1 − 6}(t_{N + 1}). It verifies whether the deformation associated to φ_{1 − 6}(t_{N + 1}) have brought the element beyond its yield surface. If it is the case, a correction is applied to φ_{1 − 6}(t_{N + 1}) using a return mapping method, so that corrected angles |$\varphi ^E_{1-6}(t_{N+1})$| keep stresses within the yield surface. This correction is equivalent to a measure of plasticity: total angles φ_{1 − 6}(t_{N + 1}) are separated into ‘elastic angles’ |$\varphi ^E_{1-6}$|⁠, and their plastic correction or ‘plastic angles’ |$\varphi ^P_{1-6}$|⁠, as in:

After the correction step, a linear relationship transforms angles |$\varphi _{1-6}^E$| into the deviatoric Cauchy stress |$\boldsymbol{S}^D$| of the element. This action can be summarized by a linear operator F_ij [detailed steps can be found in Lei et al. (2020)]:

By analogy with the volumetric stress glut, the deviatoric stress-glut |$\boldsymbol{\Pi }^D$| will thus be defined as the difference between ideal stresses ‘if all deformations associated to φ_{1 − 6} were elastic’ and true stresses:

Because F_ij is linear, the deviatoric stress glut is simply obtained by computing the stresses associated to plastic angles:

In the same way that a residual volumetric strain ε_P is created in the plastic regime of the volumetric model in Section 3.2.1, the element might also accumulate a final plastic angle |$\varphi ^P_{1-6}$| in each mode of shear deformation.

Note that this measure of plasticity requires us to compute the deformation of a reference element, with its own fixed reference axes. The computation steps hidden within the function F_ij, aim at mapping the deformation of this reference element onto the global simulation space, taking into account the position of each simulation element. Contrary to the volumetric stress-glut, which is a measure of the trace of the stress-glut tensor, and is consequently independent of any change of reference frame, the deviatoric stress-glut is sensitive to geometrical changes brought for instance by the rotation of a simulation element (see section 2 of the Supporting Information for an example). This sensitivity of deviatoric stress-glut to rotation can be problematic. Indeed, even if inelastic deformation has stopped and plastic angles |$\varphi ^P_{1-6}$| are constant, the components of the stress-glut tensor will change over time as long as the position of the reference frame is changing with respect to the global simulation frame. Long-term variation of the stress-glut may appear, even though the exchange of forces with simulation elements may have ended. A method to quantify the effect of rotation on the source is proposed in the next section.

3.3 Final numerical representation of the source

Computing the stress glut is one piece of the solution to the problem of representing the impact seismic source. Here, we state the protocol established for the full determination of the seismic point-source:

• A numerical simulation of an impact is run. The simulation domain must be large enough to encompass the entirety of the inelastic source region. The simulation is stopped when the displacement recorded on sensors outside the inelastic source region has stabilized.

• Images of the full simulation domain and its various fields (⁠|${\boldsymbol{\Pi }}^V$|⁠, |${\boldsymbol{\Pi }}^D$|⁠, |${\boldsymbol{v}}$|⁠...) are saved at regular time intervals of typically 0.1 ms for a metre-size crater. A shorter interval (∼1e − 6 s) is used in the first few milliseconds of the impact to correctly capture the very fast exchange of momentum between the impactor and target.

• The volume integrals of eqs (14) and (18) are computed in a discrete way by summing variables over each simulation element. Because simulations are often run with a restricted domain due to computational costs (for example: a 45° or 180° cylindrical slice centred on the impact point), components of the physical fields or of the source have to be rotated or mirrored to recover the full 360° volume. In the source computation, a distinction between the ejecta (i.e. material not in contact anymore with the planet) and the target (i.e. the planet) is made: the target volume is considered to be composed of every element whose centroid coordinate z_e is below the ground surface z = 0. For example, the force component of the source is given by:

  $$\begin{eqnarray} F_i(\tau ) = \frac{1}{t_{j+1}-t_j} \left[ \left( \sum _{e=0}^{N_{E,j+1}} \rho ^0[e] v_i[e] V^0[e] \right) _{j+1} - \left(\sum _{e=0}^{N_{E,j}} \rho ^0[e] v_i[e] V^0[e] \right)_{j} \right], \end{eqnarray}$$

  (27)

  where j and j + 1 designate two successive simulation images with time t_j and t_{j + 1}, respectively, N_{E, j} is the number of element in the target domain at time t_j, V⁰[e] (resp. ρ⁰[e]) is the initial, undeformed volume (resp. density) of one of these elements, and v_i[e] the velocity of its centroid.

• To measure the influence of rotation on the integrated moment tensor, an alternative version of it is computed using the ‘corotated Cauchy stress’ |${\boldsymbol{\bar{\Psi }}}$| as a measure of stress and stress glut. This alternative stress measure infers the finite rotation of an element from the polar decomposition of its deformation gradient |${\boldsymbol{F}}[e]$| (Hoger & Carlson 1984), such that:

  $$\begin{eqnarray} {\boldsymbol{F}}[e] = {\boldsymbol{R}}[e] {\boldsymbol{U}}[e] = {\boldsymbol{V}}[e] {\boldsymbol{R}}[e], \end{eqnarray}$$

  (28)
  where |${\boldsymbol{R}}[e]$| is the rotation tensor. The corotated formulation of stresses is then:

  $$\begin{eqnarray} {\boldsymbol{\bar{\Psi }}}[e] = {\boldsymbol{R}}^T[e] {\boldsymbol{\Psi }}[e] {\boldsymbol{R}}[e]. \end{eqnarray}$$

  (29)

These post-processing steps provide three measures of the point-source of the impact: one force component F_i(τ), and two time-varying moment tensors |$\bar{M}_{ij}({\boldsymbol{\xi }}^{*}, \tau )$| and |$M_{ij}({\boldsymbol{\xi }}^{*}, \tau )$|⁠, with and without rotation correction.

3.4 Application to two impact scenarios

The approach described above is tested for two impact scenarios. In both cases, a spherical basaltic bolide with a mass of 12 kg and a velocity of 1000 m s⁻¹ impacts a cohesive surface similar to Martian regolith. In scenario (A), the impactor has a vertical incidence angle, compared to 45° in scenario (B). We describe here the characteristics of the simulation geometry and material models.

3.4.1 HOSS simulation geometry design

The two simulations are conducted in 3-D. Due to the cylindrical symmetry of the vertical impact problem, the target geometry of scenario (A) consists of a cylindrical sector of 45° angle. The target geometry for the oblique impact of scenario (B) is a half-cylinder with 180° angle to accommodate its plane symmetry. The impactor is a section of a sphere with a r_i = 10 cm radius and with a 45° or 180° aperture, placed at the centre of rotation and at the top of the target cylinder.

To ensure the correct computation of the stress glut, the target domain must be larger than the region of non-linearities close to the crater. Consequently, an initial simulation is performed with a rough mesh to estimate the size of the inelastic source region, defined as the radius R_s beyond which the stress glut vanishes. The final mesh ensures that the radius and depth of the target cylinder exceed R_s.

Here, the resulting target domain extends to 17 m. This is about a hundred times larger than the impactor, which presents a computational challenge. In an attempt to limit the total number of simulation elements, different mesh resolutions are tested. The mesh test and parameters of the final mesh for scenarios (A) and (B) are presented in the SI (section SI-3). Shock wave amplitudes in the final meshes are within 85 per cent of the ones obtained with higher-resolution meshes. These meshes provides a good compromise between accuracy, element number and load balance between parallelization domains.

3.4.2 HOSS material model design

The impactor and target materials are simulated with a model adapted from Wójcicka et al. (2020). The target model was previously implemented in HOSS in order to simulate the impact of Perseverance’s entry balance masses on Mars (Fernando et al. 2021), and it aims at representing the response of the upper ten meters of the martian surface regolith. The target deviatoric response is modelled using a Lundborg pressure-dependant strength model (Lundborg 1968), and the volumetric response with an elasto-plastic model for porous materials adapted from Froment et al. (2020), with parameters from Wójcicka et al. (2020). The impactor is made of solid basalt and simulated with a Tillotson equation of state (Tillotson 1962) and another Lundborg strength model. Parameter values and details for both materials can be found in Table 1 of the SI. Other target models could have been used for this work, but this particular model is chosen due to its high elastic velocities (v_p = 1090 m s⁻¹ and v_s = 583 m s⁻¹), which helps with simulation convergence. We leave the investigation of impact sources in other materials to a future study.

3.5 Results: source of a vertical and an oblique impact

3.5.1 Impact dynamics and stress glut

Images of the simulated stress glut fields |$\Pi ^D_{zz}$| and Π^V are shown on Fig. 2 for the vertical impact and on Fig. 3 for the oblique impact, at three different times. For both simulations, the stress glut field shows a transient stage until around 15 ms, during which the inelastic region is growing. After 15 ms, the field appears to have stabilized. With a radius of around 10 m compared to only 5 m, the region of non-zero stress glut is larger for deviatoric processes (fourth line, |$\Pi ^D_{zz}$|⁠) than for volumetric processes (third line, Π^V). This is not surprising according to the definition of material response in Section 3.1.2: the computation of the volumetric and deviatoric stress glut rely on two different material models. Volumetric stress glut appears when the pressure in the material exceeds its crushing strength P_el, while the presence of deviatoric stress glut is determined by the material’s yield surface, i.e. its shear strength and internal friction, which leads to two different spatial dimensions. Note also that the deviatoric stress glut field grows at the same speed as the region of non-zero rotational velocity (⁠|${\boldsymbol{\nabla }}\times {\boldsymbol{v}}$|⁠): this confirms that the deviatoric stress glut is inherently related to shear processes and thus to the shear (S) wave propagation. On the other hand, the volumetric stress glut follows the pressure (P) wave.

The craters formed in both simulations are shown 116 ms after the impact on Fig. 4. The vertical impactor results in a crater of about 70 cm depth, while the oblique impactor crater reaches only about 65 cm in depth. Moreover, while the vertical impactor remains at the bottom of the crater, the oblique one is subject to a rebound, and is shown to have escaped the crater at 116 ms (see Fig. 4b). This is an interesting behaviour, as the impactor rebound could enhance the impulse transferred to the target in the vertical direction. Although the 116 ms of simulation are not sufficient to capture the final crater dimension with zero gravity, we can estimate lower bounds for the final crater diameter of 1.9 m for the vertical scenario and 2 m for the oblique scenario.

3.5.2 Source parameters

The stress, stress glut and velocity fields shown in Figs 2 and 3 are integrated to obtain the point-source components, following the method of Section 3.3. In order to generalize these results, we scale the obtained point source terms in time and in amplitude. The time is scaled to the duration of the contact and compression stage for the vertical impact:

with impactor radius r_i and total velocity v_i giving τ_N = 1e − 4 s. The point source force component is scaled to the measure:

corresponding to the hypothesis that the total momentum P_i is delivered to the target within a time τ_N. We have here F_N = 1.2e8 N. Finally, the point source moments are dimensionally scaled based on the equivalence between the seismic response to a point force and to a moment tensor (see e.g. Aki & Richards 2002). This gives:

a scaling also proposed by Daubar et al. (2018) for a known impactor and material velocity. Here, M_N = 1.3e7 Nm.

Fig. 5 shows the components of the force |${\boldsymbol{F}}(\tau )$| representing the monopole of the seismic source, for both simulations. The force takes the form of a pulse with a duration of about ∼5 τ_N, which confirms that most of the impactor momentum is delivered to the target during the contact and compression stage as noted in Section 3.5. For the vertical impact, the vertical force amounts to 0.28 F_N (∼3.3 × 10⁷ N), compared to 0.19 F_N (∼2.3 × 10⁷ N) for the oblique impact. In both scenarios, the maximum force is less than F_N, which suggests that the impactor penetrates the target for about 3.5 times its radius before delivering its total momentum P_i during the contact and compression stage. We also note that for the oblique impact, the amplitude of the force in the |$\vec{z}$| direction is well predicted by projecting the force of the vertical impact at 45°, as shown by the dashed black curve in Fig. 5(a). Indeed, 0.19 ≈ sin (45°) × 0.28. However, the horizontal force reaches a maximum of only 0.13 F_N (1.5 × 10⁷ N), suggesting that the transfer of momentum is less efficient in the |$\vec{x}$| direction. A possible explanation is the lingering motion of the oblique impactor observed in Fig. 4(b) along |$\vec{x}$|⁠, which reveals that the impactor keeps part of its horizontal momentum and does not transmit it to the target. It could also be due to a difference in the generation of ejecta between both cases.

Fig. 6 displays the components of the moment tensor extracted from HOSS, either with or without correcting for element rotation. We adopt the traditional sign convention of seismology, were a positive moment signifies a compression of the material surrounding the source. The oblique and vertical moment functions present similar amplitudes reaching up to 8 M_N (∼1e8 Nm). However, some of the components differ strongly. In particular, the oblique simulation is characterized by the presence of a non-zero M_xz component. Without correcting for rotation (Fig. 6a), the M_xx, M_yy and M_zz components have rather similar shapes in both the vertical and oblique simulation. However, the correction of rotation in Fig. 6(b) has significant effects: the average amplitude of |$\bar{M}_{zz}$| is increased in the vertical simulation, while it is decreased in the oblique simulation.

As the target reaches an equilibrium, that is a state of constant deformation, we would expect the moment to converge to a constant value. Here, we note that, even after correcting for rotation in plot 6(b), the moment source time function still shows residual variations instead of converging. This long-term evolution of moments is an important aspect of the seismic source. This feature could be due to the ejecta elements still leaving the crater at large time scales, carrying with them part of the stress glut and thus artificially altering the source. Indeed, Figs 2(c) and 3(c) indicate that some elements close to the crater still have strong positive vertical velocities after 116 ms. In the absence of gravity, they continue to escape the target domain, while they would eventually settle down if it was taken into account.

3.5.3 Analysis of the source mechanism

We investigate the source mechanism further by looking into different decompositions of the moment tensor. The simplest possible decomposition consists in separating the tensor into its volumetric and deviatoric parts. This is also a meaningful operation, as it echoes the ways in which moments are integrated from a volumetric and a deviatoric stress glut in our method, themselves calculated using two different material models. This decomposition unveils three classes of source time functions, represented in Fig. 7. The volumetric moment M^V, defined as |$M^V=\frac{1}{3}(M_{xx}+M_{yy}+M_{zz})$| reaches a quasi-constant value in a short time of about 50 τ_N (Fig. 7a). The diagonal terms of the deviatoric moment tensor, defined as |$M^D_{ij} = M_{ij}- \delta _{ij}M^V$|⁠, reaches a maximum around 200 τ_N, then decrease over long time scales until they reverse sign (Fig. 7b). Finally, the non diagonal term |$M^D_{zx}$| presents a maximum around 300 τ_N, and a more tempered decrease (Fig. 7c). The first two source time functions are similar for both impact scenarios, but |$M_{zx}({\boldsymbol{\xi ^{*}}}, \tau )$| is unique to oblique impacts. The timescales of 50–300 τ_N (5–30 ms) are likely associated to the formation of the strongly damaged region below the crater, illustrated in Figs 2 and 3. The smaller extent Π^V supports a shorter formation time, and thus a shorter timescale for volumetric moments compared to deviatoric ones.

We present the maximum amplitude of the volumetric and deviatoric moments in Table 1. We note that contrary to the force component, the moment components of the oblique impacts are generally larger than the values predicted by multiplying the vertical impact source with sin (45°), and larger than M_N. Hence, no simple relationship between moments and impactor impulse appears. More importantly, the maximum amplitude of deviatoric components appear comparable to the volumetric component. This is an indication that the impact source mechanism is more complex than the isotropic explosion chosen by previous studies (Shishkin 2007; Teanby 2015; Wójcicka et al. 2020).

The decomposition of a moment tensor into basic source mechanisms, is non unique. We have presented above a decomposition into |$M^V_{ij}$| and |$M^D_{ij}$|⁠. For the vertical impact, this is equivalent to a decomposition into an isotropic (ISO) mechanism equal to M^V, and a Compensated Linear Vector Dipole (CLVD) mechanism equal to |$M^D_{ij}$|⁠, thus having a principle axis aligned with the vertical direction. However, the full moment tensor |$M_{ij} = M^V_{ij}+ M^D_{ij}$| of the vertical impact could also be viewed as a crack (C) opening in the vertical direction.

On the other hand, due to the presence of a non-diagonal component and the fact that M_xx ≠ M_yy, the moment tensor of an oblique impact cannot be represented by a simple crack or ISO + CLVD mechanism. Other, more general decompositions involving Double Couple (DC) mechanisms should instead be used. For earthquake source, common choices are ISO + CLVD + DC, ISO + DC + DC, or a ‘crack + double couple’ (CDC) mechanism (see e.g. Dahlen & Tromp 1998; Minson et al. 2007; Tape & Tape 2013; Shearer 2019). These mechanisms all define a fault plane along which the material shears. Contrary to the vertical impact case, the vertical direction is no longer a natural direction for these cracks, CLVD or DC mechanisms.

The visualization of these possible mechanism can be made clearer by placing the moment tensor |$M_{ij}({\boldsymbol{\xi ^{*}}},\tau )$| on a Lune diagram. This type of diagram was proposed by Tape & Tape (2012) to facilitate discrimination of seismic sources on Earth: it maps every possible moment tensor onto a 2-D space, giving them two coordinates (δ, γ). These coordinates are calculated from the tensor eigenvalues, thus any moment tensor has a unique position of the Lune. The proximity of the moment tensor to specific kind of sources (Isotropic, CLVD, DC) can this way be measured. On Fig. 8, we applied this representation to the vertical and oblique impact, using the original and corotated moment tensor.

For the vertical impact (Figs 8 and b), both expressions of the moment tensor lie on the top part of the Lune diagram, in a region containing purely extensive sources. Because the moment tensor of a vertical impact has two identical eigenvalues, it stays on the borders of the diagrams, close to a pure explosion. On the other hand, the moment tensor of the oblique impact possesses three distinct eigenvalues: as such, it is placed in the inner part of the diagram (Figs 8c and d). It is not purely extensive: the effect of the directivity of the impact is visible from the white region of the beachball, where the source is compressed. The arc-length between two points in the Lune diagram measures the difference between two moment tensors in matrix space, if their orientation and norm are ignored (Tape & Tape 2019). The source of the oblique impact lies at equal distance from the Isotropic, CLVD and DC source in the Lune diagram, therefore, it is equally close to each of these physical source processes. The diagram also highlights the time variability of the source: although the vertical impact appears similar to an explosion at some times during the cratering process, its moment tensor also travels through lower latitudes over time, showing the complexity brought by the deviatoric moments to the source mechanism. The rotation correction diminishes this time variability, resulting in a simpler source mechanism.

3.5.4 A possible simplification of the moment tensor

The moments represented in Fig. 6 are computed by summing the two integrals of eq. (18). In a final analysis of the seismic point source, we would like to test the hypothesis of Section 2.2.4 on a possible simplification of the expression of |$M_{ij}({\boldsymbol{\xi ^{*}}},\tau )$|⁠. To this aim, we can plot integrals of eq. (18) individually to determine their respective contribution. This was done in particular with the ‘angular momentum’ term |$\int _V\frac{\partial \rho ^0v_i}{\partial \tau }(\xi _j-\xi _j^{*})\;\mathrm{d}V$| and the ‘true stress’ term ∫_VS_ij dV. As seen on Fig. 9, taken individually, the amplitude of these two terms is of about 1 × M_N (∼10⁷ Nm), thus comparable to the total moment components (Fig. 6). However, the bottom plot of Fig. 9 shows that their sum becomes negligible a few milliseconds after the impact. We recall the exact expression of their sum:

The fact that the term on the left-hand side of eq. (33) is negligible might be associated to different explanations based on the term of the right: either the three terms compensate each other, or they are individually negligible. For example, as mentioned in Section 2.2.4, it is likely that the terms |$\left[(\rho ^0-\rho )v_i\right](\xi _j - \xi _j^{*})$| and (ρv_iv_j) are negligible except in the very first milliseconds of the impact, when the velocity of target particles are comparable to those of the impactor, leaving only the stress term S_ij. As mentioned in Section 1, the ‘angular momentum’ term is analogous to the definition of the radial moment M_rr proposed by Walker (2003) and used in Wójcicka et al. (2020) (eq. 2). Given our observation that it is deemed to be compensated by the Cauchy stress term, the question of whether eq. (2) is really representative of the source moment deserves to be debated.

In the same way, we also measured the asymmetry of the moment tensor. The asymmetric part of the moment tensor was found to be negligible compared to the symmetric part after the first millisecond of the impact: indeed, it has the same timescale as the transfer of momentum. It was therefore assumed to be zero in the rest of this work.

4 TESTS OF THE SEISMIC SOURCE MODEL

The stress glut calculation method presented in Section 3 provides an exhaustive description of impact source mechanism. In this section, we propose to test this seismic source model in different ways. We will first evaluate to which extent numerical impact seismic signals from HOSS are reproduced by our simplified point source model. This is a form of validation of the point source model against a numerical experiment. In addition, we will compare the quantitative properties of our signals, such as the seismic moment and corner frequency, to available Martian and Lunar data.

4.1 A testing approach using software coupling

4.1.1 Principle of the method

Our source model was calculated based on numerical impact simulations solving the non-linear equations of motion, with realistic material models. As such, the displacement and velocity waveforms produced by HOSS represent a ‘true’ solution to the impact wave-generation problem. We now wish to evaluate how much of this numerical impact seismic signal is explained by our application of the representation theorem and our calculation of the stress glut.

The seismic signal corresponding to the point source model is the solution of eq. (11). We propose to compare the simplified point source signals with HOSS’ signals, which are numerical solution to eqs (4) and (6). However, such approach is only meaningful if both solutions are linear at the point of comparison. To ensure this, we need to calculate point source and HOSS waveforms at a distance larger than R_s, the radius at which stress glut becomes negligible (see Section 3.4.1).

Solutions to eq. (11) are easily obtained at any distance by convolving our point source parameters with a Green’s function solution. Such Green’s functions can be calculated in any elastic media and geometry by modern wave propagation codes. Here, we use the SPECFEM3D code, based on the Spectral Element Method (Komatitsch & Vilotte 1998). However, computational costs limit the extent of HOSS simulations to the very near field. In order to enable the comparison of point source and numerical waveforms at larger distances, we prolongate the HOSS solution by coupling its wavefield to the SPECFEM3D code, following the method of Larmat et al. (2015). The HOSS-SPECFEM3D coupling method relies on passing displacement time-series extracted from HOSS to a spectral element node with the same location in SPECFEM3D. This node then acts as a source within the SPECFEM3D simulation.

In order to transfer HOSS wavefield with a high accuracy, the coupling nodes must form a dense network of points on a surface surrounding the source. In addition, HOSS input time-series must be linear, meaning that the coupling interface must be outside of the inelastic source region of radius R_s. The duration of the SPECFEM3D simulation is limited by the duration of the provided coupling time-series. Moreover, the maximum signal frequency transferable to SPECFEM3D corresponds to |$f_{max}\approx \frac{\beta }{\lambda }$|⁠, with β the elastic S-wave velocity and λ the size of a spectral element. With such coupling method, SPECFEM3D may propagate HOSS wavefield at distances 10s to 100s of times larger than HOSS initial domain.

In conclusion, our testing method is summarized by Fig. 10. On one hand, HOSS’ wavefield is recorded at a dense network of points on a box surrounding, but outside of the inelastic source region. This wavefield is passed to SPECFEM3D and propagated up to 512 m distance. These Coupled signals constitute our numerical solution to the impact wave-generation problem. On the other hand, the point source parameters of Section 3, calculated from the non-linear source region using the representation theorem, serve as a source for a different type of SPECFEM3D simulation. These Point source signals constitute our approximation of the impact seismic source. Both signals are recorded on similar stations and compared.

4.1.2 Test simulations design

The point-source and coupled signal comparison is performed for both scenarios (A) and (B) described in Section 3.4. In each case, the chosen coupling box is 10 m wide and 14 m deep, as illustrated by the dark red lines on Figs 2 and 3. SPECFEM3D simulation domain is a cube 512 m in size with elements of 2 m, which amounts to ∼16.8 × 10⁶ elements. The propagation material has the same density and elastic wave velocities as the HOSS materials in their elastic domain: ρ = 1589 kg m⁻³ (see Table 1 in the SI), v_p = 1090 m s⁻¹ and v_s = 583 m s⁻¹. It has no attenuation nor gravity. Receivers are placed every 50 m vertically below the source and on concentric circles every 50 m on the free surface.

For coupling simulations, the centre of the coupling box corresponds to the centre of the top surface of the cube (x = 0, y = 0, z = 0). The coupling box occupies 7 elements in depth and 10 elements in width, for a total of 8081 GLL points.

For point-source simulations, the source is placed in the centre of the top surface of the cube (x = 0, y = 0, z = 0). We compute an approximation of the Green’s function of the material by simulating impulsive point-sources: impulsive forces are simulated using a Dirac delta source function and moment tensors using a Heaviside source function. The results of SPECFEM3D can then be safely convolved with the source time function derived from HOSS to obtain the correct point-source signal (Komatitsch & Tromp 2002). We generate one simulation for each separate component of the source. When using a point force, the Dirac impulse is represented by a triangle function with width 2 × dt and height 1/dt. This ensures that the total momentum of the source is 1 Ns, and facilitates scaling with HOSS’ source time function. Two simulations are run, with a force in the |$\vec{x}$| and |$\vec{z}$| directions respectively. When using a moment tensor source, the source is chosen to be SPECFEM3D internal Heaviside function, with a final value of 1 Nm. Four separate simulations are run for M_xx, M_yy, M_zz and M_xz. The total point-source signal is then the sum of displacements obtained for each simulated Green’s function convolved with corresponding source time functions F_i(τ) for the force and |$\frac{\partial {M}_{ij}}{\partial \tau }({\boldsymbol{\xi ^{*}}},\tau )$| for the moments.

In order to capture the propagation of P and S wave up to 512 m distance, the simulations are ran for 1.8 s. This is longer than HOSS simulations which last ∼120 ms, therefore both the point source parameters and the coupling waveforms calculated by HOSS need to be extrapolated in time before being used by SPECFEM3D. We extrapolate the source parameters time-series by fitting them with a sum of integrable pulse functions, with name ‘Jeffreys pulse’. This pulse function was shown to successfully fit impulses recorded in laboratory impact experiments (see e.g. Daubar et al. 2018). Details and results of the fit of force |${\boldsymbol{F}}({\tau })$| and moment rate |$\frac{\partial {\boldsymbol{M}}}{\partial \tau }({\boldsymbol{\xi ^{*}}},\tau )$| can be found in Fig. S8.

4.2 Results of the test of the point-source

4.2.1 Comparison of coupled and point-source signals

For both impactor scenarios, a comparison of point-source and coupled models is presented on Fig. 11. The receiver is placed at coordinates (283,283,0) m on the surface, at a distance of 400 m from the impact point. The non-rotated source series of M(ξ, τ) of Fig. 6(a) are used to compute the point source signals. We show both the displacement time-series and associated spectra of the P and the S waves, which are calculated from two separate tapered windows of the displacement time-series using a fast Fourier transform. We point out that for this surface sensor, the window used to evaluate the S wave also contains the Rayleigh wave contribution.

At 400 m distance, we observe a good match at low frequency between the point-source and coupled signals. The horizontal displacements are better matched than the vertical displacement series. The arrival times of the P and S waves are also matched by the point-source signal in the time domain. However, the overall ratio between the point-source amplitude and the coupled waveform amplitude is ∼5 and ∼15 for the P and S waves, respectively. This discrepancy is confirmed in the frequency domain: while the point-source spectrum of the P wave displays a trend similar to the coupled spectrum, it appears shifted by a constant positive factor. In addition, the point-source spectrum of the S wave presents a significant excess of energy, of more than an order of magnitude, with respect to the coupled waveform above 2 Hz. A similar tendency is observed at smaller distances from the impact. Signals recorded at 100, 200 and 300 m can be found in the SI (Figs S9–S11). The high frequency content of the P-portion and particularly the S-portion of the signal is already found in the Green’s functions of our medium. Indeed, the response of a half-space to impulse monopoles or step-like dipoles (the well-known Lamb’s problem) produces discontinuities at the onset of the Rayleigh wave (Johnson 1974; Kausel 2013; Feng & Zhang 2018). As shown in SI-5.2, these discontinuities are associated to an increase in high-frequency amplitude for the S-portion of the signal and may induce a higher corner frequency in the signal compared to the P-portion.

Fig. 12 shows the respective amplitude of different components of the point-source source in the total P-wave signal: the vertical and horizontal force terms, the isotropic (or explosive) term M^V and the deviatoric terms M^D. For both the vertical and oblique impact, the signal of the P wave appears to be dominated by the isotropic and deviatoric moment components (plain dark blue and green curves). The amplitude produced by the deviatoric term is comparable to the isotropic term, confirming the results of Table 1. It is interesting to note that the single force components, F_x and F_z, produce a lower amplitude signal than the moment source. In the low-frequency limit, their amplitude is more than one order of magnitude lower than the summed point-source signal, for both the vertical and oblique impact simulations. Whether this partition of energy between force and moments is maintained for impacts of higher velocity remains to be investigated. Indeed, we have shown that the moment is determined by the amount of inelastic damage, while the force is related to the amount of momentum transfer. Both processes might depend differently on the impactor velocity.

In the time domain (Figs 11a and b), the vertical coupled and point-source displacement signals (u_z) present a discrepancy at long time scales. Indeed, the coupled displacement signal appears to be converging to a negative, static value of vertical displacement of about −3e − 7 m at 400 m distance. This behaviour is not observed in any of the point-source u_z components produced with HOSS source time functions. A static displacement does also exist in the u_x and u_y coupled displacement signal and is correctly reproduced by the point-source, but its amplitude is about four times smaller than in the vertical direction. The reason for this discrepancy is not clear. it might indicate a missing analytical term in the definition of the force or moment tensor, which, in their current form, do not account for the permanent deformation of the free surface due to the crater. Using a simple convolution with a step function, we find that a constant vertical force of about −5e5 N, that is two orders of magnitudes smaller than the main force pulse, would be needed to fit the residual u_z displacement.

4.2.2 Simulating an extended source

Spectral and time-domain differences observed between the coupled and point-source signals in Section 4.2.1 could be due to the limits of the point-source description in seismology. Indeed, although the source model presented in Section 4.2.1 seems to produce satisfying results at low frequency for the u_x and u_y displacements, the convolution with point-source terms appears particularly inefficient at damping the high-frequencies Rayleigh wave in the elastic Green’s function. As presented in Section 2.3, one possible explanation is that a point-source is limited in its ability to model the interferences, or antenna effect, occurring within the finite source volume itself.

This finite-source effect determines the cut-off frequency and pulse shape of earthquakes (see e.g. Aki & Richards 2002; Madariaga 2015) and explosions (Denny & Johnson 1991). Their displacement spectra usually presents an omega-square (ω⁻²) roll-off above the cut-off frequency f_c. However, slow earthquakes (Supino et al. 2020), as well as explosions in weak materials (Ford et al. 2011), are sometimes better explained using an omega-cubed model (with ω⁻³ roll-off). Here, we measure the spectral characteristics of the coupled signal using a fit to the ω-squared (Ω₂) and ω-cubed (Ω₃) model of Aki (1967) and Brune (1970). Details on the fitting method can be found in Section SI-5.3.

We find that the Ω₃ model accomplishes overall a better fit to the coupled spectra than the Ω₂ model. However, the signal being contaminated by numerical noise above the cut-off frequency, it is difficult to discriminate with confidence between the two models. Using the Ω₃ model the cut-off frequencies of the P and S wave are f_{c, P} = 19.5 ± 6.2 and f_{c, S} = 24.6 ± 8.2 Hz. Considering that the cut-off frequency is limited by interferences between the two furthest points of the source, we can estimate the source size necessary to generate f_{c, P} and f_{c, S} by D_s = c/f_c. Taking the minimal velocity and maximal cut-off frequency of our system, we find that the seismic source must be at least v_s/f_{c, S} ≈ 23 m in size. This dimension is comparable, but larger than the stress glut region in Figs 2 and 3.

In order to further investigate the effects of the point-source hypothesis, we propose to simulate an extended source using SPECFEM3D and the results of HOSS vertical impact simulation. This approach revisits the computation method of Section 3.3, but using multiple smaller source volumes distributed in space. Instead of integrating the momentum and stress glut fields over the entirety of the HOSS volume, we start by cutting the simulation space into 19652 cubes, 1 m in size, and bin HOSS elements depending on which cube they belong to. Next, the source time function terms of eqs (18) and (14) are computed within each cube, providing 19652 new point-sources positioned at the centre of the cubes.Each source time function is then extrapolated to 1.8 s using a Hanning apodization method, and stored in SPECFEM3D data files. For the nine components of the source (i.e. F_x, F_y, F_z, M_xx, M_yy, M_zz, M_xy, M_yz and M_zx), we run a separate simulation containing 19652 sources and source time functions. The resulting signal is shown on Fig. 13 in purple, compared to the coupled signal in black and the point-source signal in red.

The spectra show clearly that the extended source simulated with method (1) is unable to sufficiently damp the high-frequency energy of the P and S wave. In fact, the point-source and extended source signals appears almost identical up to the frequency resolution of 140 Hz for the P wave. The extended source signal shows a reduced damping effect above ∼60 Hz for the S wave. This means that with our current source model, interferences occur within a radius smaller than v_p/140 ≈ 7.7 m for the P wave and v_s/60 ≈ 10 m for the S wave. These size estimates are more consistent with the dimensions of the stress glut fields displayed in Fig. 2, with the volumetric and deviatoric stress glut fields covering ∼5 and  10 m in diameter, respectively.

It is clear that a larger source would be required in order to reproduce the frequency content of the coupled signal. This observation leads to the conclusion that the source extent is likely underestimated by the current stress glut model. The explanation could be that additional inelastic processes, different from plasticity but present in HOSS numerical framework, have not yet been accounted for in our numerical definition of the stress glut. We discuss this point further in Section 5.

4.3 Towards Lunar and Martian data

We further test our point source and coupled waveform models by analysing their signal at regional distance, and by comparing the associated moment tensor and source duration to existing Lunar and Martian data.

4.3.1 Signal at regional distances

Here, we investigate the evolution of SPECFEM3D waveforms at regional distance (>100 km). Variations in seismic velocities and the effects of seismic attenuation will start to affect the signal after a few kilometres of propagation. These propagation effects can be modelled on Mars in a simple way with the current knowledge on the Martian subsurface, and using parameters estimated from the analysis of the recent impacts (Garcia et al. 2022). We propose to place the receiver at r = 100 km distance from the 1000 m s⁻¹ vertical impact. Signal recorded at such distance is known to be mostly propagating within the Martian crust, which can be considered homogeneous. Thus, the wave experiences a geometric attenuation in 1/r. We use 4000 m s⁻¹ for the velocity of P waves and 2310 m s⁻¹ for S waves, as in Garcia et al. (2022). Similarly, a quality factor Q_κ = 3500 is used for the attenuation of P waves, following the interpretation of the P-wave spectrum of events S0981c and S0986c, and |$Q_\mu =\frac{4}{9}Q_\kappa$|⁠. With these parameters, the P- and S-wave spectra at r₂ = 100 km distance can be estimated from the signal at r₁ = 400 m in the following way:

The resulting spectrums of the displacements |$\tilde{u}_z$| and the velocities |$\tilde{v}_z$| associated to the P and S waves are shown on Fig. 14. The simulated P-wave low-frequency amplitudes are about one to two orders of magnitude smaller than for event S0793a ( 3e − 9 m·Hz^−1/2 at 2 Hz) and S0986c ( 2e − 9 m·Hz^−1/2 at 2 Hz), with only  3e − 11 m·Hz^−1/2 for the coupled signal at 2 Hz. This is consistent with a smaller impact crater of about 2 m diameter for this simulation, compared to 3.9 and 5.7 m observed on Mars for S0793a and S0986c. The P-wave cut-off frequency measured in the attenuated spectra using the Ω₃ model is f_c = 11 ± 6 Hz (See SI-5.3 and Fig. S14). Although smaller than the value obtained in Section 4.2.2 without attenuation (f_{c, P} = 18.8 Hz), it is also higher than observed for InSight’s detected events, that is 9.4 Hz for S0793a and 8.0 Hz for S0986c. We also note that the S-wave spectra represented in Fig. 14 are strongly contaminated by attenuation above 10 Hz, thus the apparent cut-off frequency (also about 10 Hz) at this distance is less than the true source cut-off, f_{c, S} = 23.8 Hz and the high-frequency roll-off appears stronger than for the P wave.

4.3.2 Comparison model parameters with recorded Lunar and Martian Impacts

We wish to verify the consistency of the seismic moment computed by HOSS and estimates obtained for large Lunar impacts (Gudkova et al. 2011; Daubar et al. 2018) and for Martian events S0793a, S0981c, S0986c and S1094b (Garcia et al. 2022; Posiolova et al. 2022). As a base for our comparison, we use the scaling relationship between seismic moment and impactor momentum, which has been demonstrated in several previous studies (Gudkova et al. 2015; Wójcicka et al. 2020).

The amplitude of seismic motion generated by an event depends on the value of the source seismic moment, M₀, but also on the mechanical properties of the source region (see e.g. Aki & Richards 2002). Therefore, in order to properly compare seismic moments from impact sources on different bodies with varying surface materials, as on the Moon, on Mars, and in a numerical simulation, one must first establish a common reference source material. Daubar et al. (2018) supposed that the impact-generated seismic wave was measured in a strong bedrock layer beneath the impact region, with density ρ_br and P-wave velocity v_{p, br}. A similar approach was followed by Posiolova et al. (2022) when comparing the seismic moment estimated from P-wave amplitudes using a seismic model at 50 km depth on Mars, with the seismic moments computed in a regolith or fractured basaltic material analogous to the Martian surface.

Following (Daubar et al. 2018), an impact with seismic moment M₀ happening in a layer with density ρ and P-wave velocity v_p corresponds to an equivalent moment M_{0, br} in underlaying bedrock given by :

The first fraction corresponds to an amplitude scaling term, and the second fraction to the transmission factor for seismic waves going from the unconsolidated source material to the solid bedrock.

We compare Lunar, Martian and simulated seismic moments by calculating their equivalent bedrock moment M_{0, br} using eq. (35). The obtained moments M_{0, br} are divided by v_br to obtain results in Ns, and plotted against the impactor momentum P_i in Fig. 15. The methods, parameters and references used to produce this plot are further described in Appendix  E. To the moments calculated previously for the 1000 m s⁻¹ oblique and vertical impacts, we add the results of impact simulations of Perseverance’s entry balance masses, performed by Fernando et al. (2021). The moment of this 4000 m s⁻¹ impact was calculated for a vertical impact incidence, and with the same material models and methods as this study. We exclude the results of the oblique scenario in Fernando et al. (2021) due to the differences in mesh resolution and size with the present work. Despite the significant discrepancies in the calculation methods for M₀ between Martian, Lunar and numerical studies, the graph indicates that our results are in trend with the seismic moments calculated for small Martian Impacts and with the scaling determined for Lunar impacts, for two distinct velocities.

Previous studies have also proposed a comparison of impact source duration estimates, defined as τ = 1/f_c, on the Moon and Mars (Gudkova et al. 2015; Garcia et al. 2022; Posiolova et al. 2022). We have obtained this parameter in Section 4.2.2 using the coupled waveform spectra of our impact simulations. However, again, Gudkova et al. (2015) have shown that this seismic parameter is sensitive to surface material properties, such as regolith porosity. Garcia et al. (2022) also observed that Lunar impact events are characterized by a longer source duration than similar Martian events. The authors relate the longer Lunar impact source duration to the lower seismic velocities of the Lunar regolith, which result in slower source dynamics. To be able to better compare sources properties on Mars, the Moon and in our simulations, we propose to base our comparison on another useful scaling variable for seismic sources: the estimate of the source size, |$\left(\tau \, v_s\right)$|⁠, that is the product of the source duration with the S-wave velocity at the source location. Just like the scaling of Fig. 15, this multiplication by v_s aims at reducing the bias associated to the difference in surface properties between each impact type.

Some hypotheses must still be made when choosing the values of v_s and τ for the Moon and Mars. In the same way, the estimation of the source duration is not straightforward and is affected by choices of methods for spectral estimates, and choices of models for the scattering and attenuation present in real seismic data. The different approaches of Martian and Lunar studies are presented in Appendix  E. With these limitations in mind, we plot the source size |$\tau \, v_s$| as a function of the impactor momentum P_i and impactor kinetic energy E_i on Fig. 16. On both scalings, we find that Martian and Lunar events follow a similar trend, with an overall increase of the source size with impactor energy or momentum. Source size estimates align with impactor momentum and energy with a slope of 0.13 and 0.11, respectively. The numerical simulation estimate is again consistent with Lunar and Martian data.

5 DISCUSSION

The model developed in Section 2 provides an exhaustive description of the impact seismic source. Contrary to previous studies, which proposed models based on the seismic moment of explosions, or on the seismic impulse, we introduce an expression of impact-generated displacements which combines a 9-component moment tensor and a vector force, both varying in time. In support of this model, we develop a new numerical method to compute the stress glut, an essential parameter in the seismic source. This mixed analytical/numerical approach is able to better represent the temporal and mechanical complexity of the impact phenomenon. It also offers different levels of approximations, from the point source to the extended source representation.

A key finding of this study is that the impact seismic source cannot be rigorously described with only an impulse, nor with only its seismic moment: in fact, its point source expression is a combination of both. For both vertical and oblique impacts, the deviatoric part of the moment tensor has amplitudes comparable to the volumetric part. Hence, a pure isotropic explosion in not sufficient to describe impacts. In particular, the source of an oblique impact is complex and lies between isotropic, CLVD or crack and DC mechanisms. As shown in Table 1, the force and volumetric source components of an oblique impact are well approximated by multiplying the results of a vertical impact by sin (45°). This interesting observation suggests that part of the source could be proportional to the vertical projection of the impactor momentum, as suggested for example by Wójcicka et al. (2020). The moment and force components have different time scales related to different stages of the cratering process: the very fast transfer of momentum dominates the force, while the generation of damage and escape of ejecta dominate the moments. We note that in the impact scenarios presented here, the contribution of the vector force to far-field motion is limited, and the signal appears dominated by the dipole of the source.

The method proposed in this study relies on non-linear impact simulations to compute the seismic source terms. As such, the definition of the stress glut in the numerical model is key to properly retrieve the source. At low frequency, the seismic signal obtained using this stress glut model agrees with the prolongated non-linear signal to within an order of magnitude. However, at high frequency, we note strong discrepancies between the spectra obtained from the stress glut and from the coupling. Point source signals (and generally, speaking, extended-source signals) sum the convolutions of different source time functions with the Green’s function of an ideal-elastic half-space medium. The Green’s function of a half-space is the solution to Lamb’s problem (Johnson 1974) with impulse or step sources. If the source and receiver are both close to the surface, resulting displacements contain discontinuities and are dominated by a strong Rayleigh wave (see SI-5.2). The source time functions obtained from the stress glut are not presently able to damp this high-frequency energy and suppress the Rayleigh wave. We pinpoint next the current limitations of the stress glut model and possible causes for the lack of damping.

A first limitation lies in assumptions and approximations made in the calculation of the stress glut. In Section 3, we note that element finite rotation in space affects stress glut measurements significantly. The problem of rotation effects ties in to a fundamental question in material mechanics, which is the definition of stress and strain for finite deformation. The classical stress measure in seismology is the Cauchy Stress Tensor, which describes the forces applied to a volume element in its deformed configuration. Other stress measures, like the Second Piola–Kirchhoff stress tensor, can also be used to better describe pre-stressed media (see e.g. Dahlen & Tromp 1998). However, these stress measures are only valid in cases where deformation can be considered infinitesimal, with the infinitesimal strain tensor defined as |$\varepsilon _{ij} = \frac{1}{2}(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i})$| (Dahlen & Tromp 1998; Aki & Richards 2002). In impact simulations, deformations can have a scale comparable to the typical size of the volume element, and infinitesimal strain theory is no longer applicable. Instead, finite strain theory describes deformations using the deformation gradient tensor introduced in Section 3.3, which relates the current configuration of a volume element to its initial configuration. It is unclear at this time whether the difference in stress and strain definition in classical seismology and HOSS numerical simulation could introduce some errors in the definition of equivalent forces, and how the finite rotation of the elastic material at the source might influence the seismic signal.

Secondly and more importantly, a limit of this study lies in the fundamentally different mathematical frameworks used in classical seismology, with respect to the FDEM used in HOSS. Indeed, while equivalent forces derived in eq. (8) are valid in a purely continuous world, the reality of material deformation and the approach followed by HOSS includes discontinuities. In Section 3.1.2, we explained that we purposely left out the effect of inter-element fracture and friction in our description of the source. This is a strong approximation, as we know that friction mimics a plastic process at play within granular geomaterials. Moreover, friction is a damping mechanism, or energy sink, that might contribute to the reduction of high-frequency energy within the source region. Its absence in the point-source model could explain the excess of high-frequency energy in the signals of Figs 11, 12 and 13. To measure its effect on the source, inter-element friction should be accounted for in the expression of equivalent forces. However, given that this process is considered to be a surface action in HOSS’ numerical framework, it cannot be simply represented in the form of an additional stress glut tensor. In fact, the action of friction forces would be better understood by rewriting the equations of Section 2.1.2 in the form of a mass-spring system rather than with the equations of continuum mechanics. Due to the extensive reformulation required to include this additional frictional and discrete processes in our description of the source, we leave it for a future study.

We also want to emphasize the discrepancy in vertical displacement between the coupled and point-source signals. Such difference hints at the fact that our analytical model might need to be refined in order to accommodate a permanent vertical force. Further testing of the point-source representation could involve including higher-order moments, such as quadrupoles, in the Green’s function expansion of eq. (10). Higher-order moments have indeed a strong potential to investigate finite-source effects and mechanism complexity (Stump & Johnson 1982; Jordan & Juarez 2019). We also point out the limitation of the Green’s function representation itself. Indeed, this method applies to ideal elastic media with well known boundary conditions. The definition used here assumes a plane and stable free surface at z = 0. But the impact cratering process involves a significant disruption of the true free surface. The approximation we make of a plane and constant boundary could explain why a strong impulsive Rayleigh wave is observed in synthetic Green’s function or point source signals, and not in coupled signals.

To insist on analytical developments, we recall that in Section 4.2.2, an ω⁻³ seismic source spectrum was used to model the high-frequency content of the displacement spectra. Although this model is intuitively appealing for a 3-D seismic source, it does not directly relate to impact dynamics or to their equations of motion. Another possible improvement of the source model, mentioned in Section 2.3, would be to analytically derive the effect of the source extension on the spectrum, in a similar way as Savage (1966) or Haskell (1966). A first order solution could use analytical expressions of the Green’s function for Lamb’s problem as developed for example by Johnson (1974), Kausel (2013) or Feng & Zhang (2018).

Another example of possible limitations is the omission of gravity in the analytical expressions of moments (Section 2) and in simulations (Section 3). Gravity is known to have an influence on the dynamics of crater and ejecta formation (Holsapple 1993; Froment et al. 2020): for example, it limits the duration of the crater growth and collapse (e.g., Melosh 2011). Thus, the inclusion of gravity might affect the long-term decay of the source time functions in Fig. 7. In Appendix  D, we propose a modification of equivalent forces |${\boldsymbol{\gamma ^V}}$| and |${\boldsymbol{\gamma ^S}}$| to include the non-linear effects of gravity in the stress glut theory. It is also possible to include gravity in HOSS simulations, although at an increased computational cost. Indeed, adding a constant vertical acceleration requires giving the simulated material enough time to relax to its lithostatic equilibrium, and we leave this study for future investigations. However, even without gravity, simulation initial and boundary conditions influence long-term dynamics, as evidenced by the slow drift in our modelled source time function.

Despite the frequency content differences between modelled and prolongated signals, we show that some key parameters of the modelled impact seismic source match global trends observed on Mars and the Moon. We emphasize that comparisons between Lunar, Martian and simulation data is challenging, due to the absence of direct measurements of impactor or material properties. For example, in the data presented above, impactor momentum was inferred from crater measurements (Collins et al. 2022) or from signal amplitudes (Gudkova et al. 2011) using different models. Seismic velocities are estimated from a variety of Martian and Lunar seismic models, which at this time cannot account for local variations on different parts on the planetary surface. Still, we note that the scaled seismic moments and source size estimates obtained with HOSS are in trend with Lunar and Martian results.

This is to our knowledge the first time that a trend is evidenced between the impactor energy and momentum and the estimated source size. Fig. 16 indicates that the source size estimate for Lunar and Martian impacts scales as the impactor momentum and kinetic energy to the power 0.13 and 0.11, respectively. This result can be connected to existing scaling laws for impact cratering. It is known that impact crater size or depth does not scale directly with energy |$(E_i=r_i^3 \rho _i v_i^2/2)$| or momentum |$(p_i=r_i^3 \rho _i v_i)$|⁠: instead, the pi-scaling introduced by Holsapple (1993) proposes that crater dimensions scale with a mixed point-source measure |$\mathcal {C}=r_i v_i^\mu \rho _i^\nu$|⁠, with 1/3 < μ < 2/3 and ν ≈ 1/3. For an impact in the strength regime in Lunar regolith or a dry soil, crater size scales with energy to the power 3μ/6 = 0.2. This is a power greater than the one observed in Fig. 16, which suggests that our source size estimate is not directly proportional to crater size. This observation, although preliminary, shows that further investigations of the impact seismic source on planetary bodies is needed to be able to relate it to classical scaling laws for impacts, or for seismic sources on Earth.

In particular, this study focused on only two impact scenarios with a common impactor velocity of 1000 m s⁻¹. This impact velocity is in the lower than the mean Martian impact velocity by a factor of ∼10, and than the mean lunar impact velocity by a factor of ∼20 (Le Feuvre & Wieczorek 2011). To be applicable to real impact scenarios, our model will need to be tested for a range of target materials, impactor velocities and angles closer to the observed range, as was initiated with the modelling of Perseverance entry sequence (Fernando et al. 2021).

6 CONCLUSION

We introduce an analytical model relating the mechanical fields of the meteorite impact phenomena (i.e. velocity field, plastic and elastic stress fields in the shocked material) to the seismic displacements recorded at any distance from the formed crater using the seismic Representation Theorem and the stress glut theory. A point-source model of the impact is obtained, which associates a time-varying vector force and a time-varying moment tensor exerted at the impact point. We subsequently develop a numerical method to estimate the different terms involved in the seismic source. One of these terms, the stress glut field, is calculated from the plasticity measured in 3D non-linear impact simulations performed with the HOSS software.

We test this numerical model by comparing signals produced by the point source in an elastic medium, with signals obtained by solving the non-linear equations of motion with HOSS and coupling the solution to SPECFEM3D. The comparison reveals that the modelled P- and S-wave signals agree with the coupled signal to within an order of magnitude at low frequency. Point source signals present significantly higher amplitudes at high frequency. This discrepancy is due, in part, to the impulsive Rayleigh wave present in the ideal elastic response of a half-space.

This coupling method allows us to study the respective contribution of various source terms on the impact seismic signal. We show that the source is mostly dominated by the moment tensor components, and that the deviatoric part of the moment tensor contributes significantly to the impact seismic signal. Hence, an isotropic (explosive) source does not provide a complete description of the impact source mechanism. Numerical results make it possible to study the effect of source spatial extension on the signal spectra. We show that, for a vertical, 1000 m s⁻¹ impact simulation generating a 2 m crater, finite source effects are not sufficient to explain the lack of high-frequency energy in the coupled spectrum compared to the modelled point-source spectrum. We hypothesize the absence of some stress glut terms in our numerical description, in particular terms associated to the specificities of the FDEM implemented in HOSS. Deviations from the elasto-dynamic equations due to numerical damping, interelement friction and dislocation, or the finite strain theory may be viewed as additional equivalent forces whose importance should be assessed in the future.

Despite these discrepancies, we have shown with a simple propagation model that the properties of the impact generated P wave 100 km away from the source is in line with spectral properties of small impact recorded seismically during the InSight mission. We also showed that, once scaled with the material seismic properties, the measured seismic moment and source duration agree well with measurements made on Lunar and Martian data. The comparison also reveals a possible scaling relationship between seismic source size and impactor energy and momentum. In the future, we hope to conduct a more complete validation study of the seismic source parameters, and to further investigate the scaling of these key source parameters.

The proposed model is here applied to impact phenomena, which are a rather exotic source from the point of view of Earth seismology, but the developments remain true for any other type of source, such as explosions, volcanoes or of course the earthquake sources, for which the stress glut was initially invented by Backus & Mulcahy (1976a). In fact, the stress glut field and equivalent forces can be computed for any mechanical disruption in a solid medium, as long as the right initial and boundary conditions are provided.

Acknowledgement

We would like to thank editor Carl Tape, reviewer Brian W Stump and an anonymous reviewer for their insightful remarks and questions, which helped us improve the quality of this paper. This work is InSight Contribution Number 219 and LA-UR-22-30392. This research is supported by ANR MAGIS (ANR-19-CE31-0008-08) and the Initiative d’Excellence (IdEx) Université Paris Cité (ANR-18-IDEX-0001). It used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001. ZL and ER acknowledge the support from Los Alamos National Laboratory LDRD program (grant number 20230310ER). MF, CL and ER received funding from the Center for Space and Earth Science (CSES) of LANL. CSES is funded by LANL’s Laboratory Directed Research and Development (LDRD) program under project number 20210528CR. MF received support from the Research Council of Norway (FRIPRO project 335903, AIR) for part of this research. Numerical simulations presented in this work were performed with the Hybrid Optimization Software Suite (HOSS) and the SPECFEM3D software. HOSS is available under commercial licensing, and an educational version, HOSSedu, is accessible for universities only. More information can be found on https://www.lanl.gov/hoss. SPECFEM3D is part of a family of open-source spectral element method codes designed for computational seismology. Information on, access to and support for the code can be found on https://specfem.org/ and https://github.com/SPECFEM/specfem3d#readme.

DATA AVAILABILITY

Parameters, geometries and initial and boundary conditions used in the simulations of this paper are provided in a Github repository: https://github.com/m-froment/Impacts_Stress_Glut.git.

References

APPENDIX A: EQUATIONS OF MOTION WITH SURFACE MASS LOSSES

The equation of motion of a variable-mass system is given by the generalized form of the Reynolds Transport Theorem. This theorem was introduced by Osborne Reynolds (Reynolds 1903), and is commonly known in a restricted form applicable to material volumes, that is volumes surrounding a constant set of particles in motion, whose surface moves together with the outermost particles. The generalized form accounts for cases where mass (and thus particles) is allowed to flow through the surface: the surface thus moves at a velocity |${\boldsymbol{v}}^\Sigma$|⁠, distinct from the particle velocity |${\boldsymbol{v}}$|⁠. Using the formalism of Irschik & Holl (2004) (eq. 4.6), the generalized Reynolds Transport Theorem writes:

In eq. (A1), V^t and Σ^t refer to the variable volume and surface at time t and V designates the instantaneous material volume composed of the particles of V^t. Following Irschik & Holl (2004), the derivative |$\dfrac{{\text{d}}_\Sigma }{{\text{d}}t}$| means that the measure of the total momentum within V^t takes into account the inflow and outflow of mass through Σ^t, while the derivative |$\dfrac{{\text{d}}}{{\text{d}}t}$| simply considers variation of momentum of the set of particles instantaneously present in V. Using Cauchy’s momentum equation on this material volume, the following expression of the conservation of momentum is obtained:

Here, |${\boldsymbol{S}}$| represents the true, non-linear stress exerted in the strongly shocked medium. It is different from the ideal elastic stress |${\boldsymbol{\Psi }}$| of eq. (3) and does not follow Hooke’s law of elasticity. Similarly, |${\boldsymbol{h}}^V$| stands for the non-linear volume forces applied to V^t. This global expression of the conservation of momentum can be completed by a local form on an mesoscopic volume element. To this aim, we make use of eq. (2.3d) of Irschik & Holl (2004):

and write:

Upon using the Gauss–Green–Ostrogradsky divergence theorem on the surface integrals of eq. (A4) and projecting in direction i, this yields:

Finally, true non-linear tractions on the surface Σ^t now write:

where |$f_i^\Sigma$| denote external forces applied to surface Σ in direction i and n_j is the normal to surface Σ projected in direction j.

APPENDIX B: REPRESENTATION OF A SEISMIC WAVEFIELD IN A VOLUME WITH CONSTANT MASS

Let V be a volume with surface Σ. Let |${\boldsymbol{u}}({\boldsymbol{\xi }},\tau )$| be an elastic displacement field generated by surface tractions |${\boldsymbol{\Psi }}({\boldsymbol{u}},\tau )\cdot {\boldsymbol{n}}={\boldsymbol{f}}^\Sigma$| on surface Σ with normal |${\boldsymbol{n}}$|⁠, and volume forces |${\boldsymbol{f}}^V$| within V. Let |${\boldsymbol{v}}({\boldsymbol{\xi }},\tau )$| be a second displacement field, produced by different tractions |${\boldsymbol{\Psi }}({\boldsymbol{v}},\tau )\cdot {\boldsymbol{n}}={\boldsymbol{g}}^\Sigma$| on Σ and volume forces |${\boldsymbol{g}}^V$| in V. The equations of motion for |${\boldsymbol{u}}$| and |${\boldsymbol{v}}$| are:

Betti’s Reciprocal Relation, which is valid everywhere within V, rearranges the elastic equations of motion for |${\boldsymbol{u}}$| and |${\boldsymbol{v}}$| as:

Note that this relationship is also true when both V and Σ vary in time. Classically, the next step consists in integrating both part of the equation from time τ = −∞ to time τ = +∞. We can further add the condition that |${\boldsymbol{u}}$|⁠, |$\frac{\partial {\boldsymbol{u}}}{\partial \tau }$|⁠, |${\boldsymbol{v}}$| and |$\frac{\partial {\boldsymbol{v}}}{\partial \tau }$| are all everywhere zero before a certain time τ₀ in the past, and evaluate the field |${\boldsymbol{v}}$| at a time t − τ, where t is a fixed time. In the case where V and Σ are fixed volumes and surfaces, the time integrals over the acceleration terms |$\rho \frac{\partial ^2 {\boldsymbol{u}}}{\partial \tau ^2}\cdot {\boldsymbol{v}}$| and |$\rho \frac{\partial ^2 {\boldsymbol{v}}}{\partial \tau ^2}\cdot {\boldsymbol{u}}$| cancel each other, and the following expression of Betti’s theorem is obtained:

Note however, that if volume V and surface Σ are considered to be varying in time, the simplification of the acceleration terms must be carried more carefully. For this special case, additional analytical terms appear, which were derived by Minster (1974) and Archambeau & Scales (1989). In this study, we do not account for the effects of variable mass and volume, as we will find them to be negligible for the studied impact. However, the reader can find an exact version of Betti’s relation for variable volumes and surfaces in the Appendix  C below.

In a last step, |${\boldsymbol{v}}$| is chosen to be the Green’s function of the propagation medium, |$v_i({\boldsymbol{x}},t)=G_{in}({\boldsymbol{x}},t-\tau ,{\boldsymbol{\xi }},0)$|⁠. It represents the ith component of displacement produced at time t − τ and position |${\boldsymbol{x}}$| within V by an impulse volume force located at position |${\boldsymbol{\xi }}$| and time 0 and directed towards the nth direction of space. The volume force and boundary conditions associated to |${\boldsymbol{v}}$| are:

Reintroducing this new expression of |${\boldsymbol{v}}$| in eq. (B2) gives rise to the Representation Theorem, here written in the nth direction of motion and for a constant volume V:

A list of several special-case depending on various boundary conditions can be found in Aki & Richards (2002). In the case considered here, Σ is chosen to be a free surface, which leads to |$C_{ijkl}n_j\dfrac{\partial G_{kn}}{\partial \xi _l}=0$| everywhere on Σ and cancels the first part of the last term of eq. (B5). We replace the generic forces |${\boldsymbol{f}}^V$| and |${\boldsymbol{f}}^\Sigma$| by equivalent surface and volume forces |${\boldsymbol{\gamma }}^V$| and |${\boldsymbol{\gamma }}^\Sigma$| which are the non-linear sources of motion (see eq.7 of Section 2.1.2). The expression of displacements anywhere in V is then given by eq. (B6):

APPENDIX C: REPRESENTATION THEOREM FOR A VARIABLE-MASS SYSTEM

As mentioned in the previous section (appendix  B), Betti’s Representation Theorem as found in Aki & Richards (2002) does not account for changes in volume and surfaces. We report here on the form obtained by Minster (1974) and Archambeau & Scales (1989). The starting point is Betti’s Reciprocal Relation (eq. B2), which we copy here:

In the following step, both parts of the equation are integrated between time τ = −∞ and time τ = +∞. We add the condition that |${\boldsymbol{u}}$|⁠, |$\frac{\partial {\boldsymbol{u}}}{\partial \tau }$|⁠, |${\boldsymbol{v}}$| and |$\frac{\partial {\boldsymbol{v}}}{\partial \tau }$| are all everywhere zero before a certain time τ₀ in the past, and evaluate field |${\boldsymbol{v}}$| at a time t − τ, where t is a fixed time. This time, instead of adopting the simplification of eq. (B3), we must further develop the time integral of momentum over a time-varying volume, V^τ:

We modify the double derivatives on the right and write:

Only two terms are left- on the right-hand side of eq. (C3). As was proposed by several authors (Minster 1974; Archambeau & Scales 1989), eq. (C2) can be developed using the Reynolds Transport Theorem of Appendix  A. Following eq. (2.3d) of Irschik & Holl (2004), this gives:

In eq. (C4), |${\boldsymbol{v}}^\Sigma$| represents the velocity of the moving surface Σ^τ. Similarly to the case of fixed volume and surfaces (see Appendix  B), the first integral on the right-hand side cancels upon applying the right initial conditions on |${\boldsymbol{u}}$| and |${\boldsymbol{v}}$| and their derivatives, and evaluating |${\boldsymbol{v}}$| at time t − τ. We are left with a supplementary term to Betti’s theorem, which now writes:

Introducing the Green’s function in eq. (C5), as in eq. (B5) of Appendix  B, we obtain the exact version of the Representation Theorem, here written in the nth direction of motion and for a varying volume V^τ:

We note that an additional artificial surface force |$\tilde{f}^\Sigma _i({\boldsymbol{\xi }},\tau )=\rho v_iv_j^\Sigma ({\boldsymbol{\xi }},\tau )n_j$| and an artificial surface stress |$\tilde{\Psi }_{ij}({\boldsymbol{\xi }},\tau ) = \rho \frac{\partial G_{in}}{\partial \tau }({\boldsymbol{\xi }},t-\tau ,{\boldsymbol{x}},0)v_j^\Sigma ({\boldsymbol{\xi }},\tau )$| appear in the Representation Theorem.

APPENDIX D: STRESS GLUT ON A SELF-GRAVITATING, ROTATING PLANET

Modelling wave propagation in a planet subject to its own gravity and rotation requires us to adapt the equations of motion of Section 2.1.1. Under the effect of gravity and rotation, the planet is initially in a state of hydrostatic equilibrium, which implies the existence of a pre-stress. Elastic deformations represent perturbations of this pre-stress, which the traditional Cauchy stress tensor cannot appropriately describe. Moreover, material put in motion by a seismic wave can itself redistribute mass and perturb the gravitational field, which in return acts on the wave dynamics. The equations of motions of a self-gravitating and rotating body were presented in a number of works. Lognonné & Clévédé (2002) presented a review of the theory of normal modes, starting from the gravito-elastic equations of motion. Dahlen & Tromp (1998) exhaustively address the issue of the definition of stress in continuum mechanics, as well as the derivation of appropriate linearized equations of motion and boundary conditions.

In most applications, the studied planet is considered to be hydrostatic, or quasi-hydrostatic. This implies a distribution of density which is roughly laterally homogeneous. Within such approximation, a relationship exists between the equilibrium stress, gravitational and centrifugal potentials, and the equations of motion can be made independent from the previous history of stresses.

Let us consider a planet of volume V and surface Σ. In its initial hydrostatic equilibrium, the planet has a density field ρ⁰, an initial static Cauchy stress tensor |${\boldsymbol{T}}^0$| and rotation vector |${\boldsymbol{\Omega }}$|⁠. In the following section, for the sake of concision, we adopt vector notations. At equilibrium, the equilibrium gravitational field is:

with ϕ⁰ the gravitational potential defined by Poisson’s equation using the gravitational constant G:

The centrifugal potential ψ of the rotating planet is defined by:

Finally, the equations of the hydrostatic equilibrium is:

The onset of a seismic wave perturbs the initial position of particles. Using a Lagrangian description of the motion, this perturbation can be written:

with |${\boldsymbol{x}}$| the initial position of particles and |${\boldsymbol{u}}$| the Lagrangian displacement vector. In reaction to the motion, we consider that other physical fields experience first order perturbations, such that |$q^L({\boldsymbol{x}},t) = q^0({\boldsymbol{x}},t) + q^{L1}({\boldsymbol{x}},t)$| in the Lagrangian description and |$q^E({\boldsymbol{r}},t) = q^0({\boldsymbol{r}},t) + q^{E1}({\boldsymbol{r}},t)$| in the Eulerian description. q^L1 is the first order Lagrangian perturbation of quantity q^L, and is related to the first order Eulerian perturbation by |$q^{L1} = q^{E1} + {\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}q^0$|⁠, which is a form of linearized and integrated material derivative. With these notations, the conservation of mass in V is written:

The first order Eulerian perturbation of the gravitational potential is:

associated to a first order perturbation of the gravitational field |${\boldsymbol{g^{E1}}} = -{\boldsymbol{\nabla }}\phi ^{E1}$|⁠. The full, non-linearized equation of motion is more easily written in the Eulerian form:

Upon linearizing each field as above, and neglecting second and higher-order terms, it becomes:

The Lagrangian perturbation in stress, |${\boldsymbol{T}}^{L1}$|⁠, is more useful than |${\boldsymbol{T}}^{E1}$| in continuum mechanics as it is directly related to the gradient of deformation |${\boldsymbol{\nabla }}{\boldsymbol{u}}$|⁠. Precisely, |${\boldsymbol{T}}^{L1}$| is the incremental Lagrangian Cauchy stress, which in linearized form writes |${\boldsymbol{T^{L1}}} = {\boldsymbol{C}}:\frac{1}{2}\left[{\boldsymbol{\nabla }}{\boldsymbol{u}} + ({\boldsymbol{\nabla }}{\boldsymbol{u}})^T \right]$| with |${\boldsymbol{C}}$| is the stiffness tensor. |$T^{L1}_{ij}$| is thus equivalent to the ideal stress Ψ_ij defined in the main text. Upon applying the relationship between Lagrangian and Eulerian perturbations and using the hydrostatic equilibrium, we can finally rewrite eq. (D9) as:

In tensor notation, eq. (D10) becomes:

The previous developments bring a few additional terms to the version of the equation of motion proposed in eq. (3) of Section 2.1.1. In fact, the simple volume force |$f_i^V$| of eq. (3) is now expressed as a function of displacements |${\boldsymbol{u}}$|⁠: |$f_i^V = \rho ^0g_i^{E1} -2\rho ^0\epsilon _{ijk}\Omega _j\frac{\partial u_k}{\partial t} - \partial _i \left( \rho ^0u_j\partial _j(\phi +\psi )\right) - \rho ^{E1}\partial _i(\phi + \psi )$|⁠. This dependance of |$f_i^V$| on |${\boldsymbol{u_i}}$| means that the definition of the Green’s function changes with respect to Appendix  B, as G_in in no longer solution to the same equation of motion. The newly defined Green’s function depends on a symmetric operator of gravito-elasticity and an antisymmetric operator for Coriolis forces (see e.g. Lognonné & Clévédé 2002). A demonstration of the Representation Theorem using this form of the linearized equations of motion can be found in Dahlen & Tromp (1998, section 5.3).

Changes in the equations of motion transpose to the definition of equivalent forces. As in Section 2.1.2, using ρ^E = ρ, |${\boldsymbol{\gamma }}^V$| becomes:

We identify in this expression the stress glut |${\boldsymbol{\Pi }} = {\boldsymbol{T^{L1}}}-{\boldsymbol{S}}$|⁠.

APPENDIX E: PARAMETERS AND METHODS FOR THE SCALING OF SEISMIC MOMENTS AND SOURCE SIZE ESTIMATES

In Section 4.3.2, we proposed a method to scale seismic moments M₀ obtained by several studies to a common reference M_{0, br} measured in bedrock. We also proposed a definition of the impact source size estimate based on the product of the source duration and the source shear wave velocities. We report here the references and parameter values used to produce Figs 15 and 16.

Moment scaling: Simulation, Martian and Lunar moments are converted to an equivalent bedrock moment M_{0, br}, considering a bedrock with density ρ_br = 2700 kg m⁻³ and P-wave velocity v_{p, br} = 1000 m s⁻¹. The material properties ρ and v_p of the source layer are chosen so as to match the seismic models used in each study when computing M₀. For example, the Lunar impacts are assumed to have occurred in a material with v_p = 330 m s⁻¹ and ρ = 2000 kg m⁻³, as in Daubar et al. (2018). Other values of v_p and ρ for Martian impacts and our simulations are reported in the third and fourth columns of Table E1, with associated references.

The momentum of Lunar and Martian events displayed in Fig. 15 were obtained in multiple ways by previous studies. The momentum of Lunar impact events is given by Gudkova et al. (2015). Their unscaled seismic moment M₀ is calculated from |$M_0 = v_p\, S\, (P_i)$|⁠, with S = 1.5 the ejecta amplification factor and v_p = 330 m s⁻¹. The momentum of Martian impact events S0793a, S0981c and S0986c were estimated by Garcia et al. (2022) using a statistical model of meteoroid entry (Collins et al. 2022). Their seismic moment was estimated by scaling the amplitude of simulated waveforms. These simulations were performed in a reference model with a surface sedimentary layer given by Table E1. The seismic moment of the large event S1094b was estimated by Posiolova et al. (2022) in two different materials using a scaling relationship between seismic moment and crater diameter from Wójcicka et al. (2020). In Fig. 15, the scaling relationship found by Wójcicka et al. (2020) for a material with v_p = 1088 m s⁻¹ is converted to the Lunar seismic velocities using a multiplication by a factor 330/1088.

Source size scaling: The estimate of source size also requires a choice of seismic velocities. For Lunar and Martian impacts, we adopt values of v_s found in the literature for a source depth of approximately 10 m, and 20 m for the largest Martian impact S1094b. Corresponding values are reported in Table E1). This is a first order estimate, which does not account for variation of v_s with depth, or the variation of source depth with the size of the crater.

When using real seismic data, the estimation of the source duration or cut-off frequency is not straightforward. Indeed, its value depends on the time window chosen to compute the source spectra: long time windows will contain a mixture of P and S waves with potentially different source time scales. The determination of τ is also affected by scattering and attenuation phenomena at regional distances. In Gudkova et al. (2015), source durations of large Lunar events were computed on long time windows using a fit to a function |$\hat{s}(\omega ,\tau )$| which included an attenuation model and an ω⁻³ roll-off at high frequency. To better match these results, we use the fit of the Ω₃ model to the coupled S wave as our estimate of the source duration τ_S = 1/f_{c, S}. For small Martian events (Garcia et al. 2022), the source duration was estimated from a fit of a Ω₃ model to the first ∼5 s of signal arrival, in order to isolate the P-wave information and limit contamination from seismic scattering. We suppose that these Martian results have a ±2 Hz uncertainty.

Finally, the kinetic energy of Fig. 16 is calculated from the momentum estimate, E_i = P_i × v_i/2, with typical values of v_i given in Table E1.

Author notes