An estimate of the impact rate on Mars from statistics of very-high-frequency marsquakes

The number density of impact craters on a planetary surface is used to determine its age, which requires a model for the production rate of craters of different sizes. On Mars, however, estimates of the production rate of small craters (<60 m) from orbital imagery and from extrapolation of lunar impact data do not match. Here we provide a new independent estimate of the impact rate by analysing the seismic events recorded by the seismometer onboard NASA’s InSight lander. Some previously confirmed seismically detected impacts are part of a larger class of marsquakes (very high frequency, VF). Although a non-impact origin cannot be definitively excluded for each VF event, we show that the VF class as a whole is plausibly caused by meteorite impacts. We use an empirical scaling relationship to convert between seismic moment and crater diameter. Applying area and time corrections to derive a global impact rate, we find that 280–360 craters >8 m diameter are formed globally per year, consistent with previously published chronology model rates and above the rates derived from freshly imaged craters. Our work shows that seismology is an effective tool for determining meteoroid impact rates and complements other methods such as orbital imaging.

Supplementary Material Figure 2: Martian impact rate comparison as a cumulative size-frequency distribution (SFD) plot.The left and right axes show the yearly rate per km 2 and globally, respectively.The bottom and top axes show the crater diameter and corresponding M W . Shown are values N(≥ 16 m) for Neukum et al. [49] (grey upwards triangle), Hartmann [50] (grey downwards triangle), Malin et al. [4] (grey square), and McEwen et al. [51] (grey cross), as given by Hartmann [50].Additionally shown is Teanby [30] (model 1, dotted line with grey uncertainty band).The 'Seismology' Gutenberg-Richter-derived rate is given by the red dashed line with red uncertainty derived by standard error propagation techniques (see Supplementary Section 3 for more details).The 'Cratering' results are given by the red diamonds with errorbars (computed as sqrt(N ) for each bin, as listed in Supplementary Table 1), with best-fit rate (solid red line).
1 Occurrence rate of VF events Since the wind noise recorded by InSight is prohibitive to almost all event detections during daytime in spring and summer, and to detections at all times of day during winter, a comparison of event rates between different times of Martian year, or between several years, requires to take the respective noise amplitudes into account: The catalogued events reflect only the detection rate of events, but not their occurrence rate.
We intend to answer if the data supports that the VF events occur at a constant rate, although our detection rate appears to vary strongly over time.Since any single realization of a stochastic process shows fluctuations compared to the expectation, we need to answer if the observed deviation from a stationary Poisson process is statistically significant, i.e. if it is so unlikely that it contradicts the assumption of a constant occurrence rate.
A methodology to assess the effect of the variable noise background was developed by Knapmeyer et al.
Supplementary Material Figure 3: Velocity spectrograms for events S0263a, S0421a, S0976b, and S1337a.Phase picks are marked by vertical dashed lines.These events have epicentral distances >5 °and are not currently associated with an impact crater.Note the different amplitudes for the left and right column.
[39] to demonstrate the seasonal variability of the occurrence rate of HF events.If the recorded amplitudes of events follow a power law with a slope b, and if the catalog is complete for events with recorded amplitudes exceeding a certain threshold amplitude A TH , then one can predict the number of events that would be missed if the background noise level increases: The signal-to-noise (SNR) ratio required for a certain detection dictates that A TH increases by the same factor as the noise level, and the number of events lost due to higher noise then depends on b.A detailed mathematical description of how to compute a short-term detection efficiency from noise RMS amplitudes, and how to average short time detection efficiencies to obtain daily or long-term average detection efficiencies, is given in section S4 of the supplement of Knapmeyer et al. [39].Said supplement also shows an application where the detection rate of HF events during the noisier second half of the night is predicted from the detections during the less noisy first half (and vice versa).
In the present work, we consider the 58 VF events of qualities A, B, and C that were detected up to Supplementary Material Figure 4: Velocity spectrograms for events S0981c, S0986c, S01034a, and S1160a.Phase picks are marked by vertical dashed lines.All four events have been identified as impacts [13,15].Note the different amplitudes for the left and right column.
May 11th 2022, or Sol 1228.The distributions of recorded amplitudes and SNR suggest that the catalog is representative for all events with amplitudes higher than -200 dB (displacement with respect to 1 m) and with SNR exceeding 2.9.Of the recorded VF events, 34 are thus representative for the actual VF event population (Since we know that the catalog is not complete due to the high noises at some times, it can only be representative).The cumulative number of detected events is shown as squares in Supplementary Material Figure 6a.When testing if the data makes a certain occurrence process of events plausible, we have to keep in mind that the actual event process, together with the observed noise conditions, resulted in the actual cumulative number of detections.Whatever we assume about the production of events must reproduce the observation at least during the time interval from which we derive the process parameters.
We will assume that events occur according to a stationary Poisson process.We will also assume that the Supplementary Material Figure 5: Spectral fits used to determine magnitudes for nearby meteoritic impacts S1034a and S1160a.The black dashed line marks the ambient noise before the event and is used to estimate the frequency window in which a fit can be done.The blue line marks the spectra in a 20 s window around the phase pick, manually selected to exclude wind or other obvious artefacts.The red line is a best fit of a Brune function, superposed with a Lorenz term for the local 2.4 Hz resonance [33] added to the pre-event noise.The long period flat part of this fit function minus noise is used to determine the event magnitude as described in Böse et al. [40].Note that the 'S-wave' window for S1160a is on the chirp of the event.
Supplementary Material Table 1: Binned incremental and cumulative crater counts computed from VF events, and the respective Area Time Factor (ATF) used for each bin.The final column gives the scaled cumulative values plotted in Fig. 4. only time dependent modification of detection probability is solely due to the background noise and not, for example, due to changes of the personnel in charge (the InSight team conducts a multi-stage review process to account for all kinds of operative changes).A statistically significant deviation between observed and predicted detections will point to either an incorrect rate estimation, or a non-stationary process.Absence or insignificance will not prove that our assumption is correct, but also not give a reason to doubt it.The procedure of our test is as follows: 1. Estimate the detection rate of events assuming a stationary Poisson process.The maximum likelihood estimation for the reference time window (RTW, 28 events from sol 288 to sol 1100) results in a detection rate of 0.03448 events per sol.The RTW covers the time between two major data gaps to a priori exclude possibly lost detections during the gaps.
2. Plot the observed cumulative number of detections (squares in Supplementary Material Fig. 6a) and the cumulative number expected from the rate estimate (dashed blue line in Supplementary Material Fig. 6a extrapolates beyond the limits of the RTW) 3. Compute the daily detection efficiency for the entire time as described by Knapmeyer et al. [39].The result is shown as solid red curve in Supplementary Material Figure 6a.Note the logarithmic scale.
4. Average the detection efficiency over the entire RTW.The resulting long-term average detection efficiency is 0.83227.
5. Estimate the occurrence rate of events assuming that the occurrence follows a stationary Poisson process: 0.03448/0.83227=0.04143events per sol, suggesting that 33 to 34 events actually occurred during the RTW.In the following, we assume that events occur according to a stationary Poisson process with this event rate.
6. Compute the expected cumulative number of detections from the stationary occurrence rate (0.04143) and the daily average detection efficiency (solid blue line in Supplementary Material Fig. 6a).
7. Plot as symbols in 6b the difference ("residual") between the observed cumulative number of detections (squares in Supplementary Material Fig. 6a) and the cumulative number predicted from stationary occurrence rate and daily detection efficiency (solid blue line in Supplementary Material Fig. 6a).
8. Numerically estimate by how much the predicted detection rate, estimated from the assumed stationary occurrence rate and accounting for the daily detection efficiency, scatters, i.e. estimate the distribution of the cumulative detection count at any time.To this end, generate a random sequence with the estimated occurence rate, and modulate it according to the time-dependent detection efficiency.From the synthetic cumulative detection count compute an apparent stationary rate, and the deviation of the synthetic cumulative detection count from the expectation resulting from the apparent rate.The gray lines in Supplementary Material Fig. 6b indicate the deviation from the expected value in terms of quantiles for commonly used probabilities (e.g. the 68.28 quantile indicates that 68.28% of the realizations result in values between the respective lines.Probabilities correspond to 1,2,3, and 4 standard deviations of a gaussian distribution, i.e. 0.6827, 0.9545, 0.9973, and 0.9999, respectively).9. Compare the deviation of the observed cumulative event count from a stationary Poisson process to the distribution obtained in step 8.To this end, we extract the probability to observe a deviation like the observed for the times of all events and plot it as function of time (Supplementary Material Fig. 7).We do not derive a single number as significance criterion but draw a conservative conclusion from visual inspection of the resulting graph.
We find that the deviation of the cumulative event count from the expectation for a stationary Poisson process reaches statistically significant amounts during the first martian year (sols 288 to 687) of InSight's operation when calibrating the stationary process using events from Sol 288 to Sol 1100.The actual probability of this deviation depends on the reference time window chosen, as the distribution becomes wide at long times before or after the RTW.Tests with alternative RTWs from Sol 288 to Sol 500, and from Sol 800 to Sol 1100, however suggest that the significance of the deviation is not an artefact of the window from sol 288 to 1100.A change of detection rate was also observed for the class of HF events [39] and, to a lesser extent, LF events [16,21], where the event rate apparently increased.The reason for this is unknown, but since the study of Dahmen et al. [21] is based on a machine learning technique, a change of MQS practices or personnel can be excluded as explanation of the detection rate increase.Supplementary Material Figure 7: Probability to observe a cumulative count as high as the actual observation at the times of the observations, obtained by reading probabilities p to not exceed a certain residual from Figure 6b and then plotting 1-p.Note that the observed residuals have a probability of less than 5% during much of the first martian year (i.e. up to Sol 687).
where M W is the moment magnitude, N the number of events, and a s and b s describe the productivity and slope of the distribution, respectively.Similarly recalling the impacts power law [Eq.2; e.g.49] where N is the number of impacts, D is the diameter, k the productivity, and b i the slope.

Rate uncertainty
The uncertainty of the impact rates is different between the two methods deployed in this study.For 'Seismology', three different parameters contribute to the final uncertainty: a s , b s , and c. a s and b s are derived from the four Gutenberg-Richter fits visible in Fig. 2b by calculating the mean and standard deviation for each parameter.These fits are obtained using different magnitudes of completeness and maximum distances for the area correction.a s is obtained from anchoring the fit to the event rate at the respective magnitude of completeness.This gives a s = 4.59 ± 0.16 and b s = 1.09 ± 0.06.Finally, c = (8.1 ± 3.0) × 10 8 (cf.Eq. 3, Fig. 3).Inserting expressions 2.8 into the impacts SFD N (≥ D) = kD −bi allows the derivation of the partial derivative of N with respect to a s , b s , and c.This can then be used to calculate the uncertainty of the form where σ are the uncertainty for the rate N , and parameters a s , b s , and c.Since a s and b s cannot be considered independent, this represents an indicative lower bound estimate.For the 'Cratering' method, the uncertainty in the crater size predicted for each event is determined by the uncertainties in seismic moment (σ M ) and the proportionality constant c (σ c ).These uncertainties are combined using standard error propagation methods.
The error bars in crater counts in each bin shown in Fig. 4 are given by √ N , which is a standard procedure for binned data.Eq. 2 is then fitted to the binned data using least squares algorithm in log space.This method produces uncertainties in both k and b i .These errors can be combined using standard error propagation to calculate the total uncertainty in the number of craters larger than 8 m across all of Mars, as follows: Because the value of N does not directly depend on the value of c, it is included in these calculations by rewriting Eq. 2 in terms of M and evaluating the partial derivative of N w.r.t.c.This results in uncertainty in N (≥ 8) quoted in Table 1 of ±99 craters.Results from both methods still agree with each other within the error margin.