Determination of strewn fields for meteorite falls

ABSTRACT

When an object enters the atmosphere it may be detected as a meteor. A bright meteor, called a fireball, may be a sign of a meteorite fall. Instrumentally observed meteorite falls provide unique opportunities to recover and analyse unweathered planetary samples supplemented with the knowledge on the Solar system orbit they had. To recover a meteorite from a fireball event, it is essential that recovery teams can be directed to a well-defined search area. Until recently, simulations showing the realistic mapping of a strewn field were difficult, in particular due to the large number of unknowns not directly retrieved from the fireball observations. These unknowns include the number of fragments and their aerodynamic properties, for which the masses of the fragments need to be assumed in a traditional approach. Here, we describe a new Monte Carlo model, which has already successfully assisted in several meteorite recoveries. The model is the first of its kind as it provides an adequate representation of the processes occurring during the luminous trajectory coupled together with the dark flight. In particular, the model comprises a novel approach to fragmentation modelling that leads to a realistic fragment mass distribution on the ground. We present strewn field simulations for the well-documented Košice and Neuschwanstein meteorite falls, which demonstrate good matches to the observations. We foresee that our model can be used to revise the flux of extra-terrestrial matter onto the Earth, as it provides a possibility of estimating the terminal mass of meteorite fragments reaching the ground.

methods: analytical, methods: data analysis, Earth, meteorites, meteors, meteoroids, planets and satellites: atmospheres

1 INTRODUCTION

Meteorite falls with well-documented atmospheric trajectories are of paramount importance for planetary science as they yield tangible planetary samples supported by comprehensive spatial context. With appropriate registration of a fireball entry, the pre-impact orbit of the meteoroid in the Solar system can be calculated and a possible connection to its parent body may be established (Dmitriev, Lupovka & Gritsevich 2015).

The count of currently known asteroids has passed the 1 million mark (NASA 2020; MPC 2021). It has been estimated that the Solar system has over 150 million asteroids larger than 100 m and a countless amount of smaller ones. The Meteoritical Bulletin Database holds records of 65 184 valid meteorite names (MBD 2021) including those, which ‘sample’ the parent objects that do not longer exist in the Solar system in their original form (Vinković & Gritsevich 2020). Comparing the number of known meteorites to the numbers of known asteroids implies that at present our sampling of the Solar system objects is very limited.

When an object enters the Earth's atmosphere, it produces a luminous event called a meteor (Silber et al. 2018; Vinković & Gritsevich 2020). We illustrate the phenomena and associated processes in Fig. 1. The objective of this study is to present a first of its kind model that encompasses a description of the processes that have an influence on the resulting strewn field. Naturally, these processes occur during both, the luminous and dark flight.

Figure 1.

Different stages and phenomena related to a meteoroid entering the Earth's atmosphere. The pre-atmospheric (or initial) mass of the meteoroid is the mass that a meteoroid possesses on its orbit before it starts to interact with the Earth's atmosphere. The starting mass of the meteoroid is defined here as the mass that it has at the starting point of our MC simulation, following initial ablation but prior to any fragmentation to occur. The terminal mass is the cumulative mass of meteorite fragments that land on the Earth's surface.

The trajectory of a fireball (an exceptionally bright meteor) can be resolved from photographic or video recordings. It is a geometrical triangulation task with difficulties of its own (Ceplecha 1987; Borovička et al. 2013; Lyytinen & Gritsevich 2013, 2016a; Spurný et al. 2014; Sansom et al. 2019a; Jansen-Sturgeon et al. 2020; Peña-Asensio et al. 2021). The next question to answer based on the analysis of the fireball is deciding its probable fate and if any part of the object has survived the atmospheric entry (Gritsevich, Stulov & Turchak 2012; Turchak & Gritsevich 2014; Sansom et al. 2019b; Moreno-Ibáñez et al. 2020). Most meteor events undergo complete ablation, hence in most cases nothing or only a negligible amount of matter reaches the ground. For a prominent meteorite fall, the next question arises – where can it be found?

It is possible to solve the dark flight path of a meteorite numerically, assuming a single projectile and using the data derived from the fireball trajectory supplemented with atmospheric data. Such calculations produce a single impact point on the ground. By repeating the calculation for different assumed masses, a polyline can be drawn on the map (Spurný et al. 2012; Brown et al. 2019). This is a nominal or central line of the strewn field, sometimes called the highest probability line (Spurný et al. 2014). This line is traditionally taken as a good place to start any meteorite recovery efforts.

However, the central line of the strewn field derived this way is not the highest probability line to find a meteorite fragment. Fragmentation of the projectile, shapes of individual fragments, and actual atmospheric winds during the dark flight distribute fragments over a wide area. This creates a need for the generation of probability heat maps (Devillepoix et al. 2018). Another (practical) problem is that this central line does not always fall in favourable terrain for the search of meteorites.

Monte Carlo (MC) methods are used when complex physical phenomena with several degrees of freedom need to be simulated. MC simulations became a widely used approach to solve different kinds of complex problems in physics and have also been proposed to simulate the dark flight of meteorites.

The main argument towards using the MC methods for strewn field simulations is that the exact atmospheric condition profile for the fall trajectory of each individual fragment is not possible to obtain. In addition to that simulations should also account for great variations in shapes and sizes of the fragments.

Our first dark flight Monte Carlo (DFMC) simulations were done in 2010 January. The model was tested with several falls where the trajectory of a fireball and distribution of meteorites were known. The most illustrative case study was the Košice meteorite fall that occurred on 2010 February 28. This well-documented fall produced data of 218 different sized fragments from 0.3 g up to 2.37 kg in addition to the fragments recovered by private meteorite hunters (Borovička et al. 2013; Gritsevich et al. 2014c).

The first DFMC simulation leading to the successful field campaign for a new fall was made in 2014 when two pieces of the Annama meteorite were found inside the predicted strewn field in Russia (Gritsevich et al. 2014a, b).

In this study, we describe how to determine a realistic strewn field map using the MC methods and the real corresponding atmospheric measurements for the time, heights, and location of the fall. As in our earlier study of the luminous part of the trajectory (Lyytinen & Gritsevich 2016b), the atmospheric data comes from the GFS, Global Forecast System, and from the ECMWF, the European Centre for Medium-Range Weather Forecasts, and is provided by national weather services. The performance of the DFMC model is demonstrated using the example of the Košice meteorite fall by comparing the simulated strewn field with the actual positions of the recovered fragments. Possible implications of the model are then discussed by also studying the example of the Neuschwanstein meteorite fall.

2 METHODS

2.1 General considerations

MC simulations are based on repeating calculations multiple times using phenomena related equations and rules with variations in initial parameters. Plausible variations in parameters are produced using appropriate distributions. In our model, two probability distributions are used. Most variations are achieved by using the normal distribution (equation 1) as a probability function and some – by using a variation of the exponential distribution (equation 2), which for continuous random variable |$\chi $| can be defined as

$$\begin{equation*} f\left( x \right) = \frac{1}{{\sigma \sqrt {\left( {2\pi } \right)} }}\,\,{\mathrm{ e}^{ - \frac{1}{2}{{\left( {\frac{{\chi - \mu }}{\sigma }} \right)}^2}}}, \end{equation*}$$

(1)

$$\begin{equation*} f\left( x \right) = {\left( {1 - {\rm{\mu }}} \right)^i}{\rm{\,\,}}\left( {1 - \chi } \right) + \chi . \end{equation*}$$

(2)

The parameters of the normal distribution are the mean, or the expectation of the distribution (μ), which determines the location of the peak, and the standard deviation (σ), which indicates the spread of the values around the mean. In equation (2), (μ) is the location parameter, which determines the ‘location’ or shift of the distribution. The shape of the probability curve is controlled by the rate parameter (ί).

MC methods can also be used to simplify or speed up the computations. Modelling a dark flight of a supersonic projectile can be a time consuming task. The most difficult part in dark flight simulations are aerodynamic forces often involving factors like Reynolds number (Vinnikov, Gritsevich & Turchak 2016). These forces depend on the fragment shape, orientation, and character of motion. Subsequently, forces influenced by these properties are unique for each individual fragment. For example, the Magnus effect occurs when an airborne body rotates. Yet, not all meteorite fragments rotate so fast that they experience it.

Aerodynamic forces, which cannot be ignored while simulating a dark flight of a meteorite, are drag and atmospheric winds. Aerodynamic drag is a force, which slows down an object moving through the atmosphere. This force appears due to the air resistance and acts opposite to the relative motion of the object with respect to the surrounding media. If the actual shape and orientation of an object are ignored in the simulation, air drag does not modify the direction of the trajectory. Ignoring some of these more random factors affecting dark flight is an MC technique known as importance sampling.

To produce plausible DFMC simulations, there has to be a variation in fragment shapes to represent different drag forces induced by differently shaped meteorite fragments. This can be implemented by using different values for the drag coefficient and ignoring forces which may deviate the fragment from the original trajectory.

Our DFMC simulation model is based on a programme that calculates a ballistic free fall trajectory. The programme accounts for gravity, air drag, curvature of the Earth and can be further extended to also account for the Coriolis force. A new spatial position and the forces acting on a simulated fragment are recalculated after every time-step. Our time-step interval is 0.04 s due to the video frame rate of some fireball cameras. The free fall simulation is then enhanced by adding the atmospheric winds and MC variations. A horizontal force caused by winds is calculated as air drag using equation (4). Also, the mass-loss and change in cross-section area are recalculated after every time-step as described in Section 2.4.

Inputs for our DFMC simulation are parameters of the meteor trajectory (Table 1) and atmospheric data (Table 2).

Table 1.

Parameters of the start point and the lowest observed height of the fireball are used as input data for a DFMC simulation. In the parameter file for our code, there are several other inputs, like time code and limit to calculated trajectories, but they are part of the programme execution and do not belong to the algorithm of the simulation.

Symbols	Parameter	Unit	Remarks	Default error
λ0	Longitude	degrees	WGS84 (west < 0° < east)	Included in e0
φ0	Latitude	degrees	WGS84 (south < 0° < north)	Included in e0
h0	Altitude	m	Above the sea level	Included in e0
e0	Spatial error at the start point	m	An uncertainty sphere	±300 m
δ0	Direction of trajectory	degrees	0°–360° (0° = N, clockwise)	±1|${_{.}^{\circ}}$|8
γ0	Trajectory slope	degrees	0°–90° (90° = vertical)	±1|${_{.}^{\circ}}$|8
V0	Velocity	m  s−1	Velocity at the start point	±300 m s−1
a0	Deceleration limit	m  s−2	Defines acceptable deceleration at the start point	a ≥ 0.99a0
he	End height of fireball	m	Lowest observed altitude	No variations
ρm	Density of meteoroid	kg  m−3	For chondrites 3300 kg m−3	±500 kg  m−3

Table 2.

Parameters of the atmospheric data used in the DFMC model. Units are given as used in the calculations. Atmospheric data uploaded from a weather data base may also include other parameters like a dew point temperature.

Symbols	Parameter	Units	Remarks
ha	Altitude	m	Above the sea level
ta	Temperature	°C	(May be K in data)
pa	Air pressure	kPa	(May be hPa in data)
ρa	Air density	kg m−3	(Can be calculated)
δw	Wind direction	degrees	0°–360° (0° = N, clockwise)
vw	Wind speed	m s−1	(May be knots in data)

We do not go into details of the DFMC code in this work, though a flowchart of a typical execution is shown in Fig. 2.

Figure 2.

A flowchart of the DFMC simulation code presents the main stages of the code as they are usually executed. The first fragment (or the starting mass) is defined as the nominal body (see Section 2.3) for which the trajectory code execution is presented on the grey background. No MC variations are applied for the nominal body case. After the calculation of the nominal body trajectory, a trajectory of a fragment is simulated in the main loop. The main loop ends when the total mass of simulated fragments at the start point exceeds the starting mass or the maximum number of fragments allowed is reached. Letters in the flowchart are Y for yes and N for no.

2.2 Start point

The start point is the point where a simulation run begins. It is advisable to set the start point on a luminous trajectory point where the most accurate parameters can be calculated from the fireball observations. This is almost never the point of maximum fireball brightness, because images are often overexposed so that the exact projectile location cannot be determined accurately enough. The starting point should also not be too close to the end of a luminous flight because fragments that separated from the main body earlier will not be simulated.

Higher precision of the luminous trajectory reconstruction also implies higher accuracy in the DFMC model. Spurný et al. (2014) estimated that they achieve 12 m accuracy (one standard deviation) for the Benešov fireball trajectory. However, the study by Sansom et al. (2019a) revealed that a meteoroid's trajectory during luminous flight does not fit a straight line. It can deviate from a straight line by hundreds of meters or even up to 3 km.

Parameters for the start point of a simulation are taken from the fireball observations and MC variations used for most of the start point parameters are thus error margins of the trajectory calculations. Using a normal distribution (equation 1) as a probability function, the start point values for the location and flight direction for each simulated fragment can be calculated. The calculated precise location and direction for the trajectory gives the average absolute deviation or the mean deviation (µ). Error margins can be added as a standard deviation (σ) of the distribution. If error margins are not known for the start position, some default values can be used as error margins. In Table 1, we summarize the default errors used in our DFMC simulations.

Fireball trajectory calculations usually give error margins for the start point coordinates and altitude. These errors can be notable, especially in the case when observation sites are unfavourably located with respect to the fireball trajectory geometry. Naturally these error margins can be used as MC deviations for the start point parameters. However, that will cause a skewed box-like space where possible start points of simulated trajectories are located. That can cause artefacts to the strewn field prediction and that is why we combined altitude and coordinate error margins into one spatial error parameter (e₀) of the start point position (Table 1).

MC variations for the start position and direction for a simulated fragment are introduced by applying the MC algorithm as follows. First, a random number is generated. Most random number generators return a number ranging between 0 and 1. This number is set as our test value. Next, a random number is generated for the parameter constrained by the range defined by the parameter value as a mean value and the error margins of the parameter as a deviation. If the value of the normal distribution function (equation 1) for this parameter value exceeds or equals our test value, we have found a value of the parameter for this simulated fragment.

If the test value is greater than the value of normal distribution function, a new test value and parameter value will be generated and the test will be repeated. This loop will be repeated until the test is passed. The efficiency of the code can be improved by normalizing the mean value of the normal distribution function to be 1. In this case every time when the random parameter value is the same as the mean value the test will be passed.

2.3 Mass of a meteoroid at the start point

The main question related to a meteorite fall is how much meteoritic material one can expect to find on the ground. Here, we will demonstrate that this question can be answered with the DFMC simulation when MC methods are limited to simulate only a single trajectory instead of distribution map of probabilities for multiple fragment trajectories.

How do we know the total mass of the meteoroid at the start point of simulation? Existing trajectory reconstruction methods allow us to reproduce the preatmospheric masses of meteoroids relatively well, based on fireball aerodynamics, so that these agree with the diameter estimates obtained in the laboratory based on cosmogenic radionuclide analysis done on recovered meteorites (Gritsevich 2008a; Gritsevich et al. 2017; Kohout et al. 2017; Meier et al. 2017). Pre-atmospheric size (defined as the physical radius × bulk density, with a unit of g cm⁻²) can then be translated to the preatmospheric meteoroid mass involving assumptions about the bulk density and initial meteoroid shape. Further calculations of mass of the main body along the luminous trajectory become more challenging as they require the knowledge of the mass loss parameter (or ablation coefficient) as well as the knowledge or assumption of the shape change coefficient, which generally ranges between 0 and 2/3 (Gritsevich & Stulov 2006; Gritsevich 2008c; Bouquet et al. 2014; Drolshagen et al. 2020). The mass of a meteoroid can also be estimated using Bayesian filtering techniques (Sansom, Rutten & Bland 2017). In practice, mass estimates will change even when varying assumptions on the atmospheric model (Lyytinen & Gritsevich 2016b).

Alternatively the meteoroid mass can be resolved using photometry, by analysing which fraction of the kinetic energy is transformed into the light emitted by the fireball (McCrosky 1967). However, determination of the photometric mass this way roots into the assumption of constant velocity (Gritsevich & Koschny 2011) and relies on poorly assumed luminous efficiency coefficient (Gritsevich 2008b), making it barely acceptable for fireballs producing meteorites.

The starting mass, defined here as the total mass of a meteoroid at the start point of a DFMC simulation, can be solved using observed velocity and deceleration. The mass of a meteoroid can be calculated when its shape and bulk density are fixed. We ‘fix’ the shape at the start point by fixing the drag coefficient of the meteoroid to be 1.5 and, as a default, we set bulk density of a typical ordinary chondrite as 3300 kg m⁻³.

The drag coefficient (C_d) reflects the dependence of the drag force on shape, orientation, and air flow conditions. It is calculated by equation (3), where A is the shape factor, V is velocity in m s⁻¹, and t_a is the ambient temperature in °C corresponding to the altitude. The specific heat ratio γ is defined as 1.4 for air at STP (T = 273.15 K or 0 °C, P = 101.325 kPa), the gas constant R is 286 m² s⁻² K⁻¹, and the constant T = 273.15 K is used to convert the ambient temperature from the Celsius to Kelvin scale.

$$\begin{equation*} {C}_{\text{d}} = {{AV}}/(\gamma {R}({T} + {t}_{\text{a}}))^{1/2}. \end{equation*}$$

(3)

The cross-section area (S) is estimated by calculating the characteristic radius corresponding to the mass. This is needed for the air drag equation (equation 4) to relate the mass to the produced deceleration (equation 5).

$$\begin{equation*} {{{F}}_{\rm{d}}} = {\rm{ }}\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} {\rm{ }}{{{C}}_{\rm{d}}}{{\rm{\rho }}_{\rm{a}}}{{S }}{{{V}}^2}, \end{equation*}$$

(4)

$$\begin{equation*} \frac{{\mathrm{ d}V}}{{\mathrm{ d}t}}\,\, = \frac{{{F}}}{M}.\,\, \end{equation*}$$

(5)

Hence, the starting mass is the mass that matches the deceleration observed at the start point.

This starting mass is used as the total cumulative mass of fragments at the start point. We refer to this virtual, non-fragmenting meteoroid as the nominal body. The nominal body at the start point should have the slowest observed deceleration (thus also highest velocity) what is calculated from the fireball observations.

In order to cross-check that the simulation is plausible, we calculate the full bright and dark flight trajectory for the nominal body. The trajectory of the nominal body is calculated without any variations in wind speed or in wind direction, so the nominal body lands on the central line of the strewn field. In this computation, only mass-loss occurring due to ablation will apply to the nominal body during the passage through the atmosphere. This optional computation reveals how much mass can be ablated during the atmospheric passage. The nominal body is not allowed to break up, so the terminal mass estimate resulting from this trial is usually significantly bigger than any meteorite fragment resulting from the complete simulation loop with the start point parameters described in Section 2.2.

2.4 Ablation

The start point of the simulation is on the luminous path of the fireball, as this is also the observable part of the trajectory. This implies that the meteoroid is still experiencing mass-loss, in particular due to ablation. As fireballs are relatively rapid phenomena, even shortly after the fireball terminates, surviving fragments still undergo minor ablation (Fig.   1). The velocity of the projectile is high at this point and deceleration is high. Ablation will end in seconds after the fireball fades away (given there are surviving fragments). The term ablation also usually includes ablation due to spallation, but spallation is considered here as a part of the fragmentation process.

How fast the meteoroid loses its mass due to ablation can be calculated using the ablation coefficient. The ablation coefficient of a meteoroid is difficult to observe directly but it can be obtained from the fireball observations (Gritsevich 2009; Bouquet et al. 2014). Using the same ablation coefficient with meteoroids of different bulk densities, for example, for iron and carbonaceous chondrite, simulations give about the same strewn field results, but sizes and therefore masses of the simulated fragments will be different. For obtaining a more realistic mass distribution, the ablation coefficient should be corrected to account for the meteoritic material. The problem is that the meteoroid material is not known for sure, before recovering any of its fragments on the ground. The use of chondritic material as a default for the simulations is a good choice since statistically over 80 per cent of the known meteorite falls are chondrites. Below the start point, the default ablation coefficient for the chondritic material in our code is 0.0018 s² km⁻², though it is more the lower bound for chondritic material. In theory, the ablation coefficient σ can also be related to the mass-loss parameter β as σ  = 2β[(1 − μ)V_e²]⁻¹. However, the mass-loss parameter β when determined based on fireball observations is often used to account for all the mass-loss mechanisms accompanying the phenomenon (including any fragmentation), while ablation and fragmentation are modelled independently in the present work.

To calculate mass-loss during the ablation stage, we use equation 6 described in Gritsevich (2009):

$$\begin{equation*} {{M = }}{{{M}}_{\rm{e}}}{\rm{exp ( - 3\beta }}( {{{1 - }}{{( {{{V / }}{{{V}}_{\rm{e}}}})}^{\rm{2}}}} ){\rm{)}}. \end{equation*}$$

(6)

Here, M_e and V_e stand for the pre-atmospheric mass and velocity of the meteoroid, respectively, and β is the mass-loss parameter defined from the analysis of the luminous segment of the trajectory (Lyytinen & Gritsevich 2016b). In equation (6), we have adopted the value of the shape change coefficient μ = 2/3, because this assumption is consistent with ablation occurring evenly over the fragment surface. A value of μ = 2/3 was also found to be the most common result based on the analysis of the MORP fireballs by Bouquet et al. (2014) and it is nearly so for the FRIPON meteors (Drolshagen et al. 2020). That also corresponds to the assumption that the fragments have a spherical shape and they are tumbling randomly during the atmospheric passage.

By applying equation (6) twice: (1) for the (current) point with mass M and velocity V, and (2) for the intermediate (reference) point with mass M₀ and velocity V₀, we arrive at the following:

$$\begin{equation*} {{M}}={{\rm{M}}_{\rm{0}}}{\rm{exp ( - 3\beta }}{( {{{{V}}_{\rm{0}}}{{/}}{{{V}}_{\rm{e}}}})^{\rm{2}}}( {{\rm{1 - }}{{( {{{V/}}{{{V}}_{{0}}}})}^{{2}}}}){{)}}. \end{equation*}$$

(7)

As demonstrated in Fig. 1, the end of the luminous flight (termination of the fireball) precedes the end of ablation. Even after the luminous flight ends, the temperature at the stagnation point exceeds melting temperatures of most minerals found in meteorites. For example, the melting temperature range of olivine is from 1478 to 2163 K depending on the composition and is similar for plagioclase feldspar, but some minerals melt at even lower temperatures (Wenk & Bulakh 2011). We have chosen a temperature of 1800 K as the end of ablation threshold.

The final dark flight stage starts when ablation of a meteorite fragment ends (Fig. 1). After that point, the fragment size and shape do not change unless it breaks up due to thermal shock or centrifugal force induced by fast rotation.

2.5 Transition to subsonic flight

Subsonic flight refers to local aerodynamic flow conditions Ma < 1, where Mach number (Ma) is defined as the ratio of the local flow velocity to the local speed of sound. Aerodynamic flow conditions when Ma > 1 correspond to supersonic flight. In practice, the greatest degree of aerodynamic unpredictability is associated with a Mach number range 0.8 < Ma < 1.2 called transonic flight conditions (Cook 2013).

When ablation ends and a fragment enters into the dark flight stage it is still supersonic for a short period. Transition from supersonic flight to subsonic flight has an effect on the drag coefficient. When a fragment is slowing down below Mach 5 the drag coefficient begins to increase slowly. It reaches its maximum at Mach 1.5–1.8. When slowing down from there, the drag coefficient will decrease steadily so that around Mach 1 the drag coefficient is about the same as above Mach 5. The drag coefficient changes by a factor of 2 between Mach 0 and Mach 1.5. We have applied a change of drag coefficient in our DFMC code following the best-fitting function for a sphere by Carter, Jandir & Kress (2009). In our model, the drag coefficient has a maximum of 1. When Ma > 5 the drag coefficient is 0.92 and gets as low as 0.43 when Ma < 0.1.

2.6 Wind drift

The shape and location of a meteorite strewn field are the result of the atmospheric winds occurring during the luminous flight of the fireball and the dark flight free fall stage of the meteorite in the lower part of the atmosphere. Wind is the effective horizontal force during the trajectory and even more so at the lower altitudes during the free fall stage of the dark flight.

The wind speed and direction are incorporated into the DFMC model using the atmospheric data that can be randomized using the normal distribution function (equation 1). But how much does the wind speed and direction vary? An atmospheric data profile gives only one averaged value for each altitude, however wind data vary between different radiosonde measurements.

To understand the effect of wind drift on the width of the resulting strewn field, we have studied the maps of known meteorite strewn fields (manuscript in preparation). After studying 29 different meteorite strewn field maps the following conclusions have been reached:

• In many cases, it is not clear if all meteorite fragments have been recovered or was the full extent of the (actual) strewn field even searched for. For this reason, only the strewn fields with over 10 individual meteorite fragments were selected for the study and average width of their strewn field was measured.

• Many of the narrow, less than 2 km wide, strewn fields were produced by iron meteorites like Sikhote Alin. The mean width of the strewn fields of stony meteorites are 1–7 km and over half of them are 2–4 km wide. On average, the mean width of a stony meteorite strewn field is 3.8 km.

• The length of the strewn field is mainly a result of the entry angle of the meteoroid and wind direction. In most cases, meteorite fragments have spread more or less sporadically and any trace of highest probability line mentioned by Spurný et al. (2014) inside the strewn field is difficult to see.

• In many strewn fields, clusters of meteorite finds exist. They are likely the results of fragmentation events during the trajectory.

To produce strewn fields with a width comparable to the known meteorite strewn fields, variations of wind speed have to be at least ±10 per cent and probably closer to ±20 per cent. In our models, wind speed deviation is randomized for each simulated fragment so that the wind speed deviation for the fragment will stay the same from the start point to the ground. If variations are made for each position change of the falling fragment, random wind speed changes tend to cancel each other out.

Variation to wind direction has been implemented in a similar way. However, it is even more uncertain how much wind direction can vary from radiosonde data and during a free fall trajectory. We use ±10° or ±20° variation in wind direction. Variation in wind directions causes the redistribution of fragments’ masses inside the strewn field. Variations of wind speed and wind direction used are about the same as breakthrough uncertainty requirements for radiosonde measurements by World Meteorological Organization (Nash 2015), although we use a per cent value instead of meters per second for the wind speed.

Horizontal components of the winds are applied to the trajectory as a horizontal force by using equations (4) and (5). The atmospheric data or the profile of the atmospheric conditions for the meteorite fall location and time can be generated from a weather observation data base or radiosonde readings from nearby stations can be used. However, usually vertical resolution of radiosonde data is not very good and often the highest altitude of soundings are below the starting point of a DFMC simulation. Also, if the nearest weather station is too far from the fall location, atmospheric data may be simply incorrect for the fall site. In many cases, the radiosonde data and the atmospheric data generated from the weather data base have given us almost the same result. Subsequently, the use of radiosonde data is far better than not accounting for the real atmospheric data at all.

We use an atmospheric profile forecasts corresponding to the meteorite fall time and coordinates from the GFS data base. Vertical resolution of our atmospheric data is 200 m. Parameters are listed in Table 2. More sophisticated 3D model of the atmosphere using the Weather Research and Forecasting (WRF) software with the Advanced Research WRF dynamic solver has been used by the Desert Fireball Network team (Devillepoix et al. 2018). Large variations in wind speed and direction which our DFMC applies to the wind data hides possible benefits of using the 3D atmospheric model. When predicting only the central line of the strewn field then the use of the 3D model of the atmosphere may enhance precision of the result.

2.7 Account for fragmentation

Most of the small meteoroids entering the Earth's atmosphere will be vapourized. Larger meteoroids, which do not meet this fate, will be eroded by heat, subsequent ablation, and spallation. Large bodies may also break up into smaller fragments. Disintegration of a meteor body takes place when dynamic pressure exceeds the compression strength of the meteoroid (e.g. Silber et al. 2018; Tabetah & Melosh 2018). Some of the smaller fragments produced by spallation or in fragmentation events will be vaporized. Most large meteoroids will break up in the atmosphere and some break ups are catastrophic.

Several theories describing fragmentation have been developed, including the splitting time concepts (Park & Brown 2012). It is arguable why time has any major role in splitting or fragmenting the meteoroid. Splitting time can be studied statistically but in a structural sense, meteoroids, even two meteoroids made of the same composition, may have very different histories and structural strength. Structural differences can also be seen in different parts within one meteoroid. We also question the idea of the incubation process described by Park & Brown (2012) where melted rock permeates through cracks of the meteoroid body, as there is lack of petrological evidence for that. If a crack opens under the atmospheric entry conditions it will cause the meteoroid to break apart and exposed fresh surfaces will rapidly be coated by a fusion crust. Dark shock melt veins seen in many stony meteorites are not produced during the atmospheric passage but much earlier, in space (see e.g. Kohout et al. 2014a). Melt veins form by shock-metamorphosis due to high pressure. To shock-melt a rock, peak pressures have to reach above 50GPa, which by far exceeds the maximum pressure meteoroids can withstand without breaking up (French 1998; Sharp & DeCarli 2006).

In case of the Košice meteoroid it has been suggested, and some evidence supports, that the meteoroid was composed of two pieces before entering the Earth's atmosphere (Borovička et al. 2013; Gritsevich et al. 2014c). This may well be the case as most meteoroids have experienced impact(s) during their existence. Hence, assuming one solid body for a large meteoroid may not be the default case.

Even a small, space-based impact can fracture a large meteoroid body so that it will break up during the subsequent atmospheric passage. This may happen already in considerably low air pressure flow regimes characterized by the dimensionless parameter called the Knudsen number (Moreno-Ibáñez et al. 2018). Photos of recent asteroid and comet flybys have shown that fractures, lithic breccia, and porous rubble piles are probably quite common among larger asteroids. Observed very low fragmentation pressures may be a sign of the porosity of a material (Meier et al. 2017), which decreases compression strength, or peeling of unconsolidated surface regolith as seen on the photos of the asteroid Bennu (Molaro et al. 2020).

When fragmentation occurs, velocity components are added into the fragment's flight trajectory, which cause it to deviate away from the original trajectory. These velocity components are so small that atmospheric drag will zero them very fast. That is why pieces from a fragmented fireball seem to continue along their original trajectory. Sometimes deviation downwards from the original trajectory can be clearly seen because air drag slows down and gravity assists smaller trailing fragments.

After any fragmentation, ablation will decrease the fragments in size and some will completely vapourize/disappear. This will reveal how much mass and how many fragments could be expected to reach the ground. However, uncertainty of strength of the meteoritic material in a fall means that error margins are notable. Running several separate DFMC simulations is recommended to figure out an average outcome.

When the fireball experiences multiple fragmentation events, it is preferable to initiate several simulations from different starting points. If the fireball has broken up into several individual fireballs, it may be worth running separate simulations for few of them. Although the result is usually that the strewn fields will overlap with each other, these simulations may provide additional information of the mass distribution within the strewn field. One solution would be that all break ups and tailing fragments are included in the simulations. This is doable and can provide more details about the mass distribution in the strewn field.

Gritsevich et al. (2014c) provided insight of distribution laws for the Košice meteorite fall. Based on this, the bimodal lognormal distribution was initially used as the fragmentation probability function in our MC model to give a more realistic mass distribution for simulations. Later, this was replaced with the more generic distribution described by equation (2). One reason was that the observation of the distribution of fragment sizes on the ground is not necessarily the same as the distribution of fragment sizes during the luminous trajectory. Another reason was that the model is meant to be generic – able to produce the forward prediction – and not restricted only to reproducing the scenario of the Košice meteorite fall. However, when the model gives large masses (tens of kilos) it tends to produce very few small fragments. This can be adjusted by additional runs. In the case of the Košice fall, it is good to recall that besides all the search efforts it is probable that a lot of fragments have not been found or properly documented (Gritsevich et al. 2014c; Tóth et al. 2015).

In our DFMC model, fragmentation is triggered at the start point or very soon after that. Every simulated fragment undergoes some level of fragmentation. This approach was chosen to account for possible surviving fragments from early fragmentation events.

To model the fragmentation, we use the exponential distribution function described by equation (2), where µ determines the ratio of cumulative mass of the produced largest fragments to the mass of the progenitor. Variable χ is a random test value between 0 and 1. The value of the exponent ί determines the shape of the probability curve.

Depending on the scenario, there may be a need to adjust µ and ί to get a more even distribution of the fragments, since our goal is to produce a map of the strewn field indicating the locations and sizes of the survived meteorite fragments. We usually use values between 0.02 and 0.05 for µ. If the value is too high, there are lots of large fragments and fewer small fragments. We usually adopt values from 4 to 7 for ί, which has a similar effect with low values.

The mass of the newly produced fragment is introduced as follows. A random number between 0 and 1 is generated. If the number is less or equal to the result of equation (2), this random number will be the normalized mass value we use for the fragment after the fragmentation. The fragment's new mass is simply this random number multiplied by the mass of the progenitor. After mass determination of the fragment, its size has to be recalculated for the air drag equation (equation 4).

Fragmentation changes the direction and velocity of the fragment. The separation velocity depends on the size of the fragment. For a larger fragment, the separation velocity is lower than for smaller fragments. The separation velocity of the fragment is moderate. It is usually a few tens or couple of hundreds of meters per second. We use separation velocities of 10–190 m s⁻¹ with a mean value of 100 m s⁻¹. A separation direction of the fragment and its velocity is added to the original velocity and direction of the fragment.

Small fragments with masses below 1 g do not play a useful role in meteorite recovery. Occasionally they can be found, but there is no sense to go after them. The mass of the smallest recovered fragment documented in the Košice meteorite fall is 0.3 g (Gritsevich et al. 2014c). In our code, fragments weighing less than 0.3 g will be deleted. One reason to dismiss very small fragments is that they may end up in a nearly suspension-like state leading to an excessively long settling time. This usually happens when a small fragment gets ‘attached’ to a high drag coefficient value and air drag will produce the force which is about the same as the gravitational force.

A prediction of how much material has landed on the ground can be made by using the starting mass as the total cumulative mass of fragments in the model. At the start point DFMC simulation will include a number of fragments satisfying this criterion. In cases where there is evidence for, or suspect for, some fragments trailing the main mass at the start point, an adjustment upwards can be made.

In our DFMC model, the breaking of a fragment at the free fall stage is not allowed, although it does happen in reality. Fragments that break off after the ablation stage in mid-air or when impacting on the ground can be easily recognized based on the presence of fresh fractures not covered with a fusion crust.

2.8 Terminal height filter

A strewn field produced by the DFMC simulation may be tens of kilometres long with lots of large fragments. A long strewn field is a common result in the case of a shallow entry angle, but in many cases that may not be the reality. A long strewn field may be due to the fact that simulations produce unrealistically big fragments. A filter to prevent these large fragments is the terminal height filter. It can be set to the lowest altitude where the fireball was observed. The observed terminal velocity (usually 2–4 km s⁻¹ if a meteorite is going to survive) at the end of the luminous trajectory can be used as a terminal height filter. In that case, fragments that have higher velocity at the terminal height than velocity retrieved from the observations, should be removed from consideration. In cases where terminal velocity is not available from observations, we use 3 km s⁻¹ as a default threshold for the terminal height filter.

An alternative algorithm to estimating the terminal height can be based on predicting its value based on the properties of the meteor trajectory (derived, for example, based on observations). Such methods were earlier developed by Gritsevich & Popelenskaya (2008) and Moreno-Ibáñez, Gritsevich & Trigo-Rodríguez (2015, 2017).

3 RESULTS

The described DFMC model has successfully assisted in the recovery of three recent meteorite fall cases – Annama, the asteroid 2018 LA, and Ozerki (Gritsevich et al. 2014a,b; Trigo-Rodríguez et al. 2015; Kohout et al. 2017; Maksimova et al. 2020; Jenniskens et al. 2021). Also, our simulations for the Flensburg meteorite produced a match before it was revealed that the meteorite had been found. Our code has worked well for several other meteorite falls when tested after the first meteorite recoveries have been made.

3.1 Košice meteorite fall

In this study, we demonstrate our results by analysing the example of the Košice meteorite fall. The Košice meteorite, a well-documented H5 ordinary chondrite, fell on the 2010 February 28 at 22:25 ut. The fireball was fortunately caught by several surveillance cameras when actual fireball cameras from the European Fireball Network were under overcast sky (Borovička et al. 2013; Kohout et al. 2014b; Gritsevich et al. 2014c; Tóth et al. 2015; Gritsevich et al. 2017). This meteorite fall produced a large number of fragments. Fragmentation started at an altitude of 56.8 km and the last noted fragmentation was at an altitude of 21.6 km. Borovička et al. (2013) concluded that some of the meteorite fragments may have experienced up to six fragmentation events along the trajectory.

3.2 Input data

The exact GPS locations and weights for 79 meteorites out of over 218 found of the Košice fall have been published by Tóth et al. (2015) and are shown in Fig. 3. The remaining fragments coordinates were kindly provided by Juraj Tóth (private communication) and are used here for visualization purpose (Fig. 3). A coordinate fixed map of 218 meteorites has also been published in Borovička et al. (2013) where values for the fireball trajectory were also presented. These values are used as the start point parameters for the MC simulations presented in this work (Table 3).

Figure 3

The Košice strewn field. Black solid circles represent the recovered meteorite fragments (Tóth et al. 2015; Tóth, personal communication). In this set of 218 meteorite fragments, the two largest ones have masses of 2.374 kg (A) and 2.167 kg (B). The total number of the recovered fragments exceeds 218 fragments (Gritsevich et al. 2014c). Open circles are the (91) fragments simulated in this work using the trajectory parameters for the beginning height of 68.3 km from Borovička et al. (2013). Red square (N) stands for the virtual 38.8 kg nominal body, defined in Section 2.3, and the two largest red circles correspond to 2.344 kg (C) and 2.183 kg (D) simulated fragments. Vertical axis shows degrees (north) in latitude and horizontal axis shows degrees (east) in longitude.

Table 3.

The DFMC input parameters for the Košice meteorite fall. The trajectory data are taken from Borovička et al. (2013). In the DFMC code, error margins are not applied individually to the coordinates of the start point. The error in the coordinates of the start location is accounted for as a cumulative spatial error at the start point.

Symbol	Parameter	Mean value	Error margins	Remarks
λ0	Longitude	20.705°E	(±0|${_{.}^{\circ}}$|011/∼600 m)	WGS84 (error margins used as e0)
φ0	Latitude	48.667°N	(±0|${_{.}^{\circ}}$|021/∼2300 m)	WGS84 (error margins not used in simulation)
h0	Altitude	68 300 m	(±1400 m)	(Error margins not used in simulation)
e0	Spatial error at the start point	0 m	±600 m	Error used here is the error from longitude
δ0	Direction of trajectory	71°	±4°	0°–360° (0° = N, clockwise)
γ0	Trajectory slope	59|${_{.}^{\circ}}$|8	±2°	0°–90° (90° = vertical)
V0	Velocity	15 000 m s−1	±300 m s−1	Velocity at the start point
a0	Deceleration limit	8 m  s−2	a ≥ 0.99 × 8 m s−2	Deceleration at the start point
heh	End height of fireball	17 000 m	–	Lowest observed altitude rounded downward
ρm	Density of meteoroid	3300 kg m−3	±500 kg  m−3	Our default for chondrites
σ	Ablation coefficient	0.005 s2 km−2	–	From Borovička et al. (2013)

The accuracy of the two given points on the luminous trajectory of the Košice fireball has some limitations. This is due to the fact that the cameras that captured this fireball were security cameras not dedicated to the observation of fireballs. In table 4 in Borovička et al. (2013), the height for the first point is 68.3 ±1.4 km and for the terminal point it is 17.4 ±0.6 km. The coordinates have error margins ±2.3 km for the latitude and ±0.8 km for the longitude for the first point and error margins are ±1.1 and ±0.36 km for the latitude and the longitude of the terminal point.

The first given point of the trajectory at an altitude of 68.3 km was chosen as the start point for the DFMC simulation. We used the error margins given by Borovička et al. (2013) as one σ deviation for the parameters (Table 3). We use the latitude and the longitude values as given, and set the error margin for the start point. We use the error margin of the longitude as a deviation for the start point location, which in this case is ±0.6 km. Alternatively, we could use here the error margins of the latitude or height as deviation for the location, but in case of Košice those are so large that they would cause an unnecessary spread in the resulting strewn field.

The atmospheric data generated (for the coordinates 48.5°N, 21.0°E on 2010 February 28 at 21Z) for the simulation unfortunately only cover heights from 30.6 km downward. Radiosonde data from Poprad-Gánovce observatory (station 11952, at 00Z 2010 March 1) in Slovakia reach up to 27.9 km only. The upper part of the simulation is therefore following the standard atmosphere model adjusted to continue seamlessly with the generated atmospheric data. In this case, we used radiosonde data downloaded from the University of Wisconsin Atmospheric sounding page as well, since the difference between these data and the generated atmospheric data is negligible. We provide these atmospheric data in the supplementary material.

Borovička et al. (2013) provide all the luminous segment of the trajectory parameters needed for a DFMC simulation (Table 3) except the deceleration values. This is a common problem since at high altitude at the beginning of the luminous trajectory, deceleration is low and it is often difficult to numerically estimate from observations (Halliday, Griffin & Blackwell 1996). Deceleration at the start point of the simulation can be calculated by using the known velocity and height of the meteor. Also, analytical methods to calculate deceleration exist (Gritsevich 2008a, 2009); Gritsevich & Stulov (2006) and can be used to assist in computations. Other methods, like the terminal height filter described in the previous section, can additionally be applied. For this study, we adjust the ‘observed’ deceleration value in the start point so that the largest masses we get are around 3 kg which is a little bit more than the largest known Košice meteorite fragments, because our model allows only a single fragmentation event per fragment. We chose to set deceleration at the start point as 8 km s⁻². We also apply an ablation coefficient of 0.005 s² km⁻² from Borovička et al. (2013) instead of our 0.0018 s² km⁻² default value.

3.3 Initial mass estimate

With these parameters, our starting mass value at the start point is derived as 1056 kg (see Section 2.3). This is more than three times less than 3500 kg initial mass estimated by Borovička et al. (2013). On the other hand, our mass estimate is closer to the preatmospheric Košice mass estimate of 1840 kg by Povinec et al. (2015) and 1850 kg by Gritsevich et al. (2017). These differences in mass estimates can already be explained by the differences in the initial assumptions. For example, the product of the pre-atmospheric shape factor and the drag coefficient is required in the preatmospheric mass estimate approach described in Gritsevich et al. (2017). This value was set to 1.8, similar to the calculations for the Annama and the Park Forest meteoroids (Trigo-Rodríguez et al. 2015; Lyytinen & Gritsevich 2016b; Meier et al. 2017). If however, we were to recalculate this for the spherical shape assumption in this paper, the 1850 kg estimate by Gritsevich et al. (2017) could be subsequently revised by a factor of 1.21/1.8, which yields 1244 kg.

3.4 Strewn field prediction

Fig. 3 shows the result of the DFMC simulation for the Košice fall using the input trajectory parameters shown in Table 3. A good match is observed with the actual strewn field reconstructed using the coordinates of the recovered Košice fragments (Tóth et al. 2015). This simulation produced 91 fragments on the ground with a cumulative mass of 19.498 kg. The largest simulated fragment (C) in Fig. 3 is 2.345 kg. The simulation yielded another three large fragments [2.18 kg (D), 1.98 kg, and 1.73 kg] . In reality, some of those fragments have probably fragmented further into smaller pieces as our code in its present configuration allows no more than one fragmentation event per fragment. Accounting for the secondary fragmentation events is beyond the scope of the present paper. Although it would likely yield a better match between the number of the collected and simulated fragments, it would not drastically improve the quality of the produced map.

The simulation in Fig. 3 shows the typical strewn field map that would usually be used to guide any search effort in the field. If/when meteorites are found, as happened in this case, prediction of the strewn field can be seen as subject to further adjustments/improvements. For example, by accounting for the measured bulk density of fragments and other peculiarities. However, in our case (and as can be seen in the results in Fig. 3) there is no need for adjustment, hence the given input parameters for the trajectory (Table 3) are good enough to produce a realistic strewn field fit. The simulated fragments provide excellent representation of the meteorite fall as they fully cover the area on the ground where the Košice fragments have been found.

Despite the good match to the known strewn field, its extent may seem to be larger than the extent of the documented meteorite finds (Fig. 3). Our overarching aim is to be able to realistically reproduce the strewn field for any meteorite fall (rather than cover one particular case with unprecedented detail) and therefore we categorize this outcome as the most appropriate. This result is natural due to the error margins imposed by the reconstruction of the observed fireball trajectory. Future studies may be focused, however, on fitting simulation results more closely to the known distribution of meteorites on the ground by restricting some of the trajectory parameters. For example, in this study we used our default bulk density (3300 ± 500 kg m⁻³) of chondritic material, although the bulk density of the Košice meteorite has been measured as 3430 kg m⁻³ with standard deviation of 110 kg m⁻³ (Kohout et al. 2014b) and it could be used instead when seeking a perfect fit.

Although finding the perfect fit is not the aim of this study, we noticed that if we apply our default error margins (±1|${_{.}^{\circ}}$|8) for the slope angle and direction of the trajectory with the mean slope angle as 59|${_{.}^{\circ}}$|5, the fit is nearly perfect.

To estimate the total mass of the fragments on the ground, we ran a total of 21 individual DFMC simulations with the same initial parameter set (Table 3). Taking the average of the number of fragments, the weight of the largest mass and terminal mass on the ground, we get an estimate of how much material could be recovered. The simulations produced on average 64.8 fragments to survive down to the ground. The mean weight of largest meteorites on the ground in each simulation is 3.65 kg and mean cumulative terminal mass on the ground is 14.71 kg (deviation 3.03–26.50 kg). That means that only 0.3–2.5 per cent of the starting mass survived down to the ground according to our model. The average mass of a fragment that reached the ground is 0.227 kg (out of all 1360 simulated fragments). This estimate likely represents the upper bound due to the one fragmentation event per fragment limitation present in our current version of the code. However, this is the best for a current terminal mass estimate as it provides a closer match to what is expected to be found on the ground compared to the earlier terminal mass estimates obtained at the terminal point of a fireball.

4 IMPLICATIONS FOR METEORITE RECOVERY: EXAMPLE OF NEUSCHWANSTEIN

In the preceding section, we modelled the strewn field of the Košice meteorite, which presents an excellent opportunity to demonstrate performance of the model due to the large number of recovered and well documented fragments. On the other hand, the Košice fireball was poorly observed, and consequently the uncertainties on the state vector of the meteoroid were high (Table 3). These uncertainties imposed on the input data subsequently also enlarge the calculated strewn field. In practice (and even more so under conditions of constantly growing observational means), the opposite situation is more likely: a well-observed meteorite fall from a well-derived luminous trajectory of the fireball with tight margins, but a poorly observed strewn field itself, due to the small number (or absence) of recovered meteorite fragments.

To apply our model to such a meteorite fall with more accurate trajectory parameters, we simulate the Neuschwanstein meteorite fall on the 2002 April 6 at 20:20 ut. Spurný, Oberst & Heinlein (2003) provide accurate measurements for the beginning and the terminal point of the fireball trajectory (Table 4). This bright fireball was recorded by several cameras and was unusually well documented by radiometers, infrasound detectors, and nearby seismic arrays. Three rather large meteorite fragments are known to be recovered: Neuschwanstein I (1.750 kg)  was found on 2002 July 14, Neuschwanstein II (1.625 kg)  on 2003 May 27, and the largest Neuschwanstein III (2.843 kg)  fragment was found on 2003 June 29 (Heinlein 2004; Oberst et al. 2004; Gritsevich & Stulov 2008). The fragments were found on both sides of the border between Germany and Austria and the corresponding reported WGS84 coordinates of the finding sites are: (I) 47.52392°N 10.80803°E, (II) 47.53386°N 10.80817°E, (III) 47.51618°N 10.82158°E (Fig. 4). Also as before, these coordinates are not integrated with our model in any way and are only shown here in the context of the ground truth comparison with the theoretical results obtained using the DFMC model.

Figure 4.

The Neuschwanstein strewn field shown on the Google Earth map. White bullseyes represent the three recovered meteorite fragments (Heinlein 2004; Oberst et al. 2004). These are marked as (A) Neuschwanstein I (1.750 kg), (B) Neuschwanstein II (1.625 kg), and (C) Neuschwanstein III (2.843 kg). Coordinates for Neuschwanstein I and II are given in Oberst et al. (2004) and the location of Neuschwanstein III is from Heinlein (2004). Open circles are the fragments simulated in this work using the trajectory parameters for the beginning height of 84.95 km from Spurný et al. (2003). Red square (N) stands for the virtual 56.6 kg nominal body, defined in Section 2.3, and the five orange and yellow circles correspond to 8.05 kg (D), 7.72 kg (E), 2.35 kg (F), 1.29 kg (G), and 1.00 kg (H) simulated fragments. Colour codes for the simulated fragments are orange 3–10 kg, yellow 1–3 kg, green 0.3–1 kg, cyan 0.1–0.3 kg, and blue <0.1 kg. Background map by Google/GeoBasis-DE/BKG 2020.

Table 4.

The DFMC input parameters for the Neuschwanstein meteorite fall. The trajectory data are after the table 1 in Spurný et al. (2003). These parameters are for the beginning of the luminous trajectory.

Symbol	Parameter	Mean value	Error margins	Remarks
λ0	Longitude	11.5524°E	(±0|${_{.}^{\circ}}$|0009/∼70 m)	WGS84 (error margins used as e0)
φ0	Latitude	47.3039°N	(±0|${_{.}^{\circ}}$|0006/∼70 m)	WGS84 (error margins not used in simulation)
h0	Altitude	84 950 m	(±40 m)	(Error margins not used in simulation)
e0	Apatial error at the start point	0 m	±70 m	Error used here is the error from longitude
δ0	Direction of trajectory	295|${_{.}^{\circ}}$|3	±1|${_{.}^{\circ}}$|8	0°–360° (0° = N, clockwise). Default error
γ0	Trajectory slope	49|${_{.}^{\circ}}$|75	±0|${_{.}^{\circ}}$|03	0°–90° (90° = vertical)
V0	Velocity	20 950 m s−1	±40 m s−1	Velocity at the start point
a0	Deceleration limit	3.5 m  s−2	a ≥0.99 × 3.5 m s−2	Deceleration at the start point
heh	End height of fireball	16 040 m	(±30 m)	Lowest observed altitude (error not used)
ρm	Density of meteoroid	3300 kg m−3	±500 kg  m−3	Our default for ordinary chondrites
σ	Ablation coefficient	0.0018 s2 km−2	–	Our default value

Simulation parameters for the start point of the Neuschwanstein meteorite fall are listed in Table 4. As we present it here as a benchmark case, we have once again used our default bulk density value for ordinary chondrites, although the Neuschwanstein meteorite is an enstatite chondrite (EL6) with a bulk density of 3490 kg cm⁻³ and this, in principle, could be specified. Values of ablation coefficient and deceleration were not provided by Spurný et al. (2003). Hence, as an ablation coefficient we used our default value of 0.0018 s² km⁻², which in the course of the study, proved to be a reasonable estimate for enstatite chondrites. Indeed, enstatite chondrites are usually harder rocks when compared to most ordinary chondrites.

Deceleration at the start point was chosen so that simulations produce largest meteorite fragments weighting 5–7 kg on average. This range is based on the total known weight (6.218 kg)  of the recovered Neuschwanstein meteorites. Because the start point of this simulation was at the altitude of 84.95 km, deceleration has to be relatively low and a good match was obtained with a deceleration value of 3.5 m s⁻². This value corresponds to the mass at the start point of 605.4 kg, which is consistent with the initial mass range estimate of 540–672 kg given in table 7 by Gritsevich (2008a). Lower deceleration values also give a good match, but would mean a larger nominal mass at the start point. However, there is no reason to make simulations look overly optimistic in terms of how much mass could potentially be recovered. Alternatively, and when we are predicting the strewn field before any meteorites are actually found, deceleration is usually the integral part of a trajectory retrieval that could be calculated analytically (Gritsevich 2008a, 2009). In the case studied here for Neuschwanstein, the goal was to see how well our model matches the meteorite masses which were already found on the ground.

Atmospheric data for these simulations are provided in the supplementary material. The data of München-Oberschlssheim weather observatory (station 10868, at 00Z 2002 April 7) were obtained from the University of Wyoming online weather data base (http://weather.uwyo.edu/upperair/sounding.html). Radiosonde data are available up to altitude of 33189 m. Atmospheric conditions upwards from that altitude are substituted with the standard atmospheric model. Some readings are missing in the data below the mentioned altitude, but those are not relevant for the simulations.

The direction of the trajectory and its error margins were not provided by Spurný et al. (2003). Therefore, we have calculated it as the direction (or azimuth) of a line connecting the beginning and terminal points of the luminous trajectory. We first used error of ±0|${_{.}^{\circ}}$|5 for the direction, which can be calculated using the error margins given for the coordinates of the beginning and terminal points. However, the resulting strewn field appeared rather narrow. All three recovered meteorites were located on the edge of the simulated strewn field. That could be the case also in reality. However, in our demonstration in Fig. 4 we chose to use our default error margin of ±1|${_{.}^{\circ}}$|8 for the direction of the trajectory. This increase of the error margin stretches the predicted strewn field wider, as shown in Fig. 4. In addition, we find the use of our default error margin of ±1|${_{.}^{\circ}}$|8 for the direction of the trajectory sensible, since the DFMC simulations are usually done before any meteorite has been found in the field.

Additional alteration of the width of the strewn field comes from the variations in the wind direction and wind speed, some changes may result also from the changes to the trajectory direction that occur during the fragmentation events. As the case of the Neuschwanstein demonstrates, error margins in well-documented fireball trajectories can be very small [error given by Spurný et al. (2003) for the trajectory slope is only ±0|${_{.}^{\circ}}$|03] and this is why, at first glance, our default ±1|${_{.}^{\circ}}$|8 deviation in directions may not always seem well founded. However, its primary function is to make the strewn field prediction sensible – covering the full extent of the probable area where meteorite fragments may have landed.

Already the first simulation produced a good match for the Neuschwanstein case. We also got results where the larger meteorite fragment landed close to the smaller ones, as was the case of largest found Neuschwanstein III meteorite. In all cases, the DFMC simulations illustrate where additional meteorite fragments could have potentially been recovered.

We chose one simulation for the illustration in Fig. 4. It produced a couple of fragments (D and E) that are larger than any of the found Neuschwanstein meteorites. This shows an area where larger meteorites may have landed after the fall. This also demonstrates that if a second fragmentation event is allowed in the model, it could produce Neuschwanstein I and II sized meteorites in the right area. In Fig. 4, the simulated strewn field is 25 km long and 4.5 km wide in its middle part (when the 56.6 kg nominal mass is ignored). There are 95 simulated fragments on the ground giving a total mass of 24.89 kg in this simulation. This value compares well with the earlier mass range estimate of 22–53 kg, obtained for the end of the luminous flight, given in table 8 by Gritsevich (2008a).

In the series of 20 DFMC simulations, the nominal mass at the start point is 605.4 kg. The mean cumulative mass on the ground in these simulations is 19.95 kg (deviation 11.98–27.47 kg)  which means that only 2.0–4.5 per cent of the starting mass survives down to the ground. The number of fragments per simulation is, on average, 121.4 and the mean mass of the largest fragment is 5.33 kg, while the mean weight of all 2427 simulated fragments is 0.164 kg.

5 DISCUSSION

All MC methods are based on a repeated process. MC methods rely on a large number of samples. The more repetitions that are made, the closer the result is to the realistic distribution. This is because relative errors of an MC integration decrease when the number of samples increases. Usually the result of a MC simulation is presented as a distribution map of probability of the phenomena. In order to obtain a representative gradient map, one has to simulate thousands of trajectories or more.

Our DFMC simulation was first implemented in such a way. That is how the classical MC simulations are made, but MC simulations can also be used differently. MC methods can be restricted to just one part of the phenomena and let the other parameters dictate how many repetitions will be made. When simulating a meteorite strewn field we do just that. We use MC variations for individual fragments which we trace down to the ground. That means that each fragment will have its own unique trajectory. Instead of thousands of trajectories, we limit the number of trajectories to represent the number of fragments that probably fell and can be found on the ground. Size and other properties of the fragments are affected by the MC algorithm during the fragmentations and the total mass of the fragments is limited by the starting mass as calculated above.

Our trajectory tracing DFMC simulation was developed to help the Finnish Fireball Network's meteorite recovery efforts. It was an addition to the numerical solution which was used to plot the central line of the strewn field. However, a large part of Finland is covered by water bodies and in some parts of the country undergrowth may be so dense that one can walk over a meteorite without seeing it. That is why searches should be conducted in search-friendly areas where meteorites can be spotted on the ground in the first place.

Andreić (2011) estimated that errors in velocity and deceleration in the beginning and during the dark flight cause the largest errors in simulations only when the wind is not accounted for. Our simulations agree with this. The main contributor for the strewn field shape, size and location is atmospheric winds. This is evident in the way narrow and long strewn fields bend due to winds when a traditional concept of a distribution ellipse is dismissed. If fragment shape, orientation or break up played a bigger role, then the wind drift effect would not be as clear. Although odd fragment shapes, wild orientation or fast rotation may explain some meteorites found outside the main distribution area of meteorites from the same fall. Of course, these forces have to be included if trajectory tracing backwards is to be done for a known meteorite fragment but then the true shape and Reynolds number and other factors are known.

The most common shape of a strewn field is a narrow long tie shaped area as shown in Fig. 4. This area is usually bent due to atmospheric winds. Knowing the shape of a strewn field may help to locate more meteorites also from past meteorite falls. The widely used term of a distribution ellipse is outdated and may be only valid for special cases.

Our DFMC model still has room for improvements, discussed below, but the necessity of some of them is questionable. Fully accurate physical modelling may, from a scientific standpoint, seem to be what is hoped for. However, in many practical applications the key ingredient for success is appropriate division of involved physical processes into primary and secondary. To predict a meteorite strewn field accurately enough for the search effort, some details may be ignored. Magnus effect during the free fall is one of them. Although some of the mapped strewn fields may exhibit outlying fragments which are difficult to explain in any other way, they are not very common.

There are reasons to apply more realistic shape for the nominal body. One of them would be the possibility to match the starting mass of the nominal body better with observed deceleration of the meteoroid. Because of ablation, any fixed shape of the body would not be more plausible than assuming a sphere. Implementing a dynamic shape modification right after fragmentation could also be an improvement. When a fragment breaks from the larger body, its shape may be aerodynamically very poor, causing a jump in the drag coefficient. However, these fragments may break up again in a subsequent fragmentation event and/or ablation may rapidly reshape them. The shape change coefficient is sometimes introduced to account for the change in meteoroid shape along the flight (e.g. Levin 1956; Stulov, Mirskii & Vislyi 1995; Gritsevich 2009; Sansom et al. 2019b; Drolshagen et al. 2020). Hence, a test can be thought to apply the model enhanced by the variation in shape change coefficient.

One improvement to our code would be implementing a true fragmentation history of a fireball. It could be done by interpreting the light curve of the fireball. It is known that each flash or flare of a meteor is a sign of the fragmentation event. A theoretical light curve of the non-fragmenting nominal body could be adjusted to closely fit the observed light curve of the fireball and then subtracted from it. This should produce a probability curve for fragmentation. This kind of probability curve could also be thought when multiple fragmentations are imposed.

One remaining open question is the better understanding of wind during the dark flight. How much does wind speed and direction vary while a cloud of meteorite fragments falls through certain altitudes compared to the obtained/available atmospheric data?

6 CONCLUSIONS

The MC method based modelling can be successfully applied to predict strewn fields produced by meteorite falls. Such predictions work well in cases where good quality start point data (coming from the analysis of the fireball) and corresponding atmospheric data are available. Creating a probability distribution of fragments in the strewn field improves the possibilities of finding recently fallen meteorites. Using different kinds of filters gives an estimate of how much meteorite material is expected to have landed on the ground and if the labour effort required for a field search can be effective. Some DFMC simulations fail to produce any surviving meteorites or surviving meteorites are very small, even from a bright fireball. That is not a sign of a failure but rather it implies that the fireball did not produce any recoverable fragments. No prediction will make actual finding of a meteorite easier but may prevent spending time and resources in a wrong area.

This work, as well as the analysis of the maps of known meteorite strewn fields, reveals the importance of wind drift and its possible effects. Iron meteorites produce most of the narrow strewn fields less than 2 km in width. Strewn fields produced by stony meteorites are wider, yielding on average 3.8 km in width. The length of the strewn field is mainly affected by the entry angle of the meteoroid and wind direction. Clusters of meteorite fragments that exist in many strewn fields are likely a result of fragmentation events occurring along the trajectory. Shape, orientation, and character of motion of fragments do not play a big role on the shape and location of the strewn field.

In many historical meteorite falls, the assumption of a distribution ellipse of fragments has been erroneous and probably the true extent of the strewn field has not been adequately searched. Hence, the DFMC simulations may also help to find more meteorites from old falls where meteorite finds have been rare or lacking. In cases when several meteorites have already been found, knowledge of the predicted strewn field shape may significantly help to search the rest of the strewn field.

SUPPORTING INFORMATION

Atmospheric data for the Košice meteorite fall.

Atmospheric data for the Neuschwanstein meteorite fall.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

This work was supported, in part, by the Academy of Finland project no. 325806 (PlanetS). We acknowledge Hadrien Devillepoix and Morgan Hollis for their supportive and valuable comments that encouraged us to improve this paper by adding the section on the strewn field of the Neuschwanstein meteorite fall. We thank Eleanor Sansom for helping us to proof read this paper and for stimulating discussion. We thank Juraj Tóth for providing additional coordinates of the documented Košice meteorite fragments in addition to those published in Tóth et al. (2015). We thank Dieter Heinlein and Jürgen Oberst for discussing the circumstances of the known Neuschwanstein meteorite finds and difficulties of the ground search. We thank Panu Lahtinen (Finnish Meteorological Institute) and Peter Völger (Swedish Institute of Space Physics) for their help with obtaining the actual atmospheric data for the studied fireball events and valuable discussions. We thank all members of the Finnish Fireball Network and acknowledge Ursa Astronomical Association for the support with the Network coordination. We acknowledge fruitful collaboration with the colleagues at the Ural Federal University in organizing the field trips and prompt meteorite recoveries based on provided strewn field maps. The research at the Ural Federal University was supported by the Russian Foundation for Basic Research, project nos. 18-08-00074 and 19-05-00028.

This paper is dedicated to the memory of our friend and colleague Esko Lyytinen – a man of extraordinary mind with a life-long passion for science.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author. The atmospheric data obtained and used in this study, corresponding to the location and time of the Košice and Neuschwanstein meteorite falls, are appended to this article as supplementary material.

Author notes

Deceased.