The spatial flux of Earth’s meteorite falls found via Antarctic data

Extraterrestrial material that falls to Earth is heated, ablated, and can break up during its passage through the atmosphere. Any surviving rock fragments that reach Earth’s surface are termed meteorites and provide a tangible link to the evolution and composition of the solar system at different stages in its history (Joy et al., 2016; Schmitz et al., 2016; DeMeo, 2017; Heck et al., 2017;). Present estimates for the flux of extraterrestrial material falling to Earth’s surface rely upon either short-duration fireball monitoring networks (Halliday et al., 1996; Bland et al., 2012; Howie et al., 2017) or spatially limited ground-based meteorite searches from hot deserts (Bland et al., 1996; Gattacceca et al., 2011; Hutzler et al., 2016). Both approaches have their own advantages, and they each contribute to assessing potential hazards to life from larger impacts (Brown et al., 2013). However, a distinct benefit of the hot-desert-search approach is that retrieved samples enable us to quantify the transfer of solar system material to the Earth-Moon system across a wide range of sample masses without having to estimate the masses from fireball brightness.

Yet despite the benefits of controlled hot desert searches, the absolute number of searches and the number of meteorite samples collected within those searches are dwarfed by the number of controlled Antarctic searches that have taken place and the number of meteorites collected within them (MetBull, 2018). These Antarctic samples are collected from meteorite stranding zones (MSZs) (Folco et al., 2002; Harvey et al., 2014; Righter et al., 2014; Miao et al., 2018), which are blue ice areas typically located near mountainous regions of the continent. But until now, this major resource of near-pristine material has not been usable to accurately estimate the associated flux of extraterrestrial material to Earth. In part, this has been due to pairing uncertainties between meteorite samples, but it is due primarily to the complex mechanics by which meteorites are concentrated onto the surface of the MSZs (Zolensky, 1998; Bland, 2005; Zolensky et al., 2006): both from above by direct infall, and from below by the ablating ice flow releasing any entrained meteorites (Evatt et al., 2016).

To make use of the large number of collected Antarctic meteorite samples in order to estimate the flux of extraterrestrial material arriving on Earth, we constructed a mathematical model in stages and compared intermediate outputs to results from the literature. We first modeled the conservation of mass of a MSZ to determine the effective catchment area, which we then used to help determine the associated (glacial) flux of meteorites through a MSZ. With this, we were then able to determine a local extraterrestrial flux of meteorites. To compare with non-Antarctic flux estimates, we developed a model for capturing latitudinal variations of the flux, and demonstrate how well it fits to fireball observations. The resulting comparison of fluxes between geographically distinct regions also provides us with a data-driven approach to calculating Antarctic meteorite pairing factors; i.e., the number of paired-together meteorite fragments originating from the same fall. Finally, we discuss the implications of our model; in particular, how the risk to Earth from larger impact events has a strong latitudinal variation.

To overcome the barriers to estimating Earth’s present day fall flux via the untapped resource of Antarctic meteorite data, we analyzed glaciological and meteoritical data from 45 documented Antarctic MSZs, but focused upon the results from 13 systematically searched areas. These 13 areas were selected because they each have > 100 finds, a well-defined spatial extent, and together sample a range of Antarctic localities (Table 1; Fig. 1).

The first step was to determine the effective surface area of each MSZ, which is the surface area of the MSZ plus the upstream surface area of the ice sheet that feeds into it and then ablates. By assuming conservation of glacial mass (see the Supplemental Material¹), we find the effective surface area to be 2.58× larger than the actual area, on average. This implies that the areas of the upstream catchments are significantly smaller than previously assumed (Sinisalo and Moore, 2010), but consistent with a recent study (Zekollari et al., 2019). In addition, by considering the ice surface velocity for each MSZ (see the Supplemental Material), we find that meteorites have a mean surface residency time scale of ∼7.2 k.y. (neglecting the effect of wind, which removes lighter and smaller fragments much more quickly; Benoit and Sears, 1999; Folco et al., 2002). It is crucial to note that this surface residency time scale does not necessarily equate with a meteorite’s terrestrial age (i.e., the duration of time since that meteorite arrived on Earth), which is likely to be far greater for englacial samples (Kehrl et al., 2018).

Then, considering the rate of accumulation and loss of meteorites on the MSZ, and incorporating individual meteorite collection masses for each MSZ (MetBull, 2018) (some 13,200 stones in total), we can infer the associated meteorite flux rate above a minimum sample mass (see the Supplemental Material). In so doing, we find that the central estimate for the flux of meteorites >50 g to Antarctica is 45.3 [32.3–63.4] km⁻² m.y.⁻¹, where values in brackets denote the ±1 standard error.

Although meteorite fluxes are insightful, it is the fall flux (above a minimum total mass that reaches Earth’s surface) that is more commonly used. To convert our Antarctic meteorite fluxes into fall fluxes, we must take account of meteorite pairing. Pairing the large number of Antarctic meteorites into individual falls is impractical (Zolensky, 1998). To make progress, we initially used pairing estimates from the literature, of between two and six meteorites per fall, inferred from a limited sample of meteorites (Zolensky et al., 2006). In this study, we focused on minimum terminal fall masses of 50 g, irrespective of the pairing factor. In so doing, we found the central estimate of the flux of Antarctic falls >50 g to lie between 17.6 km⁻² m.y.⁻¹ (six meteorites per fall) and 32.2 km⁻² m.y.⁻¹ (two meteorites per fall).

Before we can further improve upon the above Antarctic fall flux estimates and compare them to those of other geographically distinct studies, we must consider the effect of latitude. This dependency was hypothesized by Halliday (1964) but has not been verified observationally until now. By plotting the normalized frequency of satellite fireball observations (NASA CNEOS, 2018) against latitude (Fig. 2), we observe a clear latitudinal variation: the fireball frequency per unit area decreases as the latitude rises, to ∼65% of the equatorial frequency at the poles. To model this latitudinal variation in fireball frequencies, we considered the trajectories of meteoroids emanating from the ecliptic plane under the influence of gravity, where daily and annual rotations of Earth cause longitudinal variations to average out (see the Supplemental Material; Halliday, 1964; Le Feuvre and Wieczorek, 2008).

The modeled fireball frequencies (solid line in Fig. 2) show a reasonable fit to the observations. Comparing our method to the assumption of a uniform model (see the Supplemental Material), we found that our latitudinal-variation model performs better than the uniform case (dashed line in Fig. 2) at the 95% confidence level (p = 0.012). Due to uncertainty in past satellite network coverage (Brown et al., 2002), we repeated this test on fireball data from the past 10 yr, and again found confidence at the 95% level. We also compared our result to an extant frequency-versus-latitude curve (Le Feuvre and Wieczorek, 2008) (dotted line in Fig. 2), calculated using a different methodology (aimed at improving impact-crater counting studies of solar system bodies). Again, we found that our latitudinal-variation model produces a better fit to the fireball data than the alternative at the 95% confidence level (p = 0.033).

Accounting for latitude variation is clearly essential for fairly comparing estimates from globally distinct meteorite accumulation zones. Using the result shown in Figure 2, it is now a simple matter to correct local flux estimates to their equatorial equivalents according to their latitude (see the Supplemental Material). This done, we calculated an effective pairing factor by carrying out a least-squares fit between our equatorial-equivalent fall fluxes and those obtained from previous studies (Bland et al., 1996; Halliday et al., 1996) over a minimum mass range of 10 g to 1 kg. A pairing factor of 3.18 minimized the residuals. This approach allowed us to circumvent the aforementioned issue of experimentally pairing the large number of Antarctic meteorites, and produces a metric that accounts for the sometimes subjective field-collection issue of whether or not two “nearby” samples should be considered part of the same fall. The Antarctic pairing factor is comfortably within the independently estimated range of 2 to 6, producing a local Antarctic fall flux >50 g of 25.8 km⁻² m.y.⁻¹, while the equatorial-equivalent fall flux >50 g is 38.7 [26.8, 55.9] km⁻² m.y.⁻¹. Likewise, the resulting fall flux fits comfortably within the equatorial-equivalent flux values from our key comparative studies (Bland et al., 1996; Halliday et al., 1996), which are 31.5 km⁻² m.y.⁻¹ and 51.0 km⁻² m.y.⁻¹, respectively. Although not directly comparable due to the 2 m.y. time period analyzed, Drouard et al. (2019)’s meteorite flux value >10 g of 234 ± 16 km⁻² m.y.⁻¹ lies close to the upper standard error limit of our estimate (149 [101–219] km⁻² m.y.⁻¹) when both are converted to equatorial equivalents, potentially due to wind effects.

Our equatorial-equivalent Antarctic-derived fall flux as a function of minimum mass is shown in Figure 3. The results show a well-established logarithmic decay in the number of falls with minimum mass (Huss, 1990), with a clear flattening below ∼20 g (i.e., consistent with samples lost to wind removal; Schutt et al., 1986; Folco et al., 2002). Also plotted are the equator-equivalent fall fluxes from the literature (Bland et al., 1996; Halliday et al., 1996). It is evident that all of these flux estimates are reassuringly close to one another. However, using a larger number of collection sites allows us to move beyond qualitative error estimates (Halliday et al., 1996) or simple ranges defined by two or three sites (Bland et al., 1996), and instead provide standard errors. Specifically, the best modern estimates were still subject to error factors of 2 to 3, while earlier studies differed by a factor of 7 (accumulation studies) or an order of magnitude or more (eyewitness falls) (Bland et al., 1996). Thus, ours appears to be a substantially more robust estimate than previous studies have been able to produce, particularly when weight is given to the fact that our estimates rely upon a much larger fraction of the global meteorite collection (22% and 13,200 samples versus ∼250 events) (see the Supplemental Material).

Accounting for the latitudinal variation allowed us to expand our Antarctic results to a global setting. Using the best-fit pairing estimate, we determined a global expected fall flux of 17,600 [12,200–25,400] yr⁻¹ for masses >50 g. Going further, we calculated the expected recoverable mass flux; that is, the mass flux derived by integrating the equatorial-equivalent Antarctic flux (Fig. 3) between the lowest observed collection mass and the highest expected collection mass during a systematic search (see the Supplemental Material; Huss, 1990). We estimate this equatorial-equivalent mass flux to be 36.6 [25.1, 53.4] kg km⁻² m.y.⁻¹, a value which is 31% lower than an equatorial-equivalent estimate of 52.8 kg km⁻² m.y.⁻¹ from the combined results of Halliday et al. (1996) and Bland et al. (1996). The corresponding total mass flux to Earth is 16,600 kg yr⁻¹. Interestingly, by focusing on the integrated mass flux between 10 g and 1 kg, we found equatorial-equivalent values to be 8.11 [5.69, 11.6] kg km⁻² m.y.⁻¹, while the corresponding total mass flux to Earth is 3680 [2580, 5250] kg yr⁻¹—a value that is in agreement with previous estimates and yet better constrained (4380 kg yr⁻¹, Halliday et al., 1996; 2840–7150 kg yr⁻¹, Bland et al., 1996). This divergence at masses >1 kg may indicate a lower flux in the 1–100 kg range than previously modeled (Bland and Artemieva, 2006), or a limitation in determining mass flux of larger falls from Antarctic MSZs where the fall return period is greater than the residency time.

The emphasis in our work has been to use available glaciological and meteoritical data sets to calculate the extraterrestrial flux rate. However, the core relation (see the Supplemental Material, Equation S4) can be inverted to predict the expected meteorite number density of a potential MSZ when given the extraterrestrial flux rate. In fact, even without a flux rate estimate, a ranking of potential MSZ number densities can be obtained. Given the support required to operate remotely in Antarctica, to be able to predict new areas of the continent that are likely to hold significant accumulations of meteorites is extremely beneficial. Indeed this was the original purpose of our model. In planning the first United Kingdom–led Antarctica meteorite collection mission, we had predicted that an unvisited blue ice area near the Shackleton Range (East Antarctica) to have a meteorite number density of 10.8 [7.5, 15.6] km⁻². This mission to the Outer Recovery Ice Fields area was carried out during austral summer 2018–2019, finding a total (unclassified) meteorite number density of 7.1 km⁻², close to the lower 1 standard error estimate (Joy et al., 2019).

The latitude component of our modeling also has applications and implications beyond our immediate study. For example, it is possible to apply our latitude model to crater-counting studies on other planets. Yet at a societal level, it is the latitudinal variation of impact frequencies on Earth that likely carries the most importance—principally due to the independence of these results on meteorite mass. Consequently, our latitude model can be used as a likelihood weighting, which would help to geographically quantify threats and reassess flux estimates for large impactors, calculations that have predominantly been limited to analysis of the North American and European cratons. Finally, this work shows that equatorial regions face an appreciably greater risk of larger impactors—an increase of 12% relative to uniform latitudinal variation assumption—while at higher latitudes, the risk is reduced by 27%. As such, it may give additional reason for placing long-term contingency facilities at higher latitudes, such as the Global Seed Vault in Svalbard, Norway.

This work was supported by The Leverhulme Trust (UK) grant RPG-2016–349; Royal Society (UK) grant RS/UF140190; Science and Technology Facilities Council (STFC) grant ST/R000751/1; and Engineering and Physical Sciences Research Council (EPSRC) grant EP/R014604/1. We are grateful for useful discussions at the Royal Astronomical Society (UK) meeting held in November 2018 and at the European Geosciences Union General Assembly 2019. We are grateful to Vinciane Debaille and an anonymous reviewer for the time they have taken to review the manuscript and for their constructive suggestions. We also thank Benjamin Tipping and Philip Cartwright, who both assisted with preliminary work on the latitudinal aspect of this study.