Popigai and Chicxulub craters: multiple impacts and their associated grabens

More advanced data (gravity field model EIGEN 6C4 including the GOCE gradiometry data instead of EGM 2008) and a more sophisticated method (using a set of gravity aspects instead of gravity anomalies and the radial second derivative of the disturbing potential only) enable a deeper study of various geological features. Improved techniques were applied to study the impact craters Chicxulub and Popigai. We confirm our results from 2010, extend them, and offer more complicated models, namely by means of the gravity strike angles. Both craters are interpreted to be double or multiple craters. The probable impactor azimuth was from NE (to SW) for Chicxulub and SE (to W) for Popigai. The formation of both the craters seems to be associated with impact-induced tectonics that triggered the development of impact grabens.

1 Motivation

In 2010, we published (Klokočník et al., 2010, in this journal) the results of our tentative analysis of the gravity data for two areas of proven, huge impact craters Chicxulub and Popigai. The analyses were based on the global combined gravity model EGM 2008 (Pavlis et al., 2008a, b, 2012) to degree and order $(\mathrm{d}/\mathrm{o})\mathrm{2159}$, which was the top for modeling of the gravity field of the Earth at that time for precision and resolution. We suggested (in Klokočník et al., 2010) that Chicxulub (northern Yucatán, México) may be a double crater and Popigai (northern Siberia, Russia) may be a multiple crater (see the labels “Popigai I–IV”, according to the Rajmon, 2009, in the catalog on © Google Earth). Craters II–IV are not yet proven to be impact crater candidates. Altogether they would create a hypothetical Popigai crater family, a catena (see the Supplement S3). We analyzed gravity anomalies and second radial components of the Marussi tensor. Now, we work with a set of gravity aspects (incorporating those two functions, too).

Since that time we have analyzed many geological features on the Earth, the Moon, and Mars; we make use of new gravity models and other data sources (see below), and we have summarized our results in three books and about 15 papers (e.g., Klokočník and Kostelecký, 2015; Klokočník et al., 2017, 2018, 2020a, b, 2021, 2022a, b, 2023a, b).

We (and the readers) are well aware that solely the gravity data are not unambiguous to detect ground density anomalies (causative bodies); we always need and seek additional data (geological, geophysical, seismic, topography, archeology). Therefore, with our data and method, we offer only a step toward possible field explorations and subsequent interpretations but not a confirmation of the structures.

With increasing precision, accuracy, resolution, and reliability of our knowledge about gravity and magnetic fields, we can test diverse applications impossible before. This is not only about data available; it is also about methodology. The traditional gravity anomalies are no longer sufficient. We apply a set of gravity aspects instead (Sect. 2 and Supplement S1).

Another impetus for this study was the similarity of the two craters in the sense that they are directly associated with close linear structures that seemingly have nothing to do with the impact event but occur close to the craters. This could be a coincidence of two different genetically independent geological phenomena or a sign of the existence of a trench being modified by the impact event. The impact shock affected the entire region, and previously existing faults or fault zones that were in an extensional tectonic regime were activated to form an impact graben.

This paper is a revival concerning our previous findings about Popigai and Chicxulub. Now we have (in comparison with Klokočník et al., 2010) better tools (the set of the gravity aspects) and a better gravity model (EIGEN 6C4, with GOCE data); thus, we can support or reject the older results. We offer new and hopefully more convincing results in favor of a double or multiple character of both the Chicxulub and Popigai craters. These results are not in conflict with the known geology of the areas or with information from magnetic intensities (not studied here; see, e.g., Hildebrandt, 1998, Urrutia-Fucugauchi et al., 2022, or Mendes et al., 2023).

Many figures which may help the reader are gathered in the Supplement in Si: there is a tutorial in S2, with tests about artifacts – S3 for Popigai and S4 for Chicxulub. The theory is shortly repeated in S1. Link: https://www.asu.cas.cz/~jklokocn//CHIC-POP24_supplements/ (last access: 11 February 2025).

2 Notes on theoretical preliminaries

The gravity (gravitational) aspect (descriptor) is a functional or function of the disturbing gravitational field potential T_ij. We work with the gravity anomaly (or disturbance) Δg, the Marussi tensor (Γ) of the second derivatives of the disturbing potential (T_ij), two gravity invariants (I_j), their specific ratio (I), the strike angles (θ), and the virtual deformations (VDs).

The theory came mainly from Pedersen and Rasmussen (1990) and Beiki and Pedersen (2010). The theory with examples is summarized in our books (Klokočník et al., 2017, 2020a, 2022b); it cannot be (due to space reasons) repeated here (see Supplement S1). Only a few notes follow.

The gravity aspects are sensitive in various ways to the underground density contrasts (variations) due to the causative bodies, exciting the relevant gravity signals. The set of gravity aspects tells us much more about the causative body than the traditional Δg only. It informs us about the location, shape, orientation, tendency of the ground structure toward 2D or 3D patterns, and stress trends and may simulate “dynamic information” about existing tensions (although the input data are always the harmonic geopotential coefficients of a static gravity model).

For example, the strike angle θ can mathematically be the main direction of the Marussi tensor Γ of the second derivatives of the disturbing potential (the first column and first row of Γ are identically equal to zero for this preferred direction). The strike angle is, from the geophysical point of view, a direction important for description of the ground structures. It may indicate areas with a lower density or higher porosity, a “stress direction”, or both, or it may indicate the areas under a strong influence of rapid and/or intense geomorphic processes. When I=0, the values of θ may be symptomatic of a flat causative body. For more details see Beiki and Pedersen (2010) and our S1.

A usual situation is that the strike angle θ has diverse directions, as projected on the Earth's surface. The combed strike angles are the strike angles oriented roughly in one and the same direction in the given area. The theory for the “combed” strike was explained, together with relevant statistics, in Klokočník et al. (2020a). For statistical use we defined a degree of alignment of the strike angles by the “comb coefficient” (comb) as a relative value in the interval 〈0,1〉. Zero (0) means that the strike angles are “not combed” (totally disheveled, the vectors θ are in diverse directions); 1 means “perfectly combed” (perfectly kempt, the vectors of θ are oriented into one prevailing direction). If the comb is smaller than 0.55, we say that θ_i of the given region is not combed; if comb>0.65, we say that θ_i is combed, and for comb>0.99, they are perfectly aligned. The alignment may take a linear form (e.g., along a fault) or can take the shape of a halo (around craters); see the Theory section in S1 and many examples in the tutorial (Supplement S2).

Link: https://www.asu.cas.cz/~jklokocn//CHIC-POP24_supplements/ (last access: 11 February 2025).

3 Data, computation, and figures

We always start all our computations with the harmonic geopotential coefficients (Stokes parameters) of the static global gravity field models as the input data; they describe the gravitational potential of the Earth. The whole theory is prepared in such a way that we cannot use input other than the harmonic coefficients (Pedersen and Rasmussen, 1990; Beiki and Pedersen, 2010; Klokočník et al., 2017, 2020a, 2022b).

We make use of the high-resolution combined European Improved Gravity model of the Earth by New techniques (EIGEN 6C4; Förste et al., 2014), expanded to degree and order (d/o) 2190 in spherical harmonics; this corresponds to the ground resolution 5×5 arcmin or ∼9 km on the surface. The precision of EIGEN 6C4, expressed in terms of Δg, is N=10 mGal, but in many civilized land areas and over the oceans and open seas it is much better. The creators of EIGEN 6C4 did not have access to most of the recent high-resolution terrestrial gravity data on the continents, and thus they took a synthesized gravity anomaly grid based on EGM 2008 (Pavlis et al., 2008a, b, 2012). That means that the errors for high d/o terms in EIGEN 6C4 are dominated by the relevant errors in EGM 2008. To estimate the actual realistic precision for the given area of interest – not only for a general figure of 10 mGal – one needs to inspect gravity anomaly commission error maps of EGM 2008 (Pavlis et al., 2008a, b; also in S3). They account for a complete covariance matrix between the solved-for harmonic coefficients in this gravity model. Using those maps for the northern Yucátan peninsula, we get N=4–8 mGal. For Popigai in Siberia, it is ∼15 mGal.

A note about other data sources: ETOPO 1 global surface topography (Amante and Eakins, 2009) is a global 1^′ relief model of the Earth's surface that integrates land topography and ocean bathymetry from a large number of satellite and other measurements. Its precision globally should be 10 m in height, and its accuracy is ∼30 m.

There are alternative topography data files (not fully independent of the ETOPO files), like various versions of ASTER GDEM (JPL NASA and Japan), SRTM (NASA), GEBCO, and the ETOPO 2022 release. Google Earth is also helpful.

Bedmap 2 is a subglacial topography valid for Antarctica (Fretwell et al., 2013). It contains the bedrock elevation beneath the grounded ice sheet. It is given as a 1×1 km grid of height of the bedrock above sea level, but actual measurements are often much sparser. We also worked with RET 14 (Hirt et al., 2016), the degree-2190 gravity field model SatGravRET2014, given as a set of harmonic geopotential coefficients meaningful only for the continent of Antarctica (not globally!). Roughly speaking, it combines the global gravity field model EIGEN 6C4 and the Bedmap 2 topography.

The data for magnetic analysis on the Earth are the grid value from the worldwide EMAG 2 model for magnetic intensities (Maus et al., 2009). There are also gravity field models, global topography, and magnetic data for the Moon and Mars; gravity field models of the Moon and Mars already provide sufficient ground resolution for our analysis. It is about 10 km for the Earth and the Moon but only 130 km for Mars. We present examples based on these gravity models in S2.

We computed the gravity aspects over many regions of the world in a step 5×5 arcmin in latitude and longitude, corresponding to the ground resolution of 9 km. But we can also use (and use here) a 4 km resolution without any degradation of the results (we offer some results of our truncation error tests and testing of artifacts in Klokočník et al., 2021, and here in Sect. 4 and S2). This higher resolution sometimes adds new and valuable information.

The numerical stability of computations of high degree and order functions in the aspects is extremely important; it was intensively investigated and tested and is guaranteed to much higher degree and order than we need here (work done during the last 10 years or so by the co-authors of this paper, plus Sebera et al., 2013, and Bucha and Janák, 2013).

Our figures are not generated by an automat but are created manually and individually with specific scales to emphasize various features and details. We plot all the quantities in geodetic (geographic) latitudes and east longitudes.

The gravity disturbances (anomalies) are given in milligals (mGal), and the second-order derivatives are in Eötvös (E). Let us recall that 1 mGal = 10⁻⁵ m s⁻² and 1 E ≡ 1 Eötvös = 10⁻⁹ s⁻². The invariants have units I₁ (s⁻⁴) and I₂ (s⁻⁶). The strike angle θ (°) is expressed in degrees with respect to the local meridian; red means its direction is to the east and blue to the west of the meridian. Often, we plot θ in black and white only. The strikes are shown in a regular grid 5×5 arcmin; this does not have any geophysical meaning and is just a choice for plotting.

4 Artifacts

4.1 Our previous work

To avoid various misinterpretations we need to test the input data to our analyses in various ways (e.g., Klokočník et al., 2021). We did our best to avoid the artifacts, but nevertheless, they cannot be excluded (S4: slide no. > 23). The important aspects are the resolution and statistical significance.

The ground resolution (GR) of the gravity field model is derived from the maximum d/o of its spherical harmonic expansion (the definition of GR is in S2: 23); more detail can be found in Sect 4.2.

Another important factor is the signal-to-noise ratio, $R=S/N$. The “signal” S is given as the range of gravity anomalies in the area of interest. The noise N is the commission error of the gravity anomalies Δg (see figure in S1, last page) or the estimated precision of Δg of EIGEN 6C4. We need R>3 to have statistically significant results. We have the following.

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\min R=\left(\min \right(\max |-S|,\max +S\left)\right)/\left(\max N\right)\text{(2)}& {\displaystyle }& {\displaystyle }\max R=\left(\max \right(\max |-S|,\max +S\left)\right)/\left(\min N\right)\end{array}$$

With the figure from S1 defining N and the figures below and in S3 and S4 defining S, for Δg of Popigai and Chicxulub, we get $R(\min ,\max )=\mathrm{8}$–15 for Popigai and $R(\min ,\max )=\mathrm{5}$–20 for Chicxulub.

4.2 Resolution

The reader certainly knows about the “canals” or “human faces” on Mars; they disappeared with new and better observations and higher resolution (S2: slide no. 25). The adequate GR of the gravity model is an important and necessary but not sufficient condition and a limiting factor for the correct interpretation of the gravity aspects. The definition of GR is recalled in S2 (23), and we can only repeat (Sect. 3) that the GR of EIGEN 6C4 is 9 km but can be enhanced to about 4 km (see above). This provides a clear limit for any interpretation. The subglacial topography has a similar problem: data gathered from airplanes over Antarctica (Bedmap 2) are not homogeneous in latitude and longitude, are not complete, and have large gaps (Fretwell et al., 2013). Taking the resolution of the subglacial topography data Bedmap 2 as published, i.e., net 1×1 km literally, we can get pictures of the topography showing unrealistic shapes instead of real features (Klokočník et al., 2021); the artifacts look like walls or pyramids (for example, S2: 28). The problem is that the data density is in some places ∼5 km but ∼50 km in other localities; there are zones with no data at all. We have to know how well the data have covered the area of our interest (Fig. 3 in Fretwell et al., 2013).

4.3 Signal degradation and truncation error tests

A treacherous situation with artifacts can be demonstrated by using the gravity model to its maximum degree and order (d/o) in harmonic expansion, exactly as it was published. The result may be surprising. For example, a model is published to $\mathrm{d}/\mathrm{o}=\mathrm{1200}$ but recommended for use (by the authors of the model themselves) only to max $\mathrm{d}/\mathrm{o}=\mathrm{600}$. The reason is stabilization of the large matrix inversion by the Kaula rule for the higher-degree part of the model. The full model can show significant graining in the gravity aspects, leading to total damage of the signal; see S2 (in all the gravity aspects – faster degradation was observed for the gravity aspects with higher derivatives of the disturbing potential; Klokočník et al., 2021).

Figure S2 (29) shows one of our many tests, in this case for the Moon's crater Copernicus with the gravity model GRGM1200A (Lemoine et al., 2014) to $\mathrm{d}/\mathrm{o}=\mathrm{1200}$. The practically useful limit at $\mathrm{d}/\mathrm{o}\sim \mathrm{600}$ corresponds to the theoretical ground resolution of ∼10 km. This is already comparable to the Earth, to its EIGEN 6C4 gravity model to $\mathrm{d}/\mathrm{o}=\mathrm{2190}$ (Förste et al., 2014), because the Moon is smaller than the Earth. When we use GRGM1200A up to $\mathrm{d}/\mathrm{o}=\mathrm{600}$, we can see a reasonable result (S2: 29) showing all known features. When we cut at $\mathrm{d}/\mathrm{o}=\mathrm{130}$, part of useful signal is lost. When we use the model to $\mathrm{d}/\mathrm{o}=\mathrm{1200}$, graining is significant and we can interpret nothing.

Let us imagine that today we know the gravity field of the Moon only to $\mathrm{d}/\mathrm{o}=\mathrm{10}$. What information do we lose (or is “hidden”) in a comparison with the full model to $\mathrm{d}/\mathrm{o}=\mathrm{600}$? Not only is the resolution of the former much lower (expected) but sometimes artifacts are also created (see S2: 30) (expected?). Only a further gravity field improvement would eliminate such artifacts. We are now in an analogical situation with the gravity field EIGEN 6C4 to $\mathrm{d}/\mathrm{o}=\mathrm{2190}$ for the Earth. What we would lose and which artifacts might be generated with, say, the model cut at $\mathrm{d}/\mathrm{o}=\mathrm{80}$? The slides in S2 (31, 32) show the result in terms of the strike angles. Often the basic trend in both the full model and the cut model is the same, but not always; thanks to the dramatic difference in maximum d/o used, it must be expected, but in any case, it is a warning. The artifacts “lurk” and can eventually hamper our endeavor concerning the geo-interpretations. It is not probable but not excluded even for EIGEN 6C4 to 2190; the case of artifacts due to an aliasing of the gravity aspects in the Sahara is in Fig. 5a in Klokočník et al. (2021) and S2 (33) here.

5 Popigai

5.1 Introductory notes and geology

This large, proven, exposed impact crater, Popigai/Popigaj, is in Russia near Khatanga (Chatanga, port on the river), Krasnoyarsk district, Siberia (geodetic latitude and longitude of the center of the crater: φ = 71°36^′ N and λ = 110°55^′ E). It is a 100 km diameter proper crater that is ∼35 million years old (from the late Eocene epoch). It was considered for the first time to be an impact crater by Masaitis et al. (1972), especially based on petrographic observations of the various breccias. It is the largest known impact crater post-dating the Cretaceous–Tertiary boundary (e.g., among many others, Vishnevsky and Montanari, 1999; Whitehead et al., 2000; French and Koeberl, 2010; Masaitis, 1998, 2003, 2019).

The impactor is suggested to have been an H chondrite asteroid several kilometers in diameter (e.g., Schmitz et al., 2015) from the main asteroid belt. The asteroids may have approached the Earth at comparatively low speeds, passed the Roche limit, and produced a meteoritic shower. But also a multi-type asteroid shower may have been recorded, triggered by changes in planetary orbital elements due to orbital resonances (see again, e.g., Schmitz et al., 2015). There is no agreement among researchers.

The Popigai crater lies on the eastern edge of the Archean Anabar Shield, which is mainly composed of granitoids and gneisses. It is surrounded by a relatively complex envelope of Precambrian, Paleozoic, and Cenozoic rocks, which reach 1–1.5 km in thickness at the point of impact (Masaitis, 1998; Pilkington et al., 2002). It is a multi-ring structure with three concentric rims visible.

The bedrock is crushed to depths of at least 5 km according to the results of drilling and geophysical measurements. The internal structure of the crater is quite unusual and contains a number of enigmatic phenomena, such as the presence of impact breccias fused into glassy, also fragmented tagamites (breccia within breccia). Vishnevsky and Montanari (1999) propose that the contrasting sedimentology or the presence of water in some layers of the original pre-impact sedimentary succession may have triggered a whole chain of impact phenomena. A similar result could be caused by the nearly simultaneous close impact of two or more meteorite fragments.

Long-term evolving terrains always have a complex tectonic framework, or rather a sequence of tectonic regimes creating a network of faults of different ages and directions. Masaitis (1998), in his diagram of the crater, shows radial tectonics in the immediate vicinity of the crater, while in Vishnevsky and Montanari (1999), long faults of NW–SE and SW–SE directions are displayed. Somewhat unexpectedly, faults in the N–S direction predominate in the crater itself, without any apparent influence of the impacting body.

Looking at the broader tectonic framework, we see a number of significant structures based on faults of approximately the N–S direction. The latter follows the Ural Mountains, the western margin of the Central Siberian Plateau, the Verkhoyansk Chrebet (belt), and some rivers such as the Daldyn River directly in the crater and around parts of the course of the Anabar and Malaya Kuonamka rivers. Perpendicular to them, a long E–W structure visible with ETOPO 1 (Fig. 1) is located north of the crater (see the arrows from E and from W).

The impact's shock pressure instantaneously transformed graphite in the ground into diamonds (e.g., Masaitis et al., 1972; Masaitis, 1998; Deutsch et al., 2000). The aggregates of diamonds are sometimes up to 1 cm. They tend to retain the appearance of graphite or original organic aggregates. They are bound to outcrops of original rocks with an admixture of graphite or coal substance. They are absent in the central part of the crater, where the pressure and temperature were too high for diamonds to form or be preserved (see, e.g., Fig. 1 in Masaitis, 1998). Vishnevsky and Montanari (1999) presented (see Fig. 6, p. 26) a diamond occurrence map showing a more or less chaotic distribution caused by both an irregular admixture of carbon-rich impacted rocks and a complex, multiphase crater evolution. Popigai is most probably linked to ejecta horizons occurring in marine sequences of late Eocene age.

Pilkington et al. (2002) presented Δg based on the local gravimetric data showing a negative “valley” going from the main and proven Popigai crater in the SE direction (see Fig. 3a in Pilkington et al., 2002), which is an indication of the possibility that we have double or multiple craters.

Popigai may be a multiple crater (Klokočník et al., 2010), a catena like on the Moon (Fig. 2a, b and slide nos. 6–21 in S3); it was proposed in Klokočník et al. (2010) based on analysis of Δg and T_zz derived from (at that time the best) a global gravity field model of the Earth, EGM 2008. Popigai may represent one of two or three simultaneous impacts from one original asteroid. Some authors have considered asteroid showers from a single parent-body breakup (Schmitz et al., 2015).

Popigai has been designated by UNESCO as a world geological heritage site. For economic reasons, exploration work in this remote area (the joint German–Canadian–Russian expedition) ceased before 2000 (according to information from Deutsch et al., 2000).

For completeness of these records, we note that near Popigai, roughly in the SW direction, there is another crater that is not yet proven, which is independent of and bigger than Popigai. This feature, known as Kotuykanskaya, is a hypothetical impact crater. It is located around φ = 69°30^′ N and λ = 100°25^′ E (Rajmon, 2009; Klokočník et al., 2020c, and references in this paper; also S3: 14–17).

5.2 Our new gravity results for Popigai

The recent global satellite-based surface topography depicted by the ETOPO 1 model is shown in Fig. 1 (and variants in S3). There is a broad topographic low in the NE, E, and SE directions from the main Popigai crater. One reckons that the terrain may have been strongly affected by water or ice erosion and other influences since the time of the impact event – e.g., a river is flowing throughout the bottom of the crater, and the rim is disrupted significantly in two places with consequences for the gravity signal (see the reaction in the gravity signal in the following figures).

Figure 1 Popigai: ETOPO 1 topography (m), shaded relief; the red dot indicates the crater's center. An alternative projection, with contour lines, is in S3. The arrows show an E–W-going linear structure; its gravity signal is weak.

The structure is characterized by a strong gravity low of $\mathrm{\Delta }g=-\mathrm{40}$ mGal and $T_{zz}=-\mathrm{30}$ E amplitude (EIGEN 6C4). Superimposed on the gravity low is a concentric ring-shaped high, which is presently fragmented, possibly due to post-impact evolution. The central peak is clearly visible in Figs. 2 and 3 but not with too much intensity.

Figure 2 Popigai: (a) Δg (mGal) and θ (°), I<0.9 with the topography; (b) T_zz (E), θ, plus the topography from ETOPO 1. Technical note: everywhere, cool colors are lows, warm colors are highs, and green is for zero. North: everywhere up (meaning in the direction to the present-day North Pole of the rotation of the Earth and present-day continent positions).

Figure 3 Popigai: details for θ in the main, largest, and proven Popigai crater. The halo of the strike angles is combed around the crater bottom (circle) and its central area. The Popigai river (in blue) locally disrupts the halo.

Figure 2a shows Δg, Fig. 2b presents T_zz, and Fig. 3 is a zoom just for θ (I<0.9) in the main crater, with a halo of the strike angles; there is a signature of the central peak, too. The topography (Fig. 1) and the gravity aspects (Figs. 2, 3, and S3) do not correlate.

Beside the main proven crater, we clearly see more candidates for impact craters (which are a bit smaller than the main crater). They are aligned in the SE direction (Klokočník et al., 2010, 2020b; Khazanovitch-Wulff et al., 2013). This is obvious from Fig. 2a and b, from broad negative Δg, from negative belts and semicircles of T_zz, and from the strike angles θ included in the figures (also S3: 8, 9, 21). These θ values have a tendency to be directed along the long axis of the whole Popigai family (SE–NW), interrupted only locally inside the potential craters (e.g., S3: 21). We labeled these crater candidates as Popigai II–IV in Klokočník et al. (2010). Counting from the main and proven Popigai crater (Popigai I), a large circular structure is visible on the SE rim of the main crater. It can be the companion crater – what we called Popigai II in Klokočník et al. (2010). At that time, we did not have the strike angles and the virtual deformations at our disposal. With them now, we can better demonstrate that Popigai can indeed be a double or multiple crater, i.e., a catena, a rare phenomenon on the Earth (the “Popigai family”).

6 Chicxulub

6.1 Introductory notes and geology

The impact crater Chicxulub (Northern Yucátan, México) is centered beneath Chicxulub village (φ = 21°17^′ N and λ = 89°30^′ W) near the Progreso port. The crater is huge and not exposed, with a diameter 170–250 km, and is about 65 million years old. This enormous impact represents an external forcing event (of non-terrestrial origin) with far-reaching, global consequences in mass extinction, as is well-known (the KT event).

The Yucátan peninsula is a low-lying limestone platform. The crater is buried under Quaternary carbonate sediments (0.6–1.0 km thick) lying over Tertiary sandstone and volcanic rocks. The northern (nearly) half of the now-buried crater is in shallow waters of the Sea of Campeche (of the Gulf of Mexico), which then falls, at the northern end of the Campeche Bank, to deep areas in the Campeche Escarpment (fault).

The origin of the impactor in the solar system is not yet clear. Bottke et al. (2007) proposed that the Chicxulub impactor could have originated from a moderately young asteroid family, Baptistina. Located in the inner main belt of asteroids, this cluster is favorably positioned to deliver large objects (>5 km) to the terrestrial planets. A recent analysis of Nesvorný et al. (2021) claims that the crater was produced by the impact of a carbonaceous chondrite and suggests that the impactor came from a main belt asteroid that quite likely (≃60 % probability) originated beyond 2.5 au. Some authors have discussed a comet as the impactor (e.g., Desch et al., 2021; now the minority opinion); the comet would come from the Oort cloud. The impactor might also be a binary asteroid, but it is rare (as we know) for two asteroids to produce two craters. The asteroids must be sufficiently separated (s/c “wide binaries”). Two closer impactors can produce one crater, one elongated crater, or two overlapping craters (Miljikovic et al., 2013).

The impactor's direction has been studied, among others, by Hildebrand et al. (2003), who noted that “the impact direction was towards the northeast based on the asymmetries preserved in various of Chicxulub's structural elements in addition to the vergence observed in the central uplift: compressional structures outside the crater rim, the rim uplift, compressional deformation preserved in the slumped blocks, morphology of the peak ring, off centre position of the central uplift in the collapsed disruption cavity (CDC), elongated CDC, and initiation of slumping of Cretaceous stratigraphy off the Yucátan platform”.

The literature about the Chicxulub crater is really rich and includes Alvarez et al. (1979, 1980), Smit and Hertogen (1980), Hildebrand (1998), Hildebrand et al. (2003), Ramos (1975), Campos-Enríquez et al. (2004), Gulick et al. (2008), Gulick and the Expedition 364 scientists IOPD (2016), Goderis et al. (2021), and Urrutia-Fucugauchi et al. (2022). We recall the important role of terrestrial gravity data in its study. Hildebrandt (1998) used not only terrestrial gravity anomalies (measured for oil and gas prospects) but also horizontal second derivatives to enhance resolution, but they were not aware of the concept of gravity aspects. According to Klokočník et al. (2010), Chicxulub may be a double crater; it was suggested after the analysis of Δg and T_zz based on EGM 2008 (compare with Fig. 2 in Hildebrand et al., 2003).

Strong impacts like this one have global effects; regionally, enormous pressure can trigger many post-impact activities and features. Let us recall Donofrio (1998), who wrote that “seventeen confirmed impact structures occur in petroliferous area of North America, nine of which are being exploited for commercial hydrocarbons … Disrupted rocks in proximity to impact structures, such as Chicxulub in the Gulf of Mexico off Yucátan, also contain hydrocarbon deposits”. James et al. (2002), p. 40, wrote that “… there are several craters that host fossil fuels, with the submarine Chicxulub impact crater …” and “… a total of 21 craters have oil/gas/hydrocarbon/coal resources, of which 19 host oil and gas”. The reader can see slide nos. 17–18 in S4.

A rapid burial of Chicxulub by Cenozoic sediments contributes to its preservation but also limits its study. The direct, surficial, or submarine geological study of Chicxulub is impossible because the structure is buried by several hundred meters to 1 km of porous Tertiary limestones (Ramos, 1975). A 2016 drilling project revealed a central ring composed of originally deep-seated, coarse-grained granite (Morgan et al., 2016). It is important because analogously we can expect rocks from depths of >10 km in, for example, lunar craters, as Kring (2016) reports for the crater Schrödinger. The concentric structure of Chicxulub is surrounded by a ring of cenotes. It indicates finely fractured and more permeable zones in the extensive cave systems developed. At the surface it manifests as cenotes, i.e., collapsed cave ceilings (Perry et al., 1995). In the wider surroundings, karst phenomena are known on the northeastern margin of the Yucátan in the Holbox tectonic zone, but here they are much more likely connected to the broad active arc that encircles Cuba from the north and trends toward the Yucátan (the Pinar Zone and Oriente Fault Zone).

6.2 Our new gravity results for Chicxulub

The gravity anomalies around Chicxulub are shown in Figs. 5–7. The radial component T_zz is in S4 (21). T_xx, T_yy, and T_zz are in S4 (22). The invariants I₁ and I₂ are in S4 (23), and their ratio I is in S4 (24), with VDs in Fig. 8 and S4 (24, 26, 27). The strike angles θ are in S4 (25); they are also underlying several other figures with the gravity aspects. We do not forget the ring of cenotes (sinkholes, originally potable water sources used by the Maya; S4: 9–11, 27, and 31).

The surface topography from ETOPO 1 (Fig. 4) does not correlate with the gravity aspects (the crater is not visible on the surface); even Ticul Fault and Ticul Sierra (hills) do not correlate with gravity.

Figure 4 Northern Yucátan, México, showing the flat lowland and shallow-sea area of the buried Chicxulub impact crater with ETOPO 1 topography (without any gravity aspect added) in an enhanced shaded relief scale (compare to S4: 8–10, 26). Tertiary sediments cover the impact crater; only semicircular “shadows” NE of the Ticul fault, due to the cenote rings, are visible in the plain terrain (here and in S4: 8). Black line: the coast.

We can see the positive T_zz at the central peak and along the rims and negative T_zz in between the rings. The strike angles are combed inside the crater and clearly laid down along the rims (analogy to Vredefort, S2: 9), so they also correlate with the ring of cenotes (S4: 10). Outside the central crater, the prevailing direction is SW to SE. The strike angles, combed around Chicxulub to halos and following the craters' rims, are conspicuous on land. The ring of cenotes agrees well with the halo created by the strike angles along the outer, most compact ring of the crater. Cenotes then continue like a cluster on the east edge of the crater (S2: 9). It is good to note that with the gravity aspects we can register the effect of all cenotes together, but not the individual features of course.

Tertiary sedimentary layers of the flat northern Yucátan outside the crater have, as expected, linear and also highly combed θ. A contrast of the density of sediment or a changed porosity (with respect to surrounding rocks) is high enough to be gravitationally distinguishable. The cenotes as well as oil and gas deposits near Yucátan, although epigenetic, are not there by chance (e.g., Grieve, 2005, p. 21) but as a consequence, direct or indirect, of the impact event.

Our figures demonstrate a halo around the central part (minimum of two rings). The strike angles are also strongly linearly combed far from the crater, mainly SW to NE (due to the local high porosity around and the cenotes outside the rims of Chicxulub E).

Figure 8 presents the virtual deformations (VDs), with red for dilatation and blue for compression. The VDs perfectly depict the bottom of the crater, its central peak, the rings, and the combed areas around.

We newly analyzed the negative “southern gravity anomaly” (located S to SW of the main crater) in the N–S direction; we call this feature the “tail” (see Figs. 5–7, S4: 15–17, 25, and 26). The prevailing standard opinion is that this is a pre-impact feature (e.g., Gulick et al., 2008; Urrutia-Fucugauchi et al., 2022).

Figure 5 Northern Yucátan, México, showing the area of the Chicxulub impact crater using the gravity model EIGEN 6C4 to maximum degree and order ($\mathrm{d}/\mathrm{o})\mathrm{2190}$, with a 4 km resolution. The gravity anomalies Δg (mGal), are shown together with the gravity strike angles θ (°), I<0.9. Black lines: coast and state borders. The strike angles as a parameter of the gravity anisotropy tensor Γ reveal up to three ringed structures of the Chicxulub basin. The combed strike angles correlate with oil and gas deposits (it continues to the SW to Campeche offshore oil fields), also with rims and a (semi)ring of the cenotes (on land). These are sinkholes (karst features) in the local limestone sediments; they were used by the Maya as a source of drinkable water. They represent one of the post-impact effects. The second radial derivative T_zz and other gravity aspects (including the combed strike angles with the comb statistics) are shown in S4 (25).

Figure 6 Gravity anomalies Δg (complete EIGEN 6C4 to $\mathrm{d}/\mathrm{o}=\mathrm{2190}$) (mGal) and strike angles (°) with ETOPO 1 topography. Circles are for Chicxulub I (proven) and II (hypothetical) impact craters; the ellipse is for the N–S-elongated depression of graben type (the “tail”), a part of the impact event.

Figure 7 Gravity anomalies (full EIGEN 6C4) (mGal) as contour lines and the ETOPO 1 topography as shaded relief. For more figures see S4: 5, 11, 12, 19, 20, and 28. Negative Δg values are in blue.

The tail, trench-like structure, or N–S-elongated depression of the graben type has a negative gravity anomaly. Linearly combed strike angles in the same direction (Figs. 5–6 and S4) indicate a syngenetic feature with the impact crater(s). The VDs in Fig. 8 show the best the whole linear feature, the hypothetical impact graben, connected with the impact craters. The southern tail would be its southern end.

Figure 8 The virtual deformations (VDs) (–) (compression in blue, dilatation in red) with EIGEN 6C4 on a 4 km grid. Black line: the coast. Yellow lines: hypothetical impact graben including a “tail” roughly in the N–S direction.

This tail is replicated in the younger relief uplifted SW of the Ticul fault (see ETOPO 1 morphology, Fig. 4). Extending the trench axis southward (Fig. 5), another linear depression (dark and light green) is encountered in a nearly perpendicular direction, trending northward and forming a “V”-like shape. For both these structures, we suggest that the influence of the impact on pre-existing geological structures may have been at work.

7 Discussion

7.1 Popigai

1. Beside the main proven crater, we clearly see more candidates for the impact craters; they are lined in the NW–SE direction (as we observed in Klokočník et al., 2010, and denoted as Popigai II, III, and IV). Here we confirm these our previous findings (Figs. 2, 3, and S3). The area SE of the main crater has negative values of Δg and T_zz and aligned strike angles θ.

2. Topography (ETOPO 1) and the gravity aspects do not correlate well. This indicates a partial smoothing of the impact features by erosion and filling of the impact-made depressions, in this case of both craters (including the hypothetical Popigai II crater) and the hypothetical crater in the NW- to SE-running impact trench.

3. The strike angles are combed into a halo around the main proven crater, Popigai I, and partly overlap the aligned but fragmented strike angles for Popigai II (Fig. 2a, b and S3).

4. A long, E–W-oriented structure, visible with ETOPO 1 (Fig. 1) and with the gravity aspects (Fig. 2a and b), is located north of the main crater, Popigai I (see the arrows in Fig. 1). We have mentioned in Sect. 5.1. that unexpectedly the N–S faults dominate in the crater without any apparent influence of the impacting body (Masaitis, 2003), while other tectonic schemes (Mashchak and Naumov, 2005) found evidence of the expected radial tectonics. Masaitis (1998), in his diagram of the crater, shows radial tectonics in the immediate vicinity of the crater. In Vishnevsky and Montanari (1999), however, long faults of NW–SE and SW–SE directions are displayed.

   The NW–SE linear structure connecting the craters is of particular interest to us because in the SE from Popigai I we can observe a long and wide depression to the distance of ∼400 km (Fig. 2a and b). This type of image is repeatedly encountered in most geological interpretations of the gravity data, typically, e.g., for ancient Nile valleys or lake basins covered by Saharan aeolian sands or hidden under Antarctic glaciers. We therefore assume that a depression filled with younger sediments extends south of the Popigai craters. According to analogies with other terrestrial structures, the thickness of the fill could be 1 km or more.

5. Given the close spatial association of the circular impact structure (the crater) with the linear NW–SE running “basin”, we guess the linear feature could be an original tectonic belt that was reactivated in extensional mode after the impact and subsequently filled with sediments in a dynamically evolving Cenozoic landscape. It could have been formed or influenced by Neogene movements related to the Tethys belt, but also by the periglacial regime of the Siberian north. Long-term evolving terrains always have a complex tectonic framework, or rather a sequence of tectonic regimes creating a network of faults of different ages and orientations.

Looking at the broader tectonic framework, one can see a number of significant structures based on faults of roughly N–S or NW to SE direction. They follow the Ural Mountains, the western margin of the Central Siberian Plateau, the Verkhoyansk Chrebet, and other rivers such as the Daldyn River directly in the crater and around. The NW–SE linear depression resembles an impact graben, i.e., a “trench modified by impact”. The heart of an impact graben is a pre-impact geological structure, activated by the impact energy to form a graben. This is not a new concept, as we observe basaltic rock eruptions in the extensional pressure regime in impacts on the Moon and Mars (e.g., Wichman, 1993; Spudis, 1993; Dasgupta et al., 2019; Zhang et al., 2023) and in some large terrestrial craters (Sudbury). In contrast, in the compressional tectonic regime, impact horsts are formed, such as those observed on the uplifted crater rims. The two stress–release or compression–extension regimes are complementary, usually perpendicular or oblique to each other. Especially in inhomogeneous terrestrial conditions (except perhaps in stable Archean blocks), meteorites strike areas with already existing regional stress fields. The stresses are then activated in specific directions by the enormous kinetic energy of the impactor.

Mashchak and Naumov (2005) stated the following: “… Thus, the 35-Ma-long post-impact modification history of the Popigai crater is determined by the superimposition of the regional tectonics on the long-term relaxation movements. As a whole, the late modification stage tectonics is found to have only an insignificant effect to the Popigai crater, so that both the original structure and the crater topography have been retained in a good state”. The existence of the circular structure of Popigai II and closely associated trench evokes the possibility of the almost concurrent formation of an impact crater (or craters) and an impact NW–SE-oriented graben.

7.2 Chicxulub

1. Topography (ETOPO 1) and the gravity aspects (namely Δg, T_zz, VD, and the invariants) do not correlate.

2. The majority of cenotes agree with the innermost ring (or the second ring, when counting the central ring around the central peak as the first ring), having positive Δg, T_zz, and strike angles combed into a halo.

3. The “southern tail” with negative Δg and T_zz and with the strike angles θ, aligned in the S–N direction, seems to be an inseparable part of one impact event (this impact may consist of several explosions). The strike angles, continuing from the main crater from its halo to the south, have a stream flowing from the halo to the S–N tail, slowly changing its direction from NW–SE to N–S; it looks like one common feature (the crater and the tail together).

4. Besides the main proven crater, we predicted (in Klokočník et al., 2010) another smaller crater in the NE direction. Accounting for all new gravity aspects, this still remains possible (see the circle in Fig. 6). Moreover, after a careful inspection, one can distinguish several more small circular features (in Fig. 5) near the Chicxulub crater (namely SW of it), which might also be impact craters, scattered around the primary. But this is just speculation.

Christeson et al. (2009) and others have argued that the gravity signal near Chicxulub is associated with the pre-existing Cretaceous basin proposed for this location (Gulick et al., 2008) rather than with an additional crater or craters. Our tools (Δg and T_zz) and EGM 2008 (predecessor of EIGEN 6C4) in 2010 were not sufficient to solve the problem. Moreover, we always wish to rely upon additional geological, geophysical, and other data, when available. In the meantime, with the gravity aspects, our tools improved and our experience with the gravity aspects increased. It is specifically the strike angle θ that proved to be very inspiring for diverse geo-applications in the case that stresses are present. The combed strike angles around Chicxulub create a halo (which is an expected and usual phenomenon for impact craters and basins, similarly as on the Moon or Mars), from which on its south side, a flow of θ changes direction to the south. There is no interruption, no jump, and no separation as we should observe between two separate geological features, telling us that the crater itself and its southern tail belong to one and the same body.

Previous studies have suggested asymmetries in the Chicxulub crater (e.g., Hildebrand et al., 2003; Gulick et al., 2008). This might be used to estimate the direction of the impactor in the atmosphere. However, seismic data show significant variations in the composition of the target rocks around the impact site. It is unclear whether the angle of impact or target material heterogeneity is responsible for the asymmetry (e.g., Collins et al., 2008).

Similarly, as for the Popigai family, we are interested in the linear structures near Chicxulub, namely in a trench-like structure running NW–SE of the main Chicxulub crater (Fig. 8). It is replicated in younger relief uplifted SW of the Ticul fault (see ETOPO morphology; Fig. 4). Extending the trench axis southward in Fig. 5, another linear depression (dark and light green) is encountered in a nearly perpendicular direction, trending northward and forming a V-like shape. For both of these structures, we suggest that a reviving influence of the impact on the pre-existing geological structures may have been at work.

Similarly to the Popigai crater family in Figs. 1–3, we can see in Figs. 4–8 how the circular impact structure is followed by a tectonic trench. In both craters, its direction roughly corresponds to the orientation of the surrounding geological structures. Thus, we assume that the faults, fault zones, or generally weakened structures already existed in these places before the impact. According to the gravity aspects, where the crater and the adjacent trench have a similar signal, we believe that the impact activated these structures form what we call an impact graben. However, both craters were rapidly filled with younger sediments, thus burying both the circular impact structure and the linear trenches.

Another interesting view is offered by Fig. 8, which shows the virtual deformation. Let us focus on the broad, blue lines that emanate from both arms of the crater to the NW. At the easternmost line, we observe a continuation along the shelf towards the edge of the continental slope, giving the impression of a valley formed in some impact-weakened zone. The western (marked in blue) zone is much longer, partially overlapping with the structures shown earlier (see Figs. 5–7).

For both crater formations, we consider the existence of and a relationship between their circular (crater) and linear components (graben-like structures). However, there is a different post-impact geological evolution for the linear trench-like structures, as they naturally become erosional pathways, and as such, they are subject to both down-cutting into the bedrock and filling with younger sediments in different ways at both locations (Siberia, Yucátan).

7.2.1 Astronomical note

The gravity aspects themselves cannot determine whether the impactor was a single body or a binary asteroid before its impact on the Earth. Both are possible. As noted above (Sect. 5.1), there is a possibility of breakup of one body (a single asteroid) in the atmosphere or a “flying cluster” of bodies encountering the atmosphere or wide binaries. To create a double crater, components of a binary asteroid must have a big distance or, in other words, a large separation (hundreds of kilometers). These are s/c wide binaries in contrast to close binaries. The velocity of asteroids in the solar system is much higher than the velocity of a point rotating on the Earth's surface. Thus, close binary asteroids can quickly hit one place twice and create only one crater.

As for the direction of the impactor, sometimes we can deduce this direction from the strike angles (Klokočník et al., 2020b), looking at how they are combed. As a good guide, we refer to Steinheim–Ries (S2: 18). Geologists know that the impactor(s) came roughly from the west, first creating the smaller Steinheim, then the bigger Ries. We can verify it independently using the strike angles; they are combed in the ∼ WE direction and skirt around both craters, creating a fragmented halo around Ries (disturbed by post-impact activities). For Popigai, we can expect an impactor coming from the S–E to NW, producing the small(er) crater(s) and finally the biggest, already proven one (Figs. 2–3, S3: 6–21).

7.2.2 Geological note

For the Popigai family, the gravity aspects decoded an arrangement of two or more possible craters in a line, which evokes a lunar catena. Smaller craters will have a weak and fragmented gravity record or they may have disappeared due to post-impact processes. For Chicxulub, our results suggest the existence of two craters (Fig. 6). Boreholes to the bottom of a shallow sea north of the NW Yucátan peninsula are known but not in the places where we would need them for the test of the second, smaller, unproven crater. More boreholes and seismic profiles at specific localities (Fig. 6) might clarify the situation (S4: 32).

8 Conclusion

We confirm and extend our results from Klokočník et al. (2010), which were based on analysis of Δg and T_zz derived from the gravity model EGM 2008. Now we work with the gravity aspects (including Δg and T_zz) and with the EIGEN 6C4 model. Thus, we have (in comparison with 2010) better tools (the set of gravity aspects) and a better model (EIGEN 6C4, with global gradiometric GOCE data). In turn, we are able to support or reject our older results with higher reliability and with more weight. The result is that we argue in favor of double or multiple craters – and bring further findings.

The impact affects or creates not only the circular structures but also other accompanying phenomena. These may be oriented concentrically like the cenote belts, but also as linear trenches, suggesting the existence of impact grabens. Their orientation and course depend on the regional tectonic architecture and stress fields prevailing at the time of the impact event and after it.

Popigai (Figs. 1–3) is probably a multiple crater, a catena (the smaller craters are located SE of the main proven crater). We consider at least Popigai II to be a proven crater by our new method and data. SE–NW is the probable direction of the impactor. The broad and long negative gravity anomaly in the SE direction of the main crater, Popigai I, indicates a close coupling between the circular impact structure and a linear depression. They are two possibilities: (1) the depression was formed by reactivation of older geological structures and is the impact graben. (2) The circular structure adjacent to the Popigai I crater in the SW gives the impression of another, perhaps shallower and more erosion-smoothed, impact crater, Popigai II. The gravity aspects at least partially suggest the possibility of a phenomenon that is uncommon on the Earth – the impact graben may actually represent a catena.

Chicxulub (Figs. 4–8) is probably a double crater; the smaller crater is located NE of the main proven crater. NE–SW is the probable direction of the impactor. The southern negative anomaly (the tail) belongs to the impact, as clearly demonstrated by the alignments of the strike angles and changes in their direction. The strike angles are combed into halos around the main crater (typical situation for all bigger impact craters), but then, on the southern side of Chicxulub I, they turn to the south (creating the tail). This tail could be the most southern end of the impact graben (Fig. 8) running NW, W to SW of Chicxulub I in the NW–SE direction.