The Chicxulub Impact Produced a Powerful Global Tsunami

Abstract

The Chicxulub crater is the site of an asteroid impact linked with the Cretaceous‐Paleogene (K‐Pg) mass extinction at ∼66 Ma. This asteroid struck in shallow water and caused a large tsunami. Here we present the first global simulation of the Chicxulub impact tsunami from initial contact of the projectile to global propagation. We use a hydrocode to model the displacement of water, sediment, and crust over the first 10 min, and a shallow‐water ocean model from that point onwards. The impact tsunami was up to 30,000 times more energetic than the 26 December 2004 Indian Ocean tsunami, one of the largest tsunamis in the modern record. Flow velocities exceeded 20 cm/s along shorelines worldwide, as well as in open‐ocean regions in the North Atlantic, equatorial South Atlantic, southern Pacific and the Central American Seaway, and therefore likely scoured the seafloor and disturbed sediments over 10,000 km from the impact origin. The distribution of erosion and hiatuses in the uppermost Cretaceous marine sediments are consistent with model results.

Introduction

The impact of an approximately 14-km diameter asteroid is implicated in the Cretaceous/Paleogene (K-Pg) mass extinction (Schulte et al., 2010) approximately 66 Ma ago. The bolide impact caused global temperature fluctuations (Schulte et al., 2010), large aerosol plumes (Bardeen et al., 2017), large plumes of soot and dust (Brugger et al., 2017), wildfires from ejecta re-entering the atmosphere (Busby et al., 2002; Morgan et al., 2013), and a massive tsunami (Kinsland et al., 2021; Matsui et al., 2002). Drilling cores from the Integrated Ocean Drilling Program (Gulick et al., 2016) and the International Continental Drilling Program (ICDP) corroborated the models (Collins et al., 2008) of the exact physical and geophysical nature of the crater and its peak ring which has facilitated detailed modeling of the impact (Morgan et al., 2016). Earlier tsunami simulations described the effects of the tsunami within the confines of the Gulf of Mexico (e.g., Matsui et al., 2002; Ward, 2012; see Ward [2021] for a more recent simulation extending beyond the Gulf of Mexico). Subsequent submarine landslides on the marine shelf (Gulick et al., 2008) could potentially increase the energy of this tsunami.

Most global tsunami simulations to date have been of tsunamis induced by underwater earthquakes, for instance, the 2004 Indian Ocean tsunami (Smith et al., 2005; Titov et al., 2005). Tsunami propagation has traditionally been simulated with shallow-water ocean models, which assume hydrostatic water pressure and a small depth-to-wavelength ratio. Such models cannot be used to simulate the complex first 10 min of the Chicxulub impact tsunami when there was large-scale deformation of the crust and the formation of a crater (Morgan et al., 2016). The crater formation and post-impact ejecta splashing back into the ocean create highly non-linear and non-hydrostatic waves. Modeling the impact tsunami requires a multi-stage simulation, with hydrocode modeling of crater formation and post-impact non-hydrostatic water waves, before hand-off of the solution to global shallow-water models. We pursue such a two-stage strategy in this paper. We discuss drawbacks of the models used, and the potential for improvements in future work in Section 5.4.

These linked models seek to depict a complex set of events associated with the asteroid impact and to predict the pathways of propagation as applied to a world with very different sea levels, ocean gateways, and continental positions and boundaries. The models do not incorporate a description of the chaotic near-field tectonic disturbances (e.g., faulting and slope failures) and the generation of smaller tsunamis by these disturbances. Did these aspects of the impact event alter the strength or the propagation pathway of the impact tsunami, or was this tsunami so powerful that these other effects were masked and overpowered? To verify the modeled strength and pathways taken by the impact tsunami we look at a global array of K-Pg boundary intervals in marine sections on land and in ocean drilling cores. In these sites, we will look for documented evidence of erosion, sediment disturbance, and/or redeposition of sediments that can be reasonably associated with the impact tsunami.

Impact Modeling Methods

We use the axisymmetric iSALE-2D hydrocode (Collins et al., 2004; Wünnemann et al., 2006) to simulate the initial formation of the Chicxulub impact tsunami. iSALE-2D has been used to simulate impact-induced tsunamis (e.g., Weiss & Wünnemann, 2007; Weiss et al., 2006; Wünnemann et al., 2010). The results of our iSALE-2D simulations were used to create initial conditions for shallow-water models to trace the tsunami throughout the world's oceans.

Motivated by impact simulations that reproduce the seismically imaged structure of Chicxulub (Collins et al., 2008) as well as the peak shock pressures and composition of the basin's peak-ring, as constrained by recent drilling (Morgan et al., 2016), we assume that the 14-km-diameter impactor had a density of 2,650 kg/m³ and struck Chicxulub at 12 km/s. Although the Chicxulub impact is thought to be oblique (45–60° from horizontal; Collins et al., 2020; Robertson et al., 2021) the axisymmetric nature of the code limits us to simulation of vertical impacts. We expect this limitation to have a minor effect on our results as the formation of the outward propagating rim wave is dominated by emplacement of slow ejecta that tends to be symmetric (e.g., Anderson et al., 2003). Our simulations have the same setup as those in Collins et al. (2008), but with a finer grid spacing and a larger domain needed to track the formation and early evolution of the tsunami (see Table S1 in Supporting Information S1 and other material in Supporting Information). We model the target as a granitic crust overlain by a 4-km-thick layer of sediments and an ocean with a constant depth of 1, 2, or 3 km (a 2-km ocean depth was used by Collins et al. (2008) for the northwestern sector of Chicxulub). With a grid resolution of 100 m, the ocean depth is resolved by 10, 20, and 30 cells, respectively, depending on assumed ocean depths of 1, 2, and 3 km. This number of grid cells is sufficient to resolve the rim wave (Bahlburg et al., 2010; Supporting Information S1). The atmosphere is not expected to significantly affect the early propagation of the tsunami. Thus, we do not include the atmosphere in our simulations. Further details of the iSALE simulations used in this paper, and their sensitivities to grid spacing, can be found in Supporting Information S1.

Results

The dimensions and formation of the crater are similar to previous work (Collins et al., 2008; Morgan et al., 2016). The results of our “fiducial” hydrocode impact simulation, with an assumed seafloor depth of 1 km and a run time of 10 min, are shown in Figure 1. About 2.5 min after contact of the projectile, a curtain of ejecta pushing water outward from the impact produced a 4.5-km-high wave (Figure 1a). After 5 min, falling ejecta continued to impart momentum to the ocean (Figure 1b). At 10 min, after all the ejecta had fallen, a 1.5-km-high wave, known as a rim wave, located 220 km from the point of impact was left propagating throughout the deep ocean (Figure 1c). Note that the majority of the wave breaking in the iSALE simulation takes place before the 10 min endpoint.

Figure 1. Formation of Chicxulub crater and the associated tsunami. Time series with material colored according to material type (crustal material is brown, sediments are yellow, and the ocean is blue). The origin marks the point of impact. Black curves mark material interfaces (e.g., sediment-crust interface). An animation of these results, from 0 to 10 min in steps of 5 s, is shown in Supporting Information Movie S1.

The axisymmetric nature of our high-resolution hydrocode model requires an ocean layer with a constant water depth. The ocean at the point of impact is estimated to be 100–200 m deep (Gulick et al., 2008), and it becomes deeper toward the northwest. Generation of the tsunami rim wave, however, is sensitive to the ocean depth at the crater rim, not at the point of impact. Paleobathymetry estimates indicate that water depth was ∼1 km where ejecta emplacement produces the initial rim-wave (50 km from basin center). At ∼150 km from the point of impact the ocean was ∼3 km deep (Figure S1 in Supporting Information S1). To test for sensitivity of the rim wave and crater shape to pre-impact ocean depth we vary the thickness of the ocean layer from 1 to 3 km. The waveforms after the first 10 min of the fiducial simulation, and after the first 10 min of iSALE simulations with different water depths, are displayed in Figure 2. These waveforms are in good agreement with the waveforms found in Bahlburg et al. (2010). Because of wave breaking and other processes that are not handled well in the shallow water models, our shallow water results will be sensitive to the transfer (“hand-off”) from the hydrocode to the shallow water model. To test for the sensitivity of the hand-off between the hydrocode and ocean model, we run a hydrocode simulation, with a larger mesh (see Supporting Information for more detail), out to 850 seconds before emplacement of the hydrocode conditions in the MOM6 model. Figure S4 in Supporting Information S1 demonstrates that handoff to the MOM6 “larger mesh” results at 600 and 850 s give nearly identical globally integrated energies. Surprisingly, the crater and rim wave structure at these early times do not depend strongly on assumed ocean depth within the range of 1–3 km (Figure 2). We do not expect this moderate dependence to hold over much deeper or shallower ocean depths. Our two-dimensional axisymmetric model with a constant depth is clearly a simplification of the bathymetry in the Gulf of Mexico. In the case of the 1 km ocean depth simulation, a sediment rim on the impact crater 10 min into the run rose above the water column, creating a ring-shaped island. The loose sediment in the rim would likely have been quickly dispersed by wave action (Bell et al., 2004). Other authors however have argued that resurge of water into the crater occurred by penetration through the raised rim and erosion allowing flow along the rim (Bahlburg et al., 2010). To test for sensitivity to this uncertainty, we model one initial condition with a sediment rim and one without. We test for sensitivity between the two runs and found the tsunami energies to be comparable (not shown). Therefore, the 1 km water depth iSALE simulation, with no sediment rim, is used for all subsequent runs.

Figure 2. Waveform and crater shape for three different runs from the iSALE-2D hydrocode. In the left-most part of the plot, the crater depths are shown. The middle and right-parts of the plot follow the change in sea level relative to the resting sea level. The crater depths are displaced by about 1 km from each other because of the differing ocean depths of the three runs.

Tsunami Propagation Modeling Methods

To simulate the global propagation of the impact tsunami, we use two different well-established shallow-water models: the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model Version 6 (Adcroft, 2017; MOM6), and Methods Of Splitting Tsunamis (Titov et al., 2016; MOST). The rim-wave has a wavelength of about 50–100 km, similar to the wavelengths of the 2004 Indian Ocean tsunami. As this is much greater than average ocean depths of about 4 km, the shallow water assumption, which assumes hydrostatic balance and is based on a comparison of wavelengths versus water depth, is well satisfied. The similarity of simulations from two different models using the same underlying shallow-water approximation and run on the same 1/10th degree grid but differing in their respective numerical implementations (more below) ensures robustness of our results. Neither of the models used here explicitly include dispersive effects. Discussion of potential effects of dispersion is provided in the section on Future Work.

Shallow-water models solve for perturbations to the resting sea surface elevations and for depth-averaged flow velocities. Flow velocities are the velocities of particles in the water, in contrast to the phase velocities of the tsunami wave propagating throughout the ocean. Errors due to this approximation are likely less than errors due to uncertainties in bathymetry. The large amplitudes of impact-generated waves lead to nonlinear dynamics during propagation, which is described only approximately by the shallow-water wave theory. Nevertheless, the long wave approximations have been successfully applied for simulating the nonlinear tsunami dynamics of propagation in shallow coastal regions and runup. Synolakis et al. (2008), for example, include an extensive discussion of verification and validation of shallow-water tsunami models with respect to field benchmarks. Their study demonstrates that large-amplitude waves can be predicted accurately with the shallow-water wave theory, providing the long-wave assumption is valid. The tsunami model benchmark efforts included a wide range of depth-integrated models (Pedersen, 2008) and initiated ongoing discussion about the proper use of the shallow-water and the Boussinesq-type models for tsunami simulations (Kirby, 2016). We address the dispersive modeling issues in the “Future Work” section.

MOM6 has been used to model tsunamis in the deep ocean, although it has not been used to forecast tsunamis. The barotropic solver in MOM6 is based on the solver in the Hallberg Isopycnal Model (HIM)/Generalized Ocean Layered Model (GOLD), which were used in the tsunami studies of Smith et al. (2005) and Kunkel et al. (2006). The results in Adcroft (2013) suggest that deep-ocean, large-scale motions are not overly sensitive to the horizontal resolution of the model. The forecasting accuracy of the tsunami calculation is not relevant for the application of the Chicxulub impact tsunami, but at 1/10th degree global resolution the arrival times are accurate to about 1%.

MOST was developed specifically for tsunami simulations (Titov & Synolakis, 1995; Titov et al., 2016). MOST has been extensively tested for various tsunami modeling applications and has been used to simulate historical tsunamis of different origins, including modeling of global tsunami propagation and local tsunami inundation impacts. MOST is now used operationally for tsunami forecasts at NOAA Tsunami Warning Centers. While MOM6 is run for all of the cases shown in this paper, MOST is run only in the fiducial case described below.

Both tsunami propagation models used the same global 1/10th degree bathymetric grid (Tables S2 and S3 in Supporting Information S1). To accurately simulate tsunami propagation, a global Maastrichtian (66 Ma) paleobathymetry is combined with the initial condition from the hydrocode results. The sources for the paleobathymetry are Müller et al. (2008) and Scotese (1997). More information about the bathymetries that we combined can be found in Text S1 of Supporting Information S1.

To continue the simulation with the tsunami propagation codes we convert the axisymmetric, constant water depth hydrocode results (see Figure S2 in Supporting Information S1) to more realistic, non-axisymmetric conditions with horizontally varying resting water depths. The hydrocode results at 600 s post-impact were used for the shallow-water model initial condition. At this time there was no more resolved falling ejecta; less voluminous and potentially fine ejecta would continue to fall after 600 s, but we do not expect that this more distal ejecta would significantly affect the rim wave. At 600 s, the waveform of perturbation sea surface heights is in approximate hydrostatic balance because the wavelength is much greater than the water depth (see Figure 2). The waveform, crater shape and velocity are isolated from the density profile. Assuming radial symmetry, the waveform is converted into a ring-shaped outward propagating wave, dependent on resting sea level, and inserted into the paleobathymetry described above. In the bathymetry the crater extended onto land where water was not present before impact. We test having a crater purely in water, without the portion of the crater that is formed over land (“Half Crater”), as well as a more complete crater that extended a full 360° onto land, (“Full Crater”), and compare energies, as discussed further below. The fiducial model employs the “Half Crater” bathymetry. More information on the blending of the hydrocode results into the paleobathymetry is given in Supporting Information S1.

To test sensitivity to the horizontal grid spacing of the shallow-water model, we run a shallow-water simulation at 1/5° grid spacing and compare snapshots of two-dimensional sea surface height perturbation fields (Figure S3 in Supporting Information S1) and energies (Figure S4 in Supporting Information S1) between this run and the nominal 1/10° run.

Results

Both shallow-water propagation models are run using the same fiducial run initial conditions and bathymetry data. Snapshots of the MOM6 and MOST sea-surface amplitudes are compared at the same times to ensure consistency of the results. The models display similar tsunami propagation patterns (Figure 3). The main dissimilarities in the model behaviors are in the later-stage wave dynamics. The differences reflect different numerical implementation of the shallow-water wave equations used in the two models. MOST is using the Godonov-type method (a Riemann solver) with a directional splitting, which emphasizes wave characteristics, and a discretization of non-linear terms in Lagrangian form. MOM6 employs vector invariant equations using an energy conserving discretization, with an emphasis on a well-behaved spectra in a turbulent cascade (not resolved or relevant to this problem). In addition, the bottom dissipation is parameterized differently in the two models. MOM6 displays more short-wavelength features after the initial, highest amplitude wave passing. Additional differences arise from different treatments of the north and south boundaries by MOM6 (reflecting boundaries) and MOST (absorbing boundaries without reflection). These model differences do not affect the leading order wave dynamics. The impact tsunami spread outside the Gulf into the Atlantic after about 1 hr from impact (Figure 3a); after 4 hr, through the Central American seaway, the waves enter into the Pacific (Figure 3b); after 24 hr of propagation, the waves cross most of the Pacific from the east and Atlantic from the west and entered the Indian Ocean from both sides (Figure 3c). The tsunami front propagates in excess of 200 m/s in deep water, in accordance with the shallow-water celerity. By 48 hr post-handoff, in other words, 48 hr after the handoff from the hydrocode to the shallow-water model, significant tsunami amplitudes have reached most of the world coastlines creating a complex amplitude pattern due to wave reflection and refraction (Figure 3d). Due to wave shoaling the tsunami heights can amplify many-fold near coastlines. The tsunami heights in open waters of the Gulf of Mexico are generally higher than 100 m. Along many North Atlantic coastal regions and some South America Pacific coastal regions the models show over 10 m offshore amplitudes. The simulations predict that around the world near-shore amplitudes exceed 1 m, with the exception of some coasts along the Indian Ocean and Mediterranean. Any historically documented tsunamis pale in comparison with such global impact. Depending on the geometries of the coast and the advancing waves, most coastal regions would be inundated and eroded to some extent. The simulations used here do not include wave runup onto land, as the model resolution of 1/10° is too low to resolve details of the inundation dynamics.

Figure 3. Comparison of two tsunami propagation models: MOST model—left column, MOM6—right column. Sea surface height perturbation in meters shown at (a) 1 hr and (b) 4 hr after impact around Gulf of Mexico, (c) 24 hr and (d) 48 hr post-handoff globally. Animations for both models are provided in Supporting Information Movies S3 and S4.

The maximum wave amplitudes and flow velocities (current speeds) at each model grid cell, over the 2-day time period of the MOST simulation, are respectively shown in Figures 4a and 4b. The largest waves and current speeds are in the Gulf of Mexico, North Atlantic, and South Pacific. Near the point of impact, the flow velocity exceeds 100 m/s. In other basins, flow velocities are up to a factor of 100 times smaller in the middle of the ocean than they are near the impact origin and along the coasts. Flow velocities above 20 cm/s are expected to cause erosion of fine-grained pelagic sediments (Lonsdale & Southard, 1974; McCave, 1984). Velocities higher than 20 cm/s are predicted in offshore areas of the North Atlantic and the equatorial region of the South Atlantic, in the Central American seaway and in most of the southern and southwestern Pacific, more than 12,000 km from the impact area. Most coastal areas of the world experienced above-20-cm/s velocities. As discussed in Supporting Information S1, tsunami propagation and flow velocities of simulations with slightly different input configurations (varying model resolution, friction coefficients, hand-off times between hydrocode and shallow-water models, crater size, etc.) are also tested for sensitivity. The energy of the tsunami is not greatly changed in any of these sensitivity tests except for the case in which the rim wave is removed.

Figure 4. (a) Maximum tsunami sea surface perturbation heights and (b) maximum flow velocity at each grid cell. Contours are shown for every meter of amplitude (saturated at 1,000 cm) and every 20 cm/s of speed. Contours of modern continents are shown for reference as gray lines. The results of the MOST model are shown here, because our MOST simulation saved output more frequently than our MOM6 simulations.

Geologic Verification of the Models Methods

Identifying the K-Pg boundary in marine sections requires some form of stratigraphic information. This is usually provided by biostratigraphic or paleomagnetic investigations. Marine sections located above present-day sea level and exposed on land usually allow a broad view of boundaries in outcrop and extensive stratigraphic data can often be collected from the section. Based on these studies and the overall preservation and exposure of the interval, one section has been named the “type section.” The stratotype section for the K-Pg boundary is at El Kef, Tunisia (XXVIIIth International Geological Congress, 1989). The boundary itself has been linked to the anomalous abundance of Iridium that was derived from the impacting body.

Close to the impact site reworked sedimentary deposits are mixed with ejecta from the impact. At intermediate distances the airborne ejecta may have arrived before the tsunami; thus, airborne ejecta with higher Iridium concentrations may lie below rip up clasts and redeposited older sediments. In distant regions, high concentrations of Iridium used to define the K-Pg boundary (Kiessling & Claeys, 2001) are thought to have settled in high concentrations over a period up to several years (Claeys et al., 2002; Toon et al., 1982). This is compared to the modeled tsunami reaching a global extent in just 2 days. To verify the strength and pathway of the modeled impact tsunami we pay particular attention to these more distal regions (Schulte et al., 2010). In these regions the effects of the tsunami should be found in the interval immediately below the K-Pg boundary itself in both marine sections on land (Supporting Information Table S5) and in scientific ocean drilling cores (Supporting Information Table S6). We take any sign of missing biostratigraphic or paleomagnetic intervals or depositional disturbance immediately below this level (e.g., erosional truncations of bedding or bioturbation features, sediment deformation, allochthonous clumps or clasts) as evidence of current activity or disturbance associated with the impact tsunami.

A few of the studied boundary sections have paleomagnetic stratigraphy. The K-Pg boundary has been found to be within the upper half of subchron C29r in Gubbio, Italy (Lowrie & Alvarez, 1977) and Agost, Spain (Canudo et al., 1991). The estimated duration of the Cretaceous part of subchron C29r is 300 kyr (Husson et al., 2011). Biostratigrapy often provides an important indication of missing section in deeper water, pelagic sections. The Abathomphalusmayaroensis Zone defines the uppermost Cretaceous foraminiferal zone in many of the older studies of the K-Pg boundary; however, Keller (1988) found that A. mayaroensis disappeared below the K-Pg boundary in the type section at El Kef. To fill this gap, Pardo et al. (1996) defined a total range biozone (Plummerita hantkeninoides) that marks the top of the Maastrichtian and lies within the lower half of subchron C29r. The uppermost Cretaceous nannofossil zone is defined by the range of Micula prinsii (Sissingh, 1977). This Zone occupies most of the lower half of subchron C29r. Other fossil assemblages have been used to evaluate the ages within the Late Cretaceous, but they have not been well documented in more than one or two complete K-Pg boundary sections. Carbon and oxygen isotope stratigraphies have been generated for several of the K-Pg boundary sections (Supporting Information Tables S5 and S6). The records of the carbon isotopes show an abrupt break at the K-Pg boundary, with the isotopes becoming sharply negative (a drop of 2‰ at El Kef; Keller & Lindinger, 1989). However, the oxygen isotopes signal is more variable and may depend on what microfossils are being measured (cf., El Kef, in Keller and Lindinger [1989], MacLeod et al. [2018], and other sections in Caravaca, in Canudo et al. [1991]; and in Zumaia, in Margolis et al. [1987]).

The advent of orbital tuning of geologic records has greatly advanced our ability to develop estimates of ages with comparable precision well back into the Cretaceous (e.g., Batenburg et al., 2012; Dinarés-Turrel et al., 2014; Husson et al., 2011 and references therein). These studies use calculated variations in the Earth's orbit as a template for matching variations in stable isotopes, color, iron content, or bed thickness; however, beyond 60 Ma only the 405 kyr eccentricity cycle is known with sufficient accuracy to be used in tuning the time scale (Laskar et al., 2011; Westerhold et al., 2012). From these tuning efforts we know that the K-Pg boundary lies at the top of the 405 kyr orbital cycle of eccentricity designated as Ma405 1 (Batenburg et al., 2012). Any effort at tuning must take place within a stratigraphic framework defined by other tools, normally a paleomagnetic stratigraphy, which in turn usually relies on a biostratigraphic framework.

For the drill sites reported in this study, we list those sites (Supporting Information Table S6) in which the K-Pg boundary interval is recovered and is fossiliferous, with stratigraphies that define the location of that boundary. Based on stratigraphic evaluations for both drilled cores and outcrop sections, we class the K-Pg boundary sections as: (1) complete, (2) apparently complete, (3) having a detectable depositional disturbance, hiatus, or disconformity, or (4) having a long erosional hiatus or non-depositional surface (Figure 5). If such long missing sections range from the Cretaceous well up into the Paleocene or even younger sections, we cannot claim that they are attributable to the impact tsunami (category 4, above), and we discount them from our analysis. If, however, the lower part of the Paleocene is present, while a part of the Upper Cretaceous is missing, we classify this as possibly caused by the impact tsunami (category 3, above).

Figure 5. Plate reconstruction and site locations at the age of the K/Pg boundary, from ODSN website (http://www.odsn.de/odsn/services/paleomap/paleomap.html) using the magnetic reference frame. Continental blocks in gray with modern continental outlines in red. Green shaded ocean areas depict approximate regions where the models of the K-Pg impact tsunami showed flow velocities in excess of 20 cm/s (see Figure 4b). Most coastal regions were indicated by the models to have experienced such high velocities, but are not shown here. Drill site locations indicated by circles; K-Pg land outcrop sites indicated by squares (see legend). Small filled circles indicate sites with hiatuses of a million years or more duration that span the K-Pg boundary and range well into the Paleogene.

Results

The devastating effects of the asteroid impact in the Caribbean and Gulf of Mexico included earthquakes, slope failures, and debris flows, all of which could have contributed to tsunami formation (e.g., Alegret & Thomas, 2005; Alvarez et al., 1992, 1995; Bourgeois et al., 1988; Bralower et al., 1998; Campbell et al., 2008; Denne et al., 2013; Keller et al., 1997, 2007; Kinsland et al., 2021; Maurrasse & Sen, 1991; Montanari et al., 1994; Sanford et al., 2016; Schulte et al., 2006, 2008; Smit et al., 1996; Stinnesbeck et al., 1997). These ancillary effects are not accounted for in the impact tsunami models, but nevertheless disrupted the K-Pg boundary. The modeled impact tsunami took principal radiation pathways directed to the east and northeast into the North Atlantic and to the southwest, through the Central American passage and into the southwestern Pacific (Figure 4). At flow speeds greater than 20 cm/s (Figure 4b) the passing tsunami could have eroded fine-grained marine sediment even on the deep seafloor (Lonsdale & Southard, 1974; McCave, 1984).

The Tethys region, the South Atlantic, the North Pacific, and the Indian Ocean basins were largely shielded from the stronger effects of the tsunami (Figure 4). This is consistent with the location of the several complete sections described from the marine outcrops around the Mediterranean, including the type section at El Kef (Figure 5). It is also consistent with the frequent recovery of complete sections at scientific ocean drilling sites in the South Atlantic Ocean and on Seymour Island in the Antarctic Peninsula, the several complete sections of the K-Pg boundary recovered in the North Pacific Ocean and on the island of Hokkaido, and the complete K-Pg boundary intervals drilled on bathymetric highs in the eastern Indian Ocean.

Looking at K-Pg boundary intervals that lay in the modeled pathway of the tsunami, the results of the comparison are also largely consistent. The drilled sections in New Jersey show gaps, rip up clasts, or tempestites in the K-Pg boundary interval. Sections studied in western Europe (Germany, Denmark, France, Bulgaria, Austria; Supporting Information Table S5) generally show biostratigraphic gaps, erosional truncations, or slumps and gravity flows in the uppermost part of the Maastrichtian section. In the North Atlantic Ocean only three sites in two areas contain what appear to be complete K-Pg boundary intervals (Figure 5). Site U1403 is the deepest site drilled on the J-Anomaly Ridge off Newfoundland. The Upper Cretaceous section is relatively thick here, lying between two southeast trending basement highs (Expedition 342 Scientists, 2012) and may represent a depocenter for sediment eroded from the nearby locations. Sites 1259 and 1260 are located on the slope of the Demerara Rise off Suriname, South America. During the Late Cretaceous their location was within a few degrees north of the equator and may have been partially shielded from the main force of the tsunami (MacLeod et al., 2007). However, farther south on the coast near Recife, Brazil, at Pernambuco, a neritic section contains a graded sandy bed, including ejecta from the asteroid impact, and is overlain by an iridium anomaly (Albertão & Martins, 1996).

Almost all the drill sites in the South Pacific basin appear to have a missing uppermost Maastrichtian section. This is true even on the southern part of the Ongtong-Java Plateau, which lies near the northern edge of higher velocities associated with the impact tsunami's modeled pathway, while two sites on the northern side of the Plateau (Sites 803 and 807) have the only complete K-Pg sections recovered in the South Pacific basin (Figures 4b and 5).

Of particular interest are the outcrops of the K-Pg boundary interval on the southeast corner of North Island and northeast corner of South Island, New Zealand. Here the olistostromal deposits at the top of the Upper Cretaceous Whangi Formation were originally explained as the result of local tectonic activity (Laird et al., 2003) or mass flow deposit (Hines et al., 2013); but considering the stratigraphic position of this deposit and its location directly in line with the modeled pathway of the impact tsunami, we feel the olistostrome is recording the effects of the impact tsunami (Figures 4–6). Hollis (2003) reviewed 16 marine sections in New Zealand that ranged in paleo water depth from inner shelf to upper bathyal and found that at least 14 of them probably had a missing or disturbed K-Pg boundary interval. However, detailed biostratigraphic control of the uppermost Maastrichtian is lacking for the remaining two sections, which raises the possibility that these sections may also be incomplete. Paleomagnetic control on the sections has not been obtained due to pervasive demagnetization (Kodama et al., 2007).

Figure 6. The percent of apparently complete marine sections containing the K-Pg boundary interval listed by ocean basin, including both marine sections found on the surrounding land and in scientific ocean drilling cores. Data do not include sites with long hiatuses (see text). The number of sections studied is shown at the base of each column (see Supporting Information Tables S5 and S6). No complete sections were found in the Caribbean (including the Gulf of Mexico). The South Atlantic and South Pacific categories include sites studied in the Southern Ocean sector of these basins.

The tsunami models indicate that many coastal regions around the globe may have been affected by the impact tsunami. However, without a detailed knowledge of the bathymetry and coastal geometry at the end of the Cretaceous, and without a higher resolution model in these areas, we cannot evaluate how accurate the models might be in such shoreline areas. Our study shows that some distant near-shore areas were strongly affected (e.g., New Jersey, New Zealand, Pernambuco), while others were not (e.g., Seymour Island, Hokkaido). Still, it is probably significant that the models show only minor coastal effects in the shielded Tethys basin (Figures 5 and 6) where all the neritic sections appear to be complete (Supporting Information Table S5).

In a similar manner, all the large, relatively shallow oceanic plateaus and rises show up in the higher velocity regions of the models (Figure 4b); however, as in the coastal regions, the resolution of the models and that of the paleo bathymetry do not allow detailed comparison of the model results with the completeness of the recovered sections. We feel it is significant that only those prominent bathymetric highs that lie outside the main pathway of the impact tsunami show a preponderance of complete K-Pg sections (Figures 4b and 5).

A summary of the studied marine sections on land and in drill cores is shown in Figure 6. As noted above, all marine sections on land around the Mediterranean lie outside the modeled >20 cm/s flow velocity contour (Figures 4b and 5) and are believed to have complete K-Pg boundary records. Also noted above, the Caribbean-Gulf of Mexico region lies within the area of very high flow velocities and have no complete, undisturbed sections. Similarly, the North Atlantic Basin is an area of high flow velocities and has only four sites (11% of sites studied) that are apparently complete (Supporting Information Tables S5 and S6). The South Pacific region with flow velocities >20 cm/s (Figures 4b and 5) have two sites (11% of sites studied) that appear to be complete.

At least 65% of the studied sections in regions where modeled flow velocities are <20 cm/s have complete sections. In regions with flow velocities >20 cm/s, 91% of the studied sections have incomplete K-Pg boundary sections. The most telling confirmation of the global significance of the impact tsunami is the highly disturbed and incomplete sections on the eastern shores of North and South Islands of New Zealand. These sites lie directly in the path of the tsunami propagation, more than 12,000 km distant from the impact location (Figures 4 and 5).

Discussion Tsunami Mechanisms

Earlier theoretical and regional simulations (e.g., Matsui et al., 2002; Ward, 2012; Wünnemann & Weiss, 2015) differ on whether the rim wave or collapse wave dominates with respect to energy. The rim wave refers to the water displaced from the impact that is pushed away from the origin (Wünnemann & Weiss, 2015). The collapse wave is the secondary process arising from the cavity collapse in the crater and water rushing into the crater (Wünnemann & Weiss, 2015). To test the relative contributions of the collapse and rim waves to the total tsunami energy, we ran a simulation (“Crater Only”) with no rim wave or velocity, such that the tsunami is solely due to the collapse wave filling in the crater. Our results agree with the conclusion of Wünnemann and Weiss (2015), that the rim wave is the source of most of the energy for this impact tsunami. Four hours after impact, the “Crater Only” case is about 13 times less energetic than the “Full Crater, With Rim Wave” case. The MOM and MOST model simulations of the Half Crater scenario showed similar energy numbers 4 hr post-handoff (3.90 × 10¹⁹ and 3.84 × 10¹⁹ J correspondingly), such that the model energy estimates appear to be robust and independent of the exact model used.

The efficiency of tsunamis can be quantified by the ratio between tsunami energy and the source energy. The efficiency of tsunami generation by the Chicxulub impact is similar to that of large earthquakes. The energy ratio for earthquake-generating tsunamis averages around 0.1% (with large variations from 0.02% to 0.8%, Tang et al., 2012), while we predict that the Chicxulub tsunami has an efficiency of 0.19% (Table S4 in Supporting Information S1). Figure S4 in Supporting Information S1 shows that the impact tsunami energy dissipates relatively quickly, relative to seismogenic tsunamis, consistent with the “Van Dorn effect” (Van Dorn et al., 1968) of faster wave energy attenuation due to large non-linearities near the source of explosion-generated tsunamis. Near-field tectonic activity, triggered by passage of a strong stress wave produced by the impact, was not included in our simulations. It is likely that any effects of earthquake generated slides and collapses would be minor relative to the primary rim wave.

Hiatus Distribution

The better preserved, thicker, carbonate-rich sections in the oceans are commonly found on bathymetric highs such as continental terraces, oceanic plateaus, rises, aseismic ridges, and seamounts. Drill sites in which the K-Pg boundary is clearly identified are usually found in such locations. These locations do have their own problems, however. Such regions of bathymetric prominence also give rise to enhanced turbulence in the waters surrounding them (Cacchione & Drake, 1986; Cacchione et al., 2002; Rudnick et al., 2003; Wunsch & Ferrari, 2004); thus, they enhance the erosional power of tsunamis and tidal waves that pass over them. The preserved sedimentary sections atop bathymetric highs usually show clear evidence of erosion and the sculpting of pelagic deposits that sit upon them. The drilling strategy often employed by scientific ocean drilling expeditions takes advantage of the stratigraphic character of these deposits to sample relatively older intervals where overburden has been removed or was never deposited, the intention being to minimize the effects of diagenetic alteration on these older sediments. At other sites, full advantage was taken of the thicker, more complete sections to study the detailed paleoceanographic history. This duality of purpose means that many sites drilled on bathymetric highs contain significant gaps in the stratigraphic record, while on some highs there are close-by sites that have recovered complete sections. In regions with modeled flow velocities <20 cm/s, several sites locate the K-Pg boundary between recovered cores (Supporting Information Table S6); thus, the amount of missing section (if any) and the exact nature of the boundary is uncertain.

In basins where almost all sites show incomplete uppermost Maastrichtian sections there are still a few deep-sea sections that appear to be complete (e.g., Sites 1259, 1260, U1403 in the North Atlantic). These may represent local bathymetric shielding from erosion or local depocenters that receive sediment which has been eroded from nearby areas. The coincidence of regions having few if any complete K-Pg boundary sections and the pathway of relatively strong tsunami flow, combined with the more common occurrence of complete K-Pg boundary sections in regions that did not have strong tsunami flow, support the results of the tsunami models. The lack of complete K-Pg boundary sections in the southern South Pacific and on the eastern shores of New Zealand strongly suggest that this tsunami was of global significance, reaching at least 12,000 km across the deep ocean. It also suggests that except for some shallow coastal regions, areas such as the Tethyan region, the North Pacific, the South Atlantic and much of the Indian Ocean basin were largely geographically shielded from the effects of the tsunami.

Comparison With Large Historical Tsunamis

To provide perspective on the size of the impact tsunami, we compare our impact tsunami model estimates with some representative large historical tsunamis. The 2004 Indian Ocean tsunami (Smith et al., 2005) is possibly the largest modern-era tsunami; it killed over 230,000 people around the Indian Ocean and was recorded around the globe (Titov et al., 2005). The 2011 Tohoku tsunami was generated by a similarly strong earthquake and has become the costliest natural disaster of all time. Offshore amplitudes of the 2004 Indian Ocean tsunami 2 hr after generation were measured to be about 0.6 m, and 2 m waves were measured about 500 km away from the epicenter of the 2011 Tohoku tsunami, at a seafloor depth of 5,700 m. These deep-ocean amplitudes led to runup at coastlines of up to 40 m (for the Indian Ocean tsunami at Sumatra Island) and 50 m (for the Tohoku tsunami at Honshu Island). The 1883 Krakatau event generated another catastrophic tsunami with explosive-type initial conditions, potentially similar to the impact generation. The Krakatau wave devastated local coastlines, killing over 30,000 (second most deadly record after the Indian Ocean tsunami) with waves that ran up to 40 m and traveled distances of up to 5 km inland, but did not generate significant waves outside Sunda Strait. All these tsunamis, among the largest in recorded history, are dwarfed by the wave amplitudes and energy of the simulated Chicxulub tsunami. The Chicxulub tsunami produces offshore amplitudes over 1 m around most of the world oceans (Figure 4a). When tsunamis reach the shallow waters of a coastline or bathymetric high, wave amplitude increases due to shoaling. Comparison of our tsunami simulations with observations and modeling of the strongest recent tsunamis of 2004 and 2011 implies that the coastal amplitudes for the Chicxulub tsunami would flood most coastlines, in a manner that would be catastrophic in modern times. The total energy of our impact tsunami simulations is compared with the energy of these large historical tsunamis in Table S4 and Figure S4 of Supporting Information S1. Energy values are calculated according to standard formula for shallow-water energy (e.g., Arbic et al., 2004; their equation 14). Figure S4 in Supporting Information S1 displays the ratios of energy in the impact tsunami simulations to the 2004 Indian Ocean tsunami, as a function of time into the ocean simulation. The energy in the impact tsunami decays faster than the energy in the 2004 Indian Ocean tsunami—another manifestation of the “Van Dorn effect.” The initial energy in the impact tsunami was up to 30,000 times larger than the energy of any historically documented tsunamis. Wave energies in the “Half Crater” simulation are about 5% less than those in the “Full Crater” simulation. The “Crater Only” simulation, without the large rim wave, still has much more energy than any other historical tsunamis. For a wide variety of sources, the portion of the source energy that goes into tsunami generation is less than 1%, with large variations from about 0.01% to 0.3% (Table S4 in Supporting Information S1). An impact- and explosion-type of tsunami generation appears to have similar efficiency in transferring energy into long wave propagation. However, impact- and explosion-generated tsunamis dissipate energy much faster during propagation. Nevertheless, the sheer amount of energy of the impactor is sufficient to generate a giant global tsunami, even if only 0.2% of the impact energy goes into the tsunami.

Future Work

The first global simulation of the Chicxulub impact tsunami demonstrates that it was much larger than any recent earthquake-generated tsunami, and that it was likely large enough to leave a mark on marine sediment records. Many uncertainties remain, and there is much room for improvement in future studies. It is well known that most impacts are oblique with 45° impact angle being most likely (e.g., Robertson et al., 2021). With sufficient computer power, high-resolution, three-dimensional hydrocode simulations of the first 10 or so minutes could be performed, thus allowing for varying water depth, non-perpendicular impact angles, and other key uncertainties in the hydrocode simulation. Generally, we would expect a slightly larger rim wave in the downrange direction and a smaller wave up range. It may be instructive to vary initial conditions of the global simulation in a parameterized way to crudely account for impact angle.

For the present paper, to make the hydrocode simulations feasible, we have employed two-dimensional hydrocode simulations. There are tradeoffs in the usage of two-dimensional hydrocodes versus shallow-water codes, such that there will never be a perfect time to perform the handoff between them. For instance, the hydrocode handles breaking waves, whereas our shallow-water codes do not. On the other hand, the two-dimensional hydrocode simulations assume a constant resting water column depth whereas the shallow-water codes allow for more realistic horizontally varying depth. In the present paper, we tested two different handoff times (850 s, vs. 600 s). Future work could explore sensitivity to handoff times in more detail.

The 15 January 2022, Hunga Tonga-Hunga Ha'apai volcano explosion has demonstrated an additional mechanism of tsunami generation from large explosive events—the low frequency air pressure wave, also known as a Lamb wave (Duncombe, 2022). While the exact tsunami-generation mechanism of the air pressure Lamb wave is not fully understood, it is clear that significant waves can be generated from such air pressure waves propagating over oceans. The full analysis of such tsunami generation is out of the scope of this paper and is a subject of future research. But based upon observations and initial modeling of the Tonga event, it is clear that the Lamb wave can be a source of significant secondary tsunamis around the world. These waves would reach world coastlines much earlier than the tsunami generated by the crater formation. The energy of the Chicxulub impact is at least 100,000 times larger than the Tonga explosion. The Lamb wave from the Tonga explosion generated tsunami waves of over a meter at some locations around the Pacific and up to half a meter at other oceans. Thus, the Lamb wave from the Chicxulub explosion can be a significant source of tsunamis in the far-field from the impact source, and will be a subject of future work.

Dispersive effects may manifest themselves in the Chicxulub tsunami propagation simulations in two ways: (a) during the long-distance propagation as different wave frequencies separate from a single front; and (b) during the evolution of the initial steep wavefront into an undular bore (Glimsdal et al., 2007). Tsunami amplitudes in shallow water wave approximation models may overpredict shorter dispersive waves or underpredict sharp frontal amplitudes experiencing fission and undular bore formation. In both cases the difference may be up to 50% of amplitudes in certain cases (see e.g., Son et al., 2011; Zhou et al., 2012, 2014). Addressing these effects is a topic for future research. Both of these processes generally lead to the decrease of amplitudes in comparison with the classic shallow-water wave theory estimates. Therefore, the non-linear shallow water approximation provides, in general, a conservative (upper-bound) estimate of potential tsunami amplitudes. The use of Boussinesq-type models may provide a better resolution of the undular bore feature of the turbulent wavefront. However, these effects involve generation of much shorter (therefore much more dissipative) wavelengths that are usually confined to a relatively small part of the wave near the bore front (see e.g., Matsuyama et al., 2007; Son et al., 2011), and therefore may have very limited effect on the global wave propagation pattern—the main goal of this study. Also, the results of Glimsdal et al. (2007) show that the Boussinesq model appears to overestimate the dispersive front effects in comparison with the full hydro code, which may be attributed to differences in resolution or to the inherent tendency of Boussinesq models to overestimate dispersion. The detailed modeling of the dispersive front of the leading tsunami with higher spatial resolution dispersive simulations would show more precise dynamics of the tsunami in the near-source area and may change the details of the maximum amplitude distribution near the source. Therefore, such studies with higher resolution dispersive models would be a natural extension of this work, especially for more precise estimates of tsunami impact within the Gulf of Mexico. However, we don't expect these details to significantly change our far-field estimates of the tsunami amplitudes and tsunami energy directionality (Zhou et al., 2012, 2014).

In the case of our modeling, we expect the dispersive effects would be, at least partially, accounted for, since one of the models (MOST) includes the physical process of frequency dispersion approximated by numerical dispersion (Burwell et al., 2007). MOST has been benchmarked against laboratory tests with highly dispersive and highly non-linear waves for wave breaking dynamics (Titov & Synolakis, 1995) and compared with dispersive models during the long-distance tsunami propagation (Zhou et al., 2012). These comparisons showed that MOST provides results closely resembling the dispersive model estimates. The consistency of MOST and MOM6 results provides confidence in the robustness of our results. However, dispersive effects as well as uncertainties such as in the details and size of the impactor, and in the paleo-bathymetry estimates should be investigated more fully in future work.

Acknowledgments

We thank Gareth Collins, Finn Løvholt, Clemens Rumpf, and two anonymous reviewers for comments and suggestions that led to significant improvements in the manuscript. B.K.A. thanks Sarah Stamps for useful conversations. M.M.R. and B.K.A. thank Mark Champe, Charles Antonelli, and Michael Messina for help with setting up MOM6 on University of Michigan computers. M.M.R., B.K.A., and J.K.A. acknowledge funding support from US National Science Foundation grants OCE-0968783 and OCE-1351837, a Research Experience for Undergraduates (REU) supplement for MMR to OCE-1351837, and the University of Michigan Associate Professor Support Fund supported by the Margaret and Herman Sokol Faculty Awards. We gratefully acknowledge the developers of iSALE-2D, including Gareth Collins, Kai Wünnemann, Dirk Elbeshausen, Tom Davison, Boris Ivanov and Jay Melosh. Some plots in this work were created with the pySALEPlot tool written by Tom Davison. The MOM6 simulations in this paper were carried out on the Flux supercomputer provided by the University of Michigan Advanced Research Computing Technical Services. Much of B.K.A.'s contributions to this paper took place while he was on sabbatical in France. B.K.A. thanks many French colleagues, especially Thierry Penduff, Rosemary Morrow, Nadia Ayoub, and Florent Lyard, for their help in procuring this sabbatical year. He Wang's contributions from GFDL are supported by NOAA's Science Collaboration Program and administered by UCAR's CPAESS under awards NA16NWS4620043 and NA18NWS4620043B. This is contribution number 5300 of NOAA-PMEL.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

Results, bathymetry, and initial conditions of our work can be found at https://doi.org/10.7910/DVN/GWOFIO.