Secondary cratering on Earth: The Wyoming impact crater field

The number of impact craters discovered on Earth so far is very small compared to that of the Moon or Mars and currently stands at 208 structures (Gottwald et al., 2020; Kenkmann, 2021). Causes for this paucity on Earth are erosion, burial, tectonic deformation, and that two-thirds of the Earth`s surface is covered with water with young oceanic crust underneath. Here, we confirm 31 small impact structures in southeastern Wyoming by the documentation of shock effects. More than 60 other possible structures have been detected that await confirmation. The partial preservation of the surface morphology of these craters allows insight into the mechanism of their origin. The craters observed in the top of the Permo-Pennsylvanian Casper Formation sandstones were buried soon after formation by soft mud and became quartz-cemented after burial. Uplift during the Laramide Orogeny in the Neogene and intense Quaternary erosion exhumed the craters.

We investigate how this local concentration of impact craters occurred. We will show that the clustering is due to the process of secondary cratering. Secondary craters are formed when highspeed ejecta blocks are excavated from a primary impact crater and collide with the surface at high velocity. Secondary craters were first recognized on the Moon (Shoemaker, 1965) and later on Mars, Mercury (Schultz and Singer, 1980), the icy satellites of Jupiter such as Europa (Schenk et al., 2004), Ganymede, and Callisto (Bierhaus et al., 2005), and Saturn's moons (Schenk et al., 2020). Until now they have not been found on Earth and other planetary bodies with dense atmospheres such as Venus or Titan. Atmospheres decelerate ejecta and could disrupt high-velocity ejecta. Here, for the first time, evidence is provided that secondary cratering has been possible on Earth. We reconstruct the physical processes that led to secondary crater formation. The field of crater structures discovered offers planetologists an unrivaled opportunity to study the processes of secondary crater formation on Earth.

The identification of secondary craters is straightforward in proximity to the source crater. Fragments impacting near the primary crater have relatively low velocities and produce irregular and non-circular shaped craters that are shallower than those of fresh primaries (Pike and Wilhelms, 1978). A best-fit for the depth-to-diameter ratio (d/D) of secondary craters on the Moon was found to be 0.11 (Pike and Wilhelms, 1978) and for Mars the mean value is 0.08 (McEwen et al., 2005). Characteristics of secondary craters are elongated crater clusters and V-shaped, herringbone patterns of ejecta accumulation produced by the interactions of simultaneous ejecta (Oberbeck and Morrison, 1973). Vickery (1986) constituted a size-frequency distribution of ejecta from large craters and showed that larger secondary craters were produced by clusters of smaller ejecta boulders. Fragments ejected with a relatively high velocity of up to a few kilometers per second land at large distances from the primary crater and form more regular and deeper craters. This makes it more difficult to distinguish them from small primary craters. However, distal secondaries still tend to be shallower and asymmetrical in cross section with the steepest crater walls in the uprange direction due to oblique impact (Gault and Wedekind, 1978). Furthermore, the alignment of secondaries and association with bright crater rays connect them with the source crater.

A prominent example of a rayed crater with abundant secondary craters is the large, 86-km-diameter, fresh Tycho crater situated in the nearside lunar highlands. The extensive system of bright rays contains dense clusters of secondary craters. Shoemaker (1965) estimated that Tycho produced 10⁵ secondary craters larger than 200 m in diameter, and Dundas and McEwen (2007) calculated at least 10⁶ secondary craters larger than 63 m diameter in its ray system. Hirata and Nakamura (2006) estimated the diameter and ejection velocity of the secondary-forming fragments from the crater size and the distance from the primary at Tycho. Zunil, a 10 km rayed crater in Elysium Planitia on Mars, has only a few secondary craters within 16 crater radii distance but a huge number of secondaries extend at least 1700 km in some directions (McEwen et al., 2005; Preblich et al., 2007).

Secondary cratering has gained considerable attention as it could affect the relative and absolute dating of planetary surfaces via crater counting. Bierhaus et al. (2005) suggested that 95% of the small craters on Europa with diameters of <1 km are secondary craters. McEwen and Bierhaus (2006) proposed that distant secondaries also govern the number of small craters on the Moon and Mars. The importance of secondary cratering was further supported by the recognition of self-secondary cratering on impact deposits within a given crater (Zanetti et al., 2017). In contrast, Neukum and Ivanov (1994) argued that the steep slope of the size-frequency distribution (SFD) of lunar craters in the size range of 0.3–1 km reflects the primary SFD of asteroids and is the result of fragmentation due to collisions in the asteroid belt (Hartmann, 2005). Ivanov (2006) and Bland and Artemieva (2006) showed that the currently observed meteoroid flux in the Earth–Moon system correlates with the SFD of small craters. Werner et al. (2009) concluded that crater count measurements are “contaminated” by secondary cratering by percentages smaller than the assumed statistical error.

Proximal secondaries on Earth are known from the Sedan nuclear explosion experiment (Roberts, 1964). The velocities of such secondary impacts are sub-sonic. Based on reflection seismic studies, Poag et al. (2004) interpreted fault-bounded depressions of 0.4–4.7 km diameter in the surroundings of the 85-km-diameter marine Chesapeake Bay impact structure, USA, as secondary craters. However, there is neither proof that these depressions are impact craters nor that they relate to the Chesapeake Bay impact event. Alleged small secondary craters have briefly been mentioned to be associated with the 8-km-diameter Bigach crater (Masaitis, 1999) and the 13-km-diameter Zhamanshin craters in Kazakhstan (Florenskiy and Dabizha, 1980), but no further information was given and no morphological expressions are present (Gottwald et al., 2020). Blocks of distal ejecta were reported from the Ries crater, Germany (Reuter, 1925). These so-called Reuter’ sche Blöcke are several decimeters in size and occur along with the finer-grained Ries Brock-Horizonte at distances of 70–200 km from the Ries crater (Buchner et al., 2007). The Acraman impact structure is also associated with a distal ejecta layer that contains 30-cm-sized blocks at 250–300 km distance from the source crater (Gostin et al., 1986). However, neither the Ries nor the Acraman impact structure have secondary craters.

Recently, Kenkmann et al. (2018) discovered a field of small impact structures (centered at 42°40′38″N, 105°28′00″W) in tilted Permian strata on the northeast flank of the Sheep Mountain anticline, which is part of the Rocky Mountains Front Range system near Douglas, Wyoming (Fig. 1). They reported a minimum extent of this crater field of 7.5 km × 1.5 km and identified more than 40 roughly circular to ellipsoidal possible impact structures therein. The confirmation of the impact origin of eight of the crater depressions was based on the documentation of planar deformation features (PDFs) and planar fractures (PFs) in quartz grains. Kenkmann et al. (2018) interpreted this cluster of craters as an impact crater strewn field that was formed by the break-up of a single meteoroid. Upon entry into the upper atmosphere, meteoroids moving at hypervelocity experience aerodynamic stresses. Once the aerodynamic stresses exceed the strength of a meteoroid, it will disrupt (Bland and Artemieva, 2006). The decelerated fragments then form a crater strewn field. The maximum theoretical spreading of strewn fields is restricted and should not exceed much more than 1 km perpendicular to the trajectory (Artemieva and Shuvalov, 2001). Of the seven crater strewn fields known so far on Earth, all are formed by iron-meteoroids. At Douglas, no iron projectile relics have been detected. Kenkmann et al. (2018) also hypothesized alternative scenarios in case further craters would be detected and the size of the crater field would become incompatible with that of an atmospheric breakup of a single meteoroid. In an abstract, Schulze et al. (2019) described additional craters that expanded the strewn field further.

Since the discovery of the craters near Douglas, additional craters have been identified on the Sheep Mountain Ridge and along other outcrops of the same strata in southeastern Wyoming. The expanded area measures ~90 km in a north-south direction and 40 km in an east-west direction. The city triangle of Casper, Douglas, and Laramie borders the crater field that is situated in Wyoming's Converse and Albany Counties (Fig. 1). We call it the Wyoming impact crater field. It includes the spatially separated occurrences at Sheep Mountain (SM), Wagonhound Ridge (WR), Mule Creek (MC), Fetterman Ridge (FR), Fetterman Road (FRX), and Palmer Canyon Road (PCR). In addition to these locations, more crater depressions have been identified (Fig. 1); however, a proof of those being impact-related is pending. The apparent diameters of the confirmed and potential crater structures range in size from 10 m to almost 80 m, with 32.9 m being the mean diameter (Fig. 2). Many craters are circular, but the ellipticity can reach up to 1.7. The mean ellipticity is 1.2 (Fig. 2). A weak negative correlation between crater diameter and ellipticity may exist but is statistically not supported.

All proven craters and crater candidates are exposed at the top of the Permo-Pennsylvanian Casper Formation near the contact with the Opeche Shale member of the Permian Goose Egg Formation and so are assumed to be of equivalent, 280 m.y., Leonardian (Kungurian) age. The target lithology is the uppermost Casper sandstone, which is a cross-bedded quartz arenitic sandstone that was deposited in a mixed eolian, fluvial to shoreface environment (Steidtmann, 1974; Cardinal and Holmes, 1984). At the time of impact, these rocks were unconsolidated to semi-consolidated wet sands. Shallow water conditions were likely present at some crater sites. The uppermost Casper Formation and the craters that were formed in these beds were buried beneath the Opeche Shale Member red beds of the Goose Egg Formation (Kenkmann et al., 2018). These mudstones were deposited in a quiescent transgression over the Casper sands and indicate an arid paralic, sabkha to lagoon depositional environment. The Casper Formation and overlying Opeche Shale Member are tilted to various degrees and are exposed on the flanks of monoclinal structures and vergent anticlines of the Laramie Mountain Range (Fig. 1). The Casper Formation and Opeche Shale Member experienced a long phase of burial. The reexposure of the original crater morphologies happened during exhumation of the rocks. Partial reproduction of the original crater shapes was possible by selective erosion because the mudstone cover was easy to erode. The craters sometimes form pedestal morphologies and appear more resistant to erosion than the normal Casper sandstone. In addition to diagenetic cementation of the Casper sandstone, this is attributed to (1) intense sealing of impact-induced fractures with quartzitic cement, (2) shock lithification (Short, 1966) of the central portion of the crater subsurface, and (3) massive chert formation in the crater surrounding that may result from possible melt agglutination at the surface (Kenkmann et al., 2018). Sealing of fractures and comminuted grain aggregates reached beyond the crater limits and was possibly assisted by elevated temperatures in the target shortly after impact. Shock lithification increased the cohesion of the rocks directly upon impact, but it was restricted to the immediate central crater subsurface (Wünnemann et al., 2016; Raith et al., 2018). Evidence for shock-lithification is given by crushed and comminuted grains that fit tightly against close-pressing neighbors, and partial interpenetration of adjacent grains (Kenkmann et al., 2018), similar to shock-lithification microstructures at the Wabar impact crater field (Gnos et al., 2013).

The Casper Formation exposures are limited in extent, comprising only 1% of surface exposures in Converse County (Douglas) and only 4% in Albany County (Laramie). Hence, substantial parts of the crater field are either still buried below the Goose Egg Formation and younger strata, or they have been eroded in the uplifted range.

The Laramie Mountains as part of the Rocky Mountains Laramide belt are characterized by thick-skinned, basement-cored uplifts surrounded by broad basins. The Laramie Mountains in the areas of this study are cored with exposed Archean metamorphic rocks. The tectonic deformation of the region took place during the Laramide Orogeny from Late Cretaceous to early Eocene (75 m.y. to 50 m.y.) (DeCelles, 2004). Estimates of the amount of shortening of the entire ~500-km-wide Laramide belt range from 40 km (Bird, 1998) to as high as 120 km (Chapin and Cather, 1983); however, across the study area major shortening in the WSW-ENE direction may not have exceeded a few kilometers and is the result of reverse faulting and folding. Minor counterclockwise rotation of crustal blocks around vertical axes, in combination with transpression, was suggested by Weil et al. (2016) for the region.

The Sheep Mountain crater location (SM), labeled the “Douglas crater strewn field” by Kenkmann et al. (2018), hosts the majority of craters, which are exposed along the northeastern flank of the NW-SE–striking Sheep Mountain anticline. Strata of the Casper Formation dip at roughly 15° to NE, parallel to the slope of the range. The craters occur stratiformly in the uppermost layers of the Casper Formation close to the contact with the overlying Goose Egg Formation. The Opeche Shale has eroded back, and further craters to the NE may be exposed in the future. The Sheep Mountain anticline shows a southwest vergence; the southwestern limb is partly suppressed by a reverse fault. Precambrian rocks are exposed along the center of the anticline. Based on the surrounding geology, it is estimated that 2000 m of overburden were eroded off the Casper during and after Laramide uplift (Fan and Carrapa, 2014).

The Sheep Mountain anticline is axially depressed southward but continues toward the northwest (Fig. 1). Here, the ridge is off-set by several faults and bends toward the east-west direction and represents the northern termination of the Laramie Mountain range. Paleozoic strata dip toward the north so that stratigraphic exposures of Pennsylvanian and Permian strata decrease in age with increasing distance from the Laramie Range. Near Box Elder and La Prele Canyons, southeast of Glenrock, further possible crater depressions are exposed along the Casper Formation/Opeche Shale Member contact (location BE in Fig. 1). Further west, at the Casper Mountain south of the city of Casper, no crater depressions have been found.

Wagonhound Ridge (WR), south of the Sheep Mountain Ridge (Fig. 1), forms an en-echelon array with Sheep Mountain and the Laramie Mountain Range further southwest. This suggests minor dextral transpressional slip in a NNE-SSW direction. The structure of the ridge is similar to Sheep Mountain, though less uplifted. The gently dipping northeastern limb provides large outcrops of the contact of Casper Formation and Opeche Shale Member and exposes several craters.

The other crater field locations occur along the southwestern slope of the basement uplift of the Laramie Mountains in transition to the Shirley Basin and its southern extension. From NW to SE, these are Mule Creek (MC), Fetterman Ridge (FR) Fetterman Road (FRX), and Palmer Canyon Road (PCR) (Fig. 1). In this region, the Casper Formation generally dips very gently to the southwest or to the south. The uppermost Casper Formation that hosts the craters is partly eroded, and degraded craters have been found on local buttes.

The crater field at the northeastern flank of the Sheep Mountain anticline contains a large number of circular, irregular-shaped, and ellipsoidal impact structures (Kenkmann et al., 2018). Many of them are now confirmed as impact craters (see Supplemental Material¹). All crater structures are exposed in the uppermost sandstone of the Permo-Carboniferous Casper Formation. They were exhumed after a long period of burial. The degree of preservation varies considerably from almost pristine to strongly degraded. The freshest structures contain steep crater walls, raised rims with overturned ejecta flaps, and remains of the proximal ejecta blankets (e.g., SM-1, SM-2, SM-4, SM-6, and SM-49) (Figs. 3A–3B). Their crater floors are covered by soil and filled with muds of the overlying Opeche shales. This makes robust measurement of depth (d) to diameter (D) ratios problematic. Estimates of the apparent d/D ratios are in the order of 0.1 and less and are determined in their present state. The freshest craters have distinct elliptical (SM-1: 60 × 41 m, R = 1.46, 327 ± 5°) to ovoid (SM-2: 29 × 24 m, R = 1.21, 322 ± 5°) (Fig. 3B) crater shapes that allow the reconstruction of impact directions. Their long axes have consistent orientations of 315–328° (see Supplemental Material). The craters have asymmetrical ejecta with well-developed overturned flaps on their NW side. This and the ellipticities define an impact direction from SE to NW. We observe irregular crater clusters (e.g., SM-8, SM-9, and SM-49) and crater chains, e.g., craters SM-3-4-5-6 (Fig. 3B), where craters partly overlap. The orientation of the crater chain SM-3-4-5-6 coincides in orientation with the long axes of elliptical craters. Occasionally straight walls of quartzitic breccia are observed near the craters, e.g., at SM-1 and SM-2, which are interpreted as possible relics of ejecta interference patterns that are also known as herringbone patterns (Figs. 3A–3B). More deeply eroded craters commonly occur upslope on Sheep Mountain at greater distance from the Opeche Shale contact, where the recessive erosion of the Opeche shale cover exposed the craters for a longer period. These crater structures are circular in outline (e.g., SM-9 (Fig. 3E: 20 × 20 m, R = 1), SM-15 (26 × 26 m, R = 1), SM-34 (35 × 35 m, R = 1), SM-36 (21 × 21 m, R = 1) (Supplemental Material)), and they show internal ring features instead of a central morphological depression. A halo of lighter quartzitic sandstone is characteristic of these craters. Generally, sandstone affected by impact appears somewhat brighter than the surrounding rock and is mostly more resistant to erosion owing to shock welding (Short, 1966; Kenkmann et al., 2018) and intense quartzitic sealing of these rocks (see Supplemental Material). This causes the pedestal appearance of some of the craters (e.g., SM-1 and SM-9) (Figs. 3A–3E).

As many as 10 possible impact structures were found in the uppermost Casper Formation, which is inclined by a few degrees to the NE, on the NE limb of the Wagonhound Ridge. Two confirmed craters show slightly elevated rims with outcrops that surround very gentle depressions filled with soil (Fig. 4A). The elevation difference between center and rim after tilt correction is less than 1 m. Crater WR-4 (58 × 47 m, R = 1.23, 329 ± 5°) and crater WR-5 (47 × 38 m; R = 1.23, 329 ± 5°) (Fig. 4A) are both moderately elliptical, and their long axes are oriented identically.

The most prominent, confirmed impact crater structure of the Mule Creek area is MC-1 that is formed in a tabular Casper sandstone target (Fig. 3C). The strata are tilted very gently toward SW (1–2°) due to the tectonic uplift of the Laramie Mountains. The crater appears as a slightly elliptical feature (56 × 44 m, R = 1.27) in a plateau, where the fractured bedrock is replaced by soil on the crater cavity. The long axis trends 284 ± 5°. The topography of this and the adjacent MC-2 (30 × 27 m, R = 1.1) structures is completely inconspicuous, and elevated rims and breccia are not preserved. It is assumed that a deeper subsurface level of the crater is exposed. The surrounding target shows distinct concentric fractures like a halo around the crater up to a radial distance of 100 m (Fig. 3C). This suggests that these rocks underneath the craters were consolidated to some degree when the impact happened. Apart from the concentric fractures that are interpreted to be impact related, a north-south and northwest-southeast–trending tectonic joint system is visible in remote sensing images. Such Laramide joints serve to initiate erosion during exhumation. The same bedrocks continue ~300 m northwest of MC-1. This tectonically jointed surface contains at least nine irregularly shaped depressions filled with soil that may belong to the crater field. A crater cluster with one confirmed crater (MC-26) has been found ~2.5 km southeast of MC-1 on Casper Formation. The state of preservation of these craters, however, is poor. The crater rims are composed of sandstone breccia sealed with chert matrix and consist of very resistant quartzitic sandstones. The microcrystalline chert shows polished wind-blown surfaces and scouring. The straight wall of brecciated quartzitic sandstone between the possible craters may represent relics of herringbone patterns.

The Fetterman Ridge crater structures are situated on erosional buttes that preserve the upper-most layers of the Casper Formation. One of the buttes, with a ~200-m-long and ~100-m-wide hilltop and elevation difference of ~25 m to the surroundings, hosts the confirmed impact craters FR-1 (30 × 22 m, R = 1.36) (Fig. 3D), FR-2 (28 × 17 m, R = 1.65), and FR-3 (10 × 10 m, R = 1). The target rock is layered and cross stratified aeolian Casper sandstone that shows intensive tectonic jointing. FR-1 and FR-2 appear elliptical in outline with their long axes trending WNW-ESE. However, the southern rims are eroded. FR-3 is irregular in outline. FR-2 shows a distinctly elevated crater rim of a few meters height in the NE and W part. The W rim shows an ejecta flap. The crater rim consists of quartzitic sandstone that is intensely brecciated at the decimeter-scale and locally contains injection veins with lithic breccia. In these veins, breccia clasts float in a matrix of crushed quartz grains. The brecciated and fragmented rocks are sealed by microcrystalline silica, which makes these rocks resistant to erosion. This, along with possible shock lithification (Kenkmann et al., 2018) in the central portion of the crater subsurface, caused the crater structures to form pedestal morphologies. The surfaces of massive chert occurrences are sculptured to ventifacts by the dominating strong southwestern winds of the region (Wyoming Water Resources Data System, 2020, http://www.wrds.uwyo.edu/). Tectonic jointing of bedrock on the flat-topped buttes led to the disintegration of strata during exhumation into meter-sized rhombic or rectangular blocks. More erosional buttes occur ~1 km further southeast and 3 km northwest of the FR 1–3 locations. This southeastern location is a ~600-m-long and ~250-m-wide plateau. It shows additional possible crater structures. However, the degree of erosion is high. This is also the case for additional buttes in the northwest.

West of Fetterman Road occur six small possible craters of 7–30 m diameter, some with distinct ellipticities (Fig. 3F). The craters are situated on the gently dipping slopes formed by the uppermost Casper Formation that dips here at 4° in the SSE direction (165°). Crater FRX-20 (25 × 15 m, R = 1.67; 296°) (Fig. 3F) has a pronounced ellipticity with the long axis trending ESE. Its rim is elevated by roughly 1 m with respect to the crateŕs interior and gently uplifted with respect to the surroundings. Adjacent to this crater, very small, more irregularly shaped crater depressions occur. Some 300 m east, a strongly elongated crater depression seems to have resulted from two closely spaced impacts. There are additional irregularly shaped crater depressions in the vicinity.

The Palmer Canyon Road crater area lies ~10 km SSE from the FRX crater field, ~3–4 km northeast of the Wheatland Reservoir Number 3 at the southern slopes of the Laramie Mountains. The Archean basement uplift is overlain by middle and upper Pennsylvanian and Lower Permian Casper Formation that gently dips here at 2–3° to the south. The Casper Formation consists of quartzitic white sandstone, but carbonates are present as well. The possible craters PCR-1 (42 × 40 m, R = 1.05) and adjacent PCR-2 (28 × 19 m, R = 1.47) are morphologically the most conspicuous crater structures (Fig. 3G) in close proximity to the Opeche Shale Member that is less obvious in color here. In total, eight possible structures have been found in the PCR area. PCR-001, a morphologically less obvious crater structure, contains several shocked quartz grains (Supplemental Material). Here, the elevated rim is composed of quartzitic breccia with microcrystalline chert fill. Chert is wind polished and scoured, and it looks like molten wax that forms cone-like morphologies (Fig. 4C). More breccia occurrences in the surroundings are interpreted as ejecta from the crater (Fig. 4E).

The country rocks at all craters are quartz arenites of the uppermost Casper Formation. At Sheep Mountain (SM), the sandstones are relatively pure in composition and comprise 95% quartz with a unimodal grain size distribution. The content of feldspar and accessory minerals increases westward. At Mule Creek (MC), tabular hematite-rich sandstones dominate with a bimodal grain-size distribution. At Fetterman Ridge (FR), the target sandstones are brighter, the quartz content is higher, and the mean grainsize distribution is unimodal, and cross-bedding stratification is ubiquitous. FR and SM samples appear more compacted than the MC samples. Commonly sandstones are sealed by quartz cement. Carbonate-filled joints and pore spaces were found mainly in the samples of the FRX area; carbonate layers occur in the PCR area. Chalcedony is frequent in brecciated samples (Figs. 3D–3E), as matrix material filling joints, in small fault zones, and it also occurs as massive layers showing flow textures. Chalcedony is further present in some samples as cement.

The Early Permian paleo-environment in SE Wyoming shows a transition from non-marine conditions in the SW toward shallow marine conditions in the NE (Steidtmann, 1974). The described lithofacies change in the study area may reflect a transition from an eolian (FR) to fluvial braided stream (MC) environment in the west to a beach and shoreface environment (SM, WR) in the east. In the vicinity of the craters, the uppermost Casper sandstone can form a blocky ejecta blanket surface. Such ejecta blocks are commonly subrounded and of decimeter size. The presence of (subsurface) water at the time of impact caused cohesion of sand but did not lead to a complete degradation of crater rims and ejecta blankets by any sort of water-assisted flow or mass movement. The facies change from the uppermost Casper sandstone to the overlying Opeche Shale Member red beds of the Goose Egg Formation is abrupt at all localities. Where exposed, the immediate contact is variegated and laminated. The first laminae are sandy and silty, and then the strata show interbedding of red and ocher mudstones and siltstones, with thin limestone and gypsum stringers indicating that they were deposited in non-marine saline lakes, pans, and lagoonal mud-flats (Benison et al., 1998).

Brittle deformation is evident both on macroscopic and microscopic scales. Casper sandstone is tectonically jointed (Fig. 3C) and contains cleavage fractures of tectonic origin. Portions of the crater walls and the preserved ejecta blankets show breccia textures that are sealed with quartz or chalcedony (Figs. 4B and 4E). However, breccia clasts appear subrounded. Occasionally, centimeter-to decimeter-sized dikes filled with clastics were observed in the crater walls. The degree of brittle deformation varies strongly. On the grain scale, many quartz grains are intensely fractured. Hertzian-type cracks form along stress chains at grain-to-grain contacts (Fig. 5D). Those quartz grains lost their original roundness upon impact loading, developed angular corners, and indented grain boundaries. Comminuted grains fit tightly against close-pressing neighbors. Many grains show irregular intra-clast-fractures that contain abundant fluid inclusions. These fractures terminate at the over-growth seams and indicate that they were formed when the rock was unconsolidated. More regular, trans-granular fractures that also fracture the cement between grains were formed when the rocks were consolidated.

Generally, shock features are relatively rare in samples taken from the craters. We did not find more than two or three grains within one thin section that contain definite shock evidence. Some of the shocked grains are isolated and surrounded by completely undeformed grains; others occur in intensely deformed parts of the rock. Table 1 and the Supplemental Material give an overview of the craters that contain shocked grains. So far, shock features have been detected in 31 craters and 71 samples. Almost all samples were taken in situ from different locations with respect to the given crater (for exact sample locations, see Supplemental Material). Shocked grains were sampled downrange, uprange, and crossrange along the inner crater wall, the crater crest, and on the ejecta flap. They were found in breccia and dike samples. A few shocked grains occur outside craters in more distal ejecta deposits a few tens of meters away from the craters. At those sites, shocked grains are commonly isolated in otherwise undeformed sand. Local reworking by wind cannot be excluded in those cases. We found both planar deformation features (PDFs) (e.g., Figs. 5A–5B) and planar fractures (PFs) (Fig. 5C) in quartz grains of various crystallographic orientations. Figure 6 displays the angles between the poles of PDF lamellae and PFs to the respective c-axis of the quartz crystal that contains the shock feature. PDFs are more frequent than PFs. PDFs typically have a spacing of a few micrometers; PFs have a larger spacing of 10–20 μm. PDFs in the basal orientation (0001) are most frequent, followed by {1013} and {1014} orientations. Several shocked grains show multiple sets of PDFs. All PDFs and PFs are decorated with fluid inclusions. Where a second set is missing, the crystallographic indexing does not always allow unique indexing. Based on calibration of shock effects, shock pressures below 10 GPa are indicated by (0001) orientations. However, some of the shock features indicate local pressure excursions above 10 GPa. The {1013} orientation, which is the second most common observed set (Fig. 6), indicates shock pressures >10 GPa, and the {1012} orientation indicates shock pressures >20 GPa (e.g., Grieve et al., 1996).

Cherty-textured silica is abundant at the various crater sites (Fig. 4D). Occurrences in veins and fillings of fractures are obviously related to diagenetic sealing. However, associated with crater SM-28 and others, we found cherty layers with a great variety of shapes, some with elongated forms and some wavy layers with flow textures that are interlayered with sand, and some form lumps that stack on top of each other (Fig. 5E). Some layers have sharp edges to the sandy country rock and resemble glass shards. Some of the cherty layers and lumps contain abundant spheres that have a particle in their center that is surrounded and coated by concentric rings of various grain sizes of fine-grained quartz. These objects of ~1 mm diameter fulfill the characteristics of accretionary lapilli (Fig. 4F) that are known from explosive volcanic eruptions (e.g., Schumacher and Schmincke, 1995) as well as from distal ejecta from large impact structures (Schulte et al., 2010). They form as diffused or suspended material in turbulent hot air enriched with particles that carry liquid films of condensed moisture. The lapilli grow by collecting fine particles that stick to the liquid on their surface when the binding forces exceed the grain dispersive forces. In impacts, Johnson and Melosh (2014) argued that molten silicate acts as the binding agent, or any other condensable material in the ejecta curtain could act as the binding agent. The deposits of ejecta curtain layers of large impact structures, such as the Chicxulub impact structure in Mexico, are composed of solid ejecta along with melt droplets, melt fragments, and accretionary impact lapilli (Schulte et al., 2010).

There are many crater structures with distinct ellipticities that can be used for the reconstruction of the impact direction (Figs. 1–2; Supplemental Material). The orientation of crater chains also allows constraint of the flight direction. Moreover, some of the craters possess structures that define up-range and downrange. Criteria (Gault and Wedekind, 1978; Kenkmann et al., 2014) for this are (1) a steeper crater wall uprange; (2) a preserved overturned ejecta flap downrange, with a preferred deposition of ejecta downrange; (3) an ovoid crater shape with the strongest curvature downrange, and (4) V-shaped herringbone patterns of ejecta pointing up-range. We used the orientation of the long axis of the craters for trajectory reconstruction when the ellipticity of a crater is larger than 1.2 and the outline of the crater is clearly defined. The extended crater field at Sheep Mountain has 10 elliptical craters and one crater chain, all of which consistently trend NW-SE (Fig. 3, Supplemental Material). The Mule Creek (MC) and Wagonhound (WR) areas each have two sufficiently defined craters with aspect ratios of at least 1.2 to derive a trajectory. Craters at Fetterman Ridge (FR) are more strongly degraded with tectonic jointing and incomplete rim preservation complicating the definite measurement of the trajectory direction. The Fetterman Road (FRX) location has one crater chain and one elliptical crater (Fig. 3F) to derive the impact direction. The southernmost area at Palmer Canyon Road (PCR) contains two elliptical craters that are suited for constraining the impact direction.

The precision of trajectory determination has an error of only a few degrees if the crater ellipticity is 1.4 or so but increases to as much as ± 5° if the crater ellipticity is 1.2 (Fig. 2). A tectonic tilt correction was applied to the trajectories. The Casper sandstone beds are rotated back to the horizontal plane around an axis that corresponds to the strike direction of the beds. As the dip of beds is commonly not more than 15°, the change in direction of the trajectories is usually less than 1°. The trajectories have not been corrected for Laramide tectonic translation and rotation because they bear too much uncertainty.

We used the suitable craters to define a fan-like corridor of trajectories for each crater field (Fig. 1). The opening angle of the corridor is given by the scatter of measurements in a particular crater field. Four such corridors have been derived (Fig. 1). Tracing back the trajectories allows the location of the primary crater to be estimated. We suspect that the primary crater is located in the area of the intersection of the corridors. The intersection of the trajectory corridors defines a polygon whose center is located at 41°28′N and 103°59′W. This polygonal area is interpreted as the most likely position of the source crater. However, tectonic shortening or rotations of the secondary crater locations may have considerably modified or shifted its actual position. Using the center of this polygon, the distance to the crater field SM is roughly 180–190 km, to WR 180 km, to MC 175 km, to FRX 165 km, and to PCR 155 km. Thus, all discovered secondary craters occur at a distance of 150–200 km from the proposed primary location.

We calculated ballistic paths of ejecta, taking into account aerodynamic drag and Earth´s curvature. Trajectories of ejecta with 1 m, 2 m, and 4 m radius were modeled with ejection angles ranging from 30° to 60° and initial speeds of 1 km/s, 2 km/s, and 4 km/s (Fig. 7), which is a reasonable assumption if the ejecta π-group scaling relationships of Housen and Holsapple (2011) are taken into account. Results are presented for the range, maximum altitude, impact energy, fraction of ejection energy, impact speed, and impact angle of the ejecta. The squared dependence of atmospheric drag on ejecta speed and the diminution of atmospheric density with altitude causes large and complex effects. The range curves show that the modeled ejecta would impact at distances between 10 km and ~700 km from the impact site, with the majority at a distance below 256 km (Fig. 7A). Many of the sets of ejecta radii, ejection speeds, and ejection angles used in this study resulted in ejecta paths that remained in the denser part of the atmosphere. For example, the smallest, slowest boulders ejected at 30° above horizontal reach at only 3.47 km altitude. The largest and fastest ejecta at 60° escape the atmosphere and travel long distances with only gravity significantly decelerating them along most of their paths (Fig. 7B). The 1 m radius ejecta all hit the ground with less than 1 GJ of energy. The 2-m-radius boulders return with between 5 GJ and 50 GJ, and the 4-m-radius ejecta return with impact energies between 64 GJ and 1000 GJ (Fig. 7C). The ratio of impact energy to launch energy is larger for the smaller starting velocities and also increases with size of the ejecta. The smaller boulders retain only ~1 percent of their initial energy because of their unfavorable ratio of surface to mass. The largest and fastest boulders in the set analyzed retain up tô30 percent of their initial kinetic energy (Fig. 7D). The impact velocity curves show that the 1-m-radius boulders return at about the speed of sound (Mach 1, 343 m/s at 20 °C and 1 atm). These are not likely to make a significant crater. Impact velocities of 1 km/s are reached for 4-m-radius boulders ejected at 2 km/s and for 2-m-radius boulders ejected at 4 km/s at an angle of 60° (Fig. 7E).

The impact angle curves show that the magnitude of that angle is larger than the initial ejection angle in all cases due to atmospheric drag and deceleration (Fig. 7E), and the impact angles are much larger for the smaller, slower boulders. The largest boulders impact at only slightly larger angles than their ejection angles, as would be expected in cases where atmospheric drag is a relatively small part of the deceleration along the path. For the larger boulders, the difference between the magnitudes of ejection and impact angles decreases with both speed and launch angle. The smaller, lower launch angle and slower boulders show a more complex impact angle behavior (Fig. 7F).

The ballistic range of 150–200 km from the primary, where we found all secondary craters, is not hit by 1-m-radius ejecta that launched at an ejection velocity of 1–4 km/s. Of the 2-m-radius ejecta, only those launched at 3 km/s or 4 km/s reached the distance where the craters are observed (Fig. 7A). Larger blocks of 4-m-radius ejecta that left the primary at a speed of 2 km/s also fall in this ballistic range. The impact velocities of these ejecta range from 500 m/s to 1000 m/s and depend on the ejection angle (Fig. 7E). The impact energies of these ejecta range from ~12–400 GJ (Fig. 7C). The dependency on the ejection angle is moderate, but higher ejection angles tend to give larger impact energies. The ratio between impact energy and initial launch energy ranges from 10% to 25% for ejecta landing at the range where the secondary craters have been found, while their impact angles range from 45° to 60°.

To link the impacting ejecta with crater formation, π-group scaling gives a first reference (Holsapple, 1993). Ejecta of 2500 kg/m³ density and 4 m radius impacting at 1000 m/s form craters of ~45 m diameter. Ejecta of 2 m radius would form a crater of roughly 25–30 m in diameter considering the obliquity of the impact. These dimensions are good representatives of the size of the mapped craters in Wyoming (Fig. 2). A more sophisticated approach to calculate the dimensions of the resulting impact craters uses the multi-material shock-physics code iSALE (Melosh et al., 1992; Ivanov et al., 1997; Wünnemann et al., 2006). We determined the expected size of the craters formed by projectiles of 1 m, 2 m, and 4 m radius, traveling at velocities between 500 m/s and 1500 m/s. Additionally, we determined the mass of material that experiences shock pressures greater than 3 GPa, which represents the minimum pressure required for the formation of shock metamorphic features in quartz (Grieve et al., 1996). All simulations were conducted using iSALE-2D, and as such, all impacts were limited to a vertical trajectory. Models included a non-porous and a porous sandstone target. Details of strength properties, equation of state, and model parameters are explained in the Materials and Methods section and in Table 2. To consider the effect of impact obliquity on the results for shocked mass, we scaled our results for crater diameter by a factor of sin^1.3θ (Schmidt and Housen, 1987; Pierazzo and Melosh, 2000). Crater diameters can be scaled by a factor of sin^0.44θ (Chapman and McKinnon, 1986) to consider impact obliquity; however, as the crater diameter is only expected to decrease by a factor of 0.86 for the most oblique case (45°) considered here, we do not directly scale our results for crater size.

The results of the iSALE simulations demonstrate that impactors of 1 m, 2 m, and 4 m radius, traveling at velocities between 500 m/s and 1500 m/s, are capable of generating craters that are 8–55 m in diameter (Fig. 8A). The depth-to-diameter ratio of these craters ranges from ~0.125 for 500 m/s impacts up to ~0.25 for 1500 m/s impacts. Impacts at 500 m/s generate no material shocked to pressures greater than 3 GPa. Impacts between 500 m/s and 1000 m/s generate between 0% and 10% of shocked material of the impactor's mass, and impacts at velocities greater than ~1000–1250 m/s (depending on target material) will generate significant, though still small, masses of shocked material (>10% of the impactor's mass) (Fig. 8B).

We discovered several crater clusters in Southeastern Wyoming located in the upper-most Casper Formation in transition to the Opeche Shale Member of the Goose Egg Formation. Due to their stratiform occurrence, we infer that all craters were formed at the same time ~280 m.y. ago in a single event. Secondary cratering appears to be the most plausible mechanism to explain clustering and crater characteristics. Observations that support the secondary cratering model are: (1) The elliptical crater morphologies allow the reconstruction of trajectories that meet in a single area. (2) The relative abundance of elliptical crater morphologies and partly irregular crater shapes suggest relatively low impact velocities, which is compatible with secondary crater formation. (3) The presence of radial crater chains and crater clusters that we observe is also known from secondary craters that formed on the Moon or Mars. (4) Commonly, small craters a few decameters in size are associated with relics of iron meteorites (Gottwald et al., 2020). The absence of meteorites in the Wyoming crater field is compatible with an ejection process from a primary crater as the cause of crater formation. As we did not observe any foreign rock types in the craters, we assume that the ejecta are also mainly composed of quartzdominated rocks such as sandstones or other felsic rocks. Ejecta boulders that cause secondary craters are believed to form by spallation from the primary crater (Melosh, 1984). This rules out that the source material is excavated from a deeper crustal level. (5) Accretionary lapilli cannot form in small craters because the shock level is insufficient to either form melt of a significant volume or to develop a hot plume above the small craters. The presence of possible accretionary lapilli mixed with the ejecta of the small craters indicates that a presumably hot ejecta plume existed and interacted with the local ejecta distribution. (Fig. 8B). (6) Linear ejecta accumulations (Figs. 2A–2B) are interpreted as the remnants of possible herringbone patterns formed by the interaction of ejecta of the primary and secondary cratering process. (7) The aforementioned calculations of impact trajectories reasonably match the conditions for the formation of craters ranging from 10 m up to 55 m at distances of 150–200 km from the proposed primary crater. (8) The small crater volumes exceeding a 3 GPa shock level are in agreement with the limited amounts of shocked minerals found in the craters.

The proposed location of the primary crater shown in Figure 1 is in the border region between Wyoming and Nebraska involving Goshen and Laramie Counties in Wyoming, and Kimball, Banner, and Cheyenne Counties in Nebraska. Tectonic shortening or large-scale block rotations during the Laramide orogeny would shift and rotate the inferred trajectories and would widen the possible area of the location of the primary crater to the entire Northern Denver basin, where 280 m.y. old Pennsylvanian-Permian strata are deeply buried beneath Triassic to Neogene beds. The Denver basin is asymmetric, bordered by the Front Range, the Laramie Range, and the Hartville Uplift to the west, and the Black Hills to the north. The structural axis of the basin trends approximately north-south from Denver, Colorado, to Cheyenne, Wyoming (McCoy, 1953; Hoyt, 1962). Stratigraphic equivalents to the Casper Formation are the Ingleside Formation near the Wyoming-Colorado border, the Fountain Formation in the Denver basin, and the Tensleep Sandstone in the northern part of the basin. In the area of interest, the equivalent unit is mapped as upper and middle Minnelusa. The basin has been intensely explored for hydrocarbons with many boreholes. The top of the crystalline rocks in the region of Wyoming and Nebraska was mapped by Cardinal and Holmes (1984) and Blackstone (1993). Minnelusa strata were found in 39 deep oil well logs, and most of these wells were deep enough to also hit the Precambrian that occurs at varying depths down to more than 3000 m. Based on a first superficial inspection of well log data, so far we have not found distorted or unexplained missing sedimentary rock sections that could give a hint at the location of the primary crater. However, the 1–35 Hawk Fee core shows some breccia layers at 3023–3066 m. We are planning to visually inspect and sample critical cores when the U.S. Geological Survey Core Research Center (Core Research Center, 2021, www.usgs.gov) re-opens after the current pandemic situation.

Does the size distribution and range of ejecta allow us to narrow down the dimensions of the source crater? Excavation dynamics of impact cratering suggest that the largest fragments ejected are most likely spall plates that originate from near the surface of the primary crater (Melosh, 1984; Vickery, 1986). The spall size scales with the radius of the primary impactor and the target tensile strength and is inversely proportional to ejection velocity. The shock level in the spalls is considered low. The spall size depends on rock mechanical properties such as target density, crushing strength, Poisson ratio, tensile strength, and Weibull parameter of the material's flaw distribution (Melosh, 1984). Impact parameters such as the impact velocity and depth of penetration scale with the projectile radius (Melosh, 1984). We considered an impact velocity of 20 km/s for the primary projectile and strength parameters listed in Melosh (1984; Table 2, case 5) and used Melosh's relationship of fragment size to ejection velocity (Melosh, 1984; Fig. 8) that is also supported by Vickery (1986) and Hirata and Nakamura (2006). Then, for the most likely ballistic scenarios (4 m radius ejecta launched at 2 km/s and 2 m radius ejecta launched at 4 km/s), we obtained an estimated projectile radius of the primary projectile that is in the order of 2.0–2.7 km.

This in turn relates to a primary transient crater diameter of 32–40 km using the scaling relationship of Schmidt and Housen (1987) and a final crater size of 50–65 km using the scaling of McKinnon and Schenk (1985) if standard impact conditions are considered (impact velocity: 20 km/s, projectile density: 3000 kgm³, Target density: 2950 kgm³, impact angle: 45°).

Using the USGS National Geological and Geophysical Data Preservation Program (2021; https://www.usgs.gov/core-science-systems/national-geological-and-geophysical-data-preservation-program), we recognized that the proposed area of the primary crater in the northern Denver basin shows a negative free-air gravity anomaly of up to -60 mgal with uneven relief. The overall shape of the gravity low appears elliptical with a NW-SE extent of ~130 km and a NE-SW extent of ~65 km in that region. However, this gravity low appears constricted in the middle part, so that it could be described by two roughly circular negative gravity anomalies of ~65 km diameter each, centered at 41°50′N/104°10′W and 41°40′/103°30′W, respectively. The northwestern gravity anomaly is framed by the arcuate Hartville High to the NW and the arcuate Morrill County High to the SW. This location has coincident missing Minnelusa Formation and a 25-m-thick salt accumulation associated with the Opeche Shale (Oldham, 1996). An in-depth search for the primary crater location and its site characterization is reserved for a future study and is beyond the scope of this work.

Using environmental scaling relationships of impact cratering (Collins et al., 2005), a 50–65 km crater would cause a fireball in the order of 60 km radius and would lead to a considerable thermal exposure of ~4 × 10⁸ J/m² with a duration of more than 10 min where the secondary craters have been found. The observed potential accretionary lapilli are strong indicators for a hot ejecta plume that expanded laterally by ~150–200 km, measured from the proposed location of the primary. Accretionary lapilli of 1.5 mm average size have been described from Albion Island, Belize, at 360 km distance from the Chicxulub source crater (Pope et al., 1999). They are exclusively known from large impact events. Johnson and Melosh (2014) showed that the accretionary impact lapilli size d_l is a function of the impactor radius R_imp and ejection velocity v_ej. According to the equation,

it is suggested that the ejection velocity of these particles has been 2.6–2.7 km/s for the assumed projectile radius of 2.0–2.5 km. At Albion Island, accretionary lapilli occur together with altered impact glass that is degraded to smectite and palagonite (Pope et al., 1999). In the Wyoming crater field, we found that the millimeter-sized possible accretionary lapilli are entrained in microcrystalline chert lumps and shard-like cherts that show abundant flow textures. If the target of the primary was dominated by silica we suggest that the chert lumps could have been silicarich melt that was transported with the ejecta along with the accretionary lapilli. The shapes of the chert lumps vary; some are sub-angular and some are rounded, which suggests that deposition occurred around the glass transition temperature, so that in some cases shard-like chert aggregates could form. We interpret that the material was glassy upon landing. However, it is unclear whether these chert formations are a primary feature or related to diagenesis.

Initially the Sheep Mountain (SM) crater field near Douglas, Wyoming, was interpreted as a large 7.5 × 1.5 km impact crater strewn field formed by violent disruption of a single meteoroid during passage through the atmosphere (Kenkmann et al., 2018). However, the size of the proposed strewn field was already unusually large in comparison to the other terrestrial impact crater strewn fields, Morasko, Odessa, Wabar, Henbury, Sikhote Alin, Kaalijärv, and Macha, which are all densely clustered (Kenkmann et al., 2018; Gottwald et al., 2020). The disruption during atmospheric traverse does not allow a separation of fragments perpendicular to the trajectory by more than 1 km (Artemieva and Shuvalov, 2001). Moreover, all other crater strewn fields are associated with iron meteorites. The scenario of a break-up of a single meteoroid is incompatible with the current size of the crater field and the other new findings and must therefore be rejected. Likewise, multiple airbursts of a single meteoroid at a very high altitude could not explain the wide distribution of small craters. A paired meteoroid, a meteoroid shower, or a tidal break-up of an asteroid prior to atmospheric entry could explain the wider extent of the crater field but not the converging trajectories and the lack of meteorites that we observe in the Wyoming crater field. An alternative explanation could be a collision in the asteroid belt that led to a phase of enhanced bombardment on Earth, similar to the Late Ordovician event (Schmitz et al., 1996; Schmieder and Kring, 2020). However, this would have led to a global increase in impact cratering rather than a local concentration of small craters. This scenario is not supported by our observations, at least at the current state of investigation.

The craters experienced a complex, 280-m.y.-long post-impact history. Sedimentation of the craters beneath shales of the Goose Egg Formation must have occurred shortly after the impact event in a low-energy depositional environment, such as a sabkha or lagoon. Otherwise, the craters would have been degraded more strongly. Progressive burial caused diagenetic cementation and compaction of the sands that transformed the sands to quartzitic sandstones. The heavily deformed rocks of the craters and their ejecta peripheries were subject to substantial sealing by quartz overgrowth during diagenesis, which made these rocks more resistant to later erosion and let them appear as pedestal morphologies. From the regional geologic context, we estimate that the craters have been buried by 2000–3000 m of late Paleozoic and Mesozoic strata. The Laramide deformation took place from 70 m.y. to 50 m.y. (Fan and Carrapa, 2014). The amount of rock uplift of individual Laramide ranges was ~0.21–0.32 mm/yr from the early Maastrichtian to the early Paleocene and increased to more thañ0.38–0.60 mm/yr during the late Paleocene–early Eocene (Fan and Carrapa, 2014). All major Laramide uplifts including the Laramie Mountains experienced a minimum rock uplift of 3 km in the later phase, associated with the tilting of Pennsylvanian-Permian beds. The Laramide ranges reached a high elevation (> 2 km) at this time, and basement rocks were widely exposed to erosion during this period. The Pennsylvanian-Permian Casper Formation began to erode in the Paleocene along the axial zone of the uplift. Since then, the erosion front slowly retreated to the surroundings of the basement uplifts, thereby exposing fresh parts of the impact event strata. This erosion process is still continuing. However, in some parts of the Sheep Mountain area, the Eocene-Oligocene White River Formation overlies the uppermost Casper Formation.

Tectonic overprinting of the region during the Laramide orogeny is obvious from strata tilting, transpressive shear zone formation, and jointing, but also from fracturing and the formation of Boehm lamellae on the grain-scale. However, a clear distinction between brittle deformation of tectonic origin and brittle deformation caused by impact is often possible on the grain-scale. We found that impact-related fracturing and shock features affected the rounded quartz grains but not the quartzitic overgrowth seams of the quartz grains (Kenkmann et al., 2018) (Figs. 2A–2C). This implies that the impact occurred in weakly consolidated and locally unconsolidated sand prior to diagenesis. In contrast, tectonic deformation typically affected fully consolidated sandstones so that the overgrowth cement also shows brittle deformation overprint in the form of transgranular and intergranular fractures.

The state of preservation of the explored craters varies considerably. While some impact craters close to the contact of the overlying Opeche Shale Member are preserved superbly, others are almost completely eroded, lack any breccia, and show very minor indicative structures such as halos and radial fracture networks. The top of the Casper Formation has been closely examined only at those sites where crater structures and ejecta were encountered. We did not systematically prove how widely shocked quartz grains occur stratiform beyond craters and their ejecta. Wind and water may have reworked and transported shocked quartz grains shortly after the impact event. A mixture of shocked quartz grains from the primary crater and from the secondary craters in the ejecta cannot be ruled out. The role of chert that occurs frequently around the craters and the possible accretionary lapilli need further geochemical analysis to unravel their origin. Our ballistic modeling includes atmospheric drag at standard atmospheric conditions. We considered only spherical ejecta and did not account for break-up effects. Modifications of atmospheric drag by thermal and blast effects and different atmospheric pressure during the Permian were neglected. The spatial resolution of the numerical model of secondary crater formation limits the resolution of local pressure excursions that are inferred from some of the shock effects. So far, the primary crater formation has not been modeled numerically.

Thirty-one impact structures could be confirmed and more than 60 additional possible craters have been identified along the narrow band of outcrops of the uppermost Casper Formation in transition to the Opeche Member shales. The inferred impact trajectories define a 30–45° corridor pointing to the possible primary crater. Finding more secondary craters in a wider range may help to better localize the primary crater. High-resolution reflection seismic studies may help to identify and visualize more crater structures that are still buried beneath the Opeche Shale Member of the Goose Egg Formation. The most significant objective for future studies, however, is searching for the primary crater. We propose a comprehensive geophysical survey of the deep northern Denver basin to focus possible deep drilling to prove the site of the primary crater. We encourage the hydrocarbon exploration industry operating in that region to actively support and accompany the search for the primary crater and to report unusual occurrences of breccia.

The discovery of the Wyoming impact crater field allows some general conclusions to be drawn even if the secondary crater field is not completely recovered. The relatively dense atmosphere and the gravity field limit the extent of secondary crater fields on Earth. Owing to atmospheric drag, ejected material should be subjected to a more intense ballistic fragmentation and a stronger deceleration than ejecta on other planetary bodies such as Mars or the Moon. Consequently, secondary craters on Earth show a strong clustering. Due to low impact velocities, many secondary craters on Earth are elliptical in outline and shallower than primary craters. The reliable identification of impact craters requires evidence of shock effects or meteoritic material. The latter is not present in secondary craters. The detection of minute amounts of shock effects is only possible if the ejecta blocks have a minimum impact velocity of ~700–1000 m/s and a diameter of at least 4 m. This implies that secondary craters can only be found in a rather limited ring around the primary crater, where these conditions are met.

Our investigation is based on remote sensing using drone and satellite imagery, geological field work, microscopic analysis, and modeling. We used a DJI Phantom 4 drone for high-resolution mapping of crater clusters. Maps have been compiled and stitched from separate drone flights using ESRI Drone to Map software and ARCMAP GIS software. We also used geo-referenced satellite images provided by the Base-Map function of ArcGIS 10.7.1 by ESRI© with spatial resolution of ~0.3 m and Google Earth Pro© for identification of crater structures and mapping purposes.

Polished thin sections of rock samples were prepared and examined with a Leica/Leitz DMR polarizing microscope equipped with a digital camera. Shock lamellae were investigated at 200× and 500× magnification. Orientations of PFs and PDFs were measured with a Leitz U-stage mounted onto a Leitz polarizing microscope. Measurements were indexed using a template for crystallographic orientations in quartz. The accuracy of U-stage measurements is estimated at ± 5°.

We used Newton's Second Law in polar coordinates from Earth's center to calculate the ballistic paths of ejecta. Polar coordinates led to simplification where the lengths of paths of ejecta covered appreciable curvature of Earth. Further, gravity always acts in the radial direction in this coordinate system. We computed

where F_r and F_θ are the forces in the radial and circumferential directions, respectively, m is mass, r is the radial coordinate, and θ is the circumferential coordinate. Each dot represents a time derivative.

Time increments in the computation varied from 0.5 seconds when far above Earth's surface to one millisecond, with the duration proportional to gravitational acceleration divided by the sum of the magnitudes of the gravitational acceleration and the acceleration due to atmospheric drag. In program development, various time increment strategies were evaluated to ensure convergence and stability. The aerodynamic drag was calculated from the Bernoulli drag equation:

where the force (F) equals half of the fluid density (ρ), times the square of the object's velocity (v), times the frontal area of the body (A), times the coefficient of drag (C_d). In this study, many of the modeled ejecta velocities were near the speed of sound, 340 m/s. The coefficient of drag was calculated for each time increment using an algebraic equivalent to the coefficient of drag in the equation stated in Carter et al. (2009). That is:

where M is the Mach number of the velocity, taken as the speed divided by the speed of sound. The density of the ejecta was assumed to be 2700 kg/m³, which is a reasonable approximate value for sediments or granite. It was assumed that the ejecta boulders are spherical. The density of the atmosphere near Earth's surface was set at 1.204 kg/m³, corresponding to 20 °C. We have not considered possible different atmospheric composition and pressure during the Permian. The exponential decrease rate of atmospheric density was, therefore, 1/8542 m^-1. The integrated drag on ejecta that reach above 50 km altitude is negligibly affected by these estimates if the calculated resultant pressure at Earth's surface is one atmosphere.

The resultant craters were modeled using the iSALE shock physics code, an extension of the SALE hydrocode (Amsden et al., 1980). Notable modifications to the original SALE code include: an elasto-plastic constitutive model, fragmentation models, various equations of state, and multiple materials (Melosh et al., 1992; Ivanov et al., 1997). More recent modifications include a modified strength model, a porosity compaction model, and a dilatancy model (Collins et al., 2004; Wünnemann et al., 2006; Collins, 2014).

Simulations in this study varied the size and velocity of an impactor with an ANEOS (Thompson and Lauson, 1972) equation of state for quartzite (Melosh, 2007) and a reference density of 2650 kg/m³. In all simulations, the impactor was non-porous. The target was composed of one of two materials. First, simulations were run with a “non-porous sandstone,” that is, a material with an equation of state for quartzite (Melosh, 2007) but with the strength parameters of a sandstone (Güldemeister et al., 2015). Porosity φ has a strong effect on the attenuation of shock pressures within a material (Wünnemann et al., 2006, 2008). Consequently, a second target material was modeled a “porous sandstone,” with the same equation of state and constitutive parameters of the “non-porous sandstone” and a porosity of 20% (Wünnemann et al., 2006; Güldemeister et al., 2015). While the target materials of the Wyoming crater field were likely water-saturated to some extent, we note that water saturation will produce an intermediate response between a dry porous and non-porous material (Pierazzo et al., 2005; Güldemeister et al., 2013), and therefore our results provide minimum and maximum boundaries on shock masses for the target conditions of the Wyoming crater field. The parameters used in the modeling are shown in Table 2.

We acknowledge thorough and thoughtful reviews by Boris Ivanov and Axel Wittmann and editorial work by Uwe Reimold. We are grateful to Dagmar Flemming, Gordon Mette, and Herbert Ickler of the University of Freiburg for preparing hundreds of thin sections in the course of this project. Michael Poelchau assisted us on one of the field trips and Tabea Schulze analyzed thin sections as part of her B.Sc thesis. We thank numerous Casper College geology students for their field documentation of craters, drone image work, and rock preparation. We gratefully acknowledge the developers of iSALE, in particular: Jay Melosh, Boris Ivanov, Gareth Collins, Kai Wünnemann, Dirk Elbeshausen, and Thomas Davison. T. Kenkmann is grateful to Michael Poelchau and Gerwin Wulf for inspiring discussions of the topic. We acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG), project KE 732/30-1.