- Email: [email protected]
comet on the Earth
ARTICLE IN PRESS
International Journal of Impact Engineering 33 (2006) 713–722 www.elsevier.com/locate/ijimpeng
Numerical simulations of hypervelocity impact of asteroid/comet on the Earth T. Saitoa,, K. Kaihob, A. Abec, M. Katayamac, K. Takayamad a
Department of Mechanical Systems Engineering, Muroran Institute of Technology, 27-1 Mizumoto-cho, Muroran 050-8585, Japan b Institute of Geology and Paleontology, Tohoku University, Aoba, Aramaki, Sendai 980-8578, Japan c CRC Solutions Corp., 2-7-5 Minamisuna, Tokyo 136-8581, Japan d Biomedical Engineering Research Organization, Tohoku University, 2-1-1 Katahira, Sendai 980-8577, Japan Available online 15 November 2006
Abstract In this study, numerical simulations of asteroid/comet impact on the earth were carried out by using the hydrocode AUTODYN-2D (Century Dynamics Inc.). The ejected mass into the atmosphere originated from different sources was recorded at different altitudes as a function of time. It is claimed that the mass extinction at the Permian–Triassic (P/T) boundary was caused by an extraterrestrial impact and the following abrupt climate/environment change. The objective of this work is to substantiate the claim with numerical analysis. It is demonstrated that the materials from different origins are splashed into the atmosphere and that the crater size is a function of the diameter as well as the kinetic energy of the impactor. The effect of impact angle to the appearance of ejected materials is also investigated using two-dimensional plane strain model. r 2006 Elsevier Ltd. All rights reserved. Keywords: Simulation; Impact; Hypervelocity; Asteroid; Comet
1. Introduction In the history of the Earth, occurrences of species mass extinction have been recorded with approximate cycles of 26.6 million years. Among them, two devastating ones took place 65 and 251 million years ago. The former is the so-called Cretaceous–Tertiary (K/T) boundary separating the mesozoic era from the cenozoic one and the latter is the Permian–Triassic (P/T) boundary in between the paleozoic era and the mesozoic era [1–3]. It is widely accepted that an extraterrestrial impact was the cause of the extensive mortality at the K/T boundary [4–6]. However, the cause of the most severe biotic crisis in the Earth’s history at the P/T boundary leading to the subsequent origination of the Holocene biota on the Earth has not been determined, yet [7]. Until recently, biological and environmental changes at the end of Permian were thought to be more gradual than those recorded at the K–T boundary. Accordingly, many earth scientists have considered Corresponding author. Tel.: +81 143 46 5354; fax: +81 143 46 5353.
E-mail address: [email protected] (T. Saito). 0734-743X/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijimpeng.2006.09.012
ARTICLE IN PRESS 714
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
different explanations such as gradual sea-level fall, climate change, oceanic anoxia, and volcanism, etc. for the end-Permian event [7]. Then some evidences indicating an abrupt end-Permian extinction of organisms in both marine and terrestrial environments were presented [2,3,7–10]. Specifically, Becker et al. [11] reported that some end-Permian sediment contains fullerenes with trapped noble gases that are indicative of an extraterrestrial source. Based on the geochemical and paleontological analyses of well-preserved sediment at southern China, Kaiho et al. [12] claim that the mass extinction at the P/T boundary was caused by an extraterrestrial impact and the following abrupt climate/environment change. They explain the severe climate change as a result of a gigantic release of sulfur from the Earth’s mantle due to the impact. As a result of chemical reactions of the sulfur with atmosphere, significant oxygen consumption and the production of strong acid rain catastrophically changed the global environment. Large-scale forest fires, shielding effect of the sunshine by dust particles floating in stratosphere also occurred. All these events attributed in the break down of food chain and the global extinction resulted. A 130-km diameter impact crater was discovered in May 2000 at Woodleigh in Australia and scientists are investigating the possibility of this crater being the evidence of the impact at the P/T boundary. Kaiho also indicated a hypothesis that the huge crater recently discovered in Siberia might be a spalling hole created with this event. The objective of this paper is to substantiate the claim presented by Kaiho et al. by showing the possibility that the Earth’s mantle is released into the atmospheric–oceanic system by an extraterrestrial impact. Numerical simulation is one of the most powerful means of investigating dynamics of such nonstationary phenomena and the impact event has been numerically investigated by many researchers [13–15]. In this study, numerical simulations of asteroid/comet impact on the Earth were carried out by using the hydrocode AUTODYN-2D (Century Dynamics Inc.). Accumulations of ejected mass from different sources into the atmosphere were recorded at different altitudes. 2. Numerical simulations 2.1. Hydrocode For the purpose of analyzing hypervelocity impact phenomena, numerical tools that are categorized as hydrocode are generally used. Wilkins originally formulated it as HEMP code in the 1950s [16]. The fundamental governing equations are based upon the classical continuum mechanics. They describe the dynamic behavior of a continuous medium with a set of partial differential equations established through the conservations of mass, linear momentum and energy. The equation of state and the constitutive equation describing material properties are solved simultaneously together with the initial and boundary conditions. Many hydrocodes such as the HEMP code have adopted finite difference schemes, while others have been formulated in the finite element method. Almost all of them, employ explicit time integration schemes in order to solve extremely time-dependent transient phenomena. 2.2. Numerical modeling AUTODYNs-2D has been used extensively for hypervelocity impact problems [17]. The numerical code has several different solvers such as Eulerian, Lagrangian, arbitrary Lagrangian Eulerian (ALE), smoothed particle hydrodynamics (SPH) and others. In this study, the Eulerian scheme was used. The code automatically handles problems with multiple materials. In order to set up simulations with a hydrocode, appropriate descriptions for the equation of state (EOS), the constitutive relations including the material strength and the failure model together with physical/ mechanical properties of each material must be specified. The Tillotson EOS is suitable for the current study since it is capable of taking the effect of shock-induced vaporization into account [18]. It becomes equivalent to the Mie–Gru¨neisen–type shock-Hugoniot EOS in the lower-pressure region (p1 TPa). The Tillotson EOS takes four different forms depending on the degree of material compression m defined as m ¼ (r/r01) and the internal energy E, where r and r0 are the current and the reference densities,
ARTICLE IN PRESS T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
715
respectively. – Region I (mX0): the pressure, p1, is calculated by bTil p1 ¼ aTil þ Zr0 E þ ATil m þ BTil m2 , 1 þ E=E 0 Z2
(1)
– Region II (mo0, EpEs): the pressure, p2, is calculated by Eq. (1) by putting BTil ¼ 0. – Region III (mo0, EsoEo E0 s): the pressure, p3 , is calculated by p3 ¼ p2 þ ðp4 p2 Þ
E Es , E 0s E s
– Region IV (mo0, E0 spE): the pressure, p4, is calculated by ( ) bTil Zr0 E 1 1 2 p4 ¼ aTil Zr0 E þ þ ATil m exp b 1 exp a 1 , 1 þ E=E 0 Z2 Z Z
(2)
(3)
where E is the current specific internal energy (energy per unit mass) and Es is the one at the sublimation point. The symbol E0 s is defined as E0 s ¼ Es+eEv with the vaporization energy Ev, and e is the fraction of the vaporization energy to be added to Es to guarantee gas-like behavior as r approaches zero. The symbol E0 is the specific internal energy at 0 1C. Note here that the values p2 and p4 in Eq. (2) are evaluated by substituting Es and E0 s into the formula for the Region II and Region IV, respectively. The variables indicated by aTil, bTil, ATil, BTil, a, b are the characteristic properties of each material. Region I represents the compressed phase of the material and extends to shock pressures of about 1000 Mbar. Region II describes material that is compressed to an energy less than the sublimation energy Es. Therefore, the energy will be released adiabatically and the pressure returns to zero. Presumably, the material keeps the state as solid. There is no provision in this EOS to describe the material at pressures less than zero. Region IV is the expansion phase of material which has been shocked to an energy E0 s sufficiently large to ensure that it will expand as a gas at very large expansions. Since details of constitutive models for the solid materials are not expected to have significant effect in the current study of extremely high-velocity impact, a simple elastic-perfectly-plastic model (VonMises model) was employed. The spall failure criterion was used for solid materials. The material instantaneously fails when the hydro static pressure reaches its negative limiting value.
3. Hypervelocity impact simulations 3.1. Axisymmetric simulations: normal collisions Two cases of the normal impact are numerically simulated. In the first case (Case 1), a 21 km diameter asteroid (specific gravity: 4) impacts on the Earth with the velocity of 20 km/s, while in the second case (Case 2), a 13 km diameter comet (specific gravity: 1) collides with the Earth with the velocity of 60 km/s. Fig. 1 shows the simulation domain. Since calculations were carried out with axisymmetic geometry, only half of the region is displayed. The space resolution is 1 km/grid point and 300 and 361 cells were used in the radial and axial directions. The simulation includes the 5 km deep ocean and the 6 km thick crust in between the atmosphere and the mantle. Water is used as the oceanic material. The oceanic sediment layer is neglected since the thickness is approximately 1 km and no significant influence is expected with the current grid resolution. Major elements of the crust are SiO2 (60%) and Al2O3 (15%) and specific gravity of 2.8–2.9 while for the mantle they are SiO2 (45%) and MgO (35%) with specific gravity of 3.3. It is assumed that cordierite (2MgO–2Al2O3–5SiO2) and steatite (MgO–SiO2) have similar properties with the crust and mantle,
ARTICLE IN PRESS 716
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
Fig. 1. Schematic of the computational domain for axisymmetric calculations.
respectively. Their mechanical properties used in the current study are taken from the matweb database (http://www.matweb.com). In this study, judging from the impact velocity, it is implicitly assumed that the large slower impactor (Case 1) is an asteroid and the smaller one with much faster velocity (Case 2 and Case A–D in the next section) is a comet. The material of the asteroid and comet are not known and they are assumed rather arbitrarily: the asteroid is pure limestone and the comet consists of 75% limestone and 25% ice. Their densities are 4000 and 3250 kg/m3, respectively. Numerical models and the physical parameters adopted in this study are summarized in Table 1. 3.2. Axisymmetric numerical results and discussions Fig. 2 shows the material distribution with velocity arrows at different times for Case 1. The times after impact are indicated in the figure. The shock wave propagating in the mantle is clearly seen in the figure at 10 s after the impact. The depth of the crater asymptotically reaches its maximum of about 70 km in about 30 s, while the diameter still increases 60 s after the impact. It is noted that the interface between the atmosphere and the mantle inside the crater becomes unstable at this late stage. Although the calculation stopped at this moment, it is expected that this is the beginning of the collapsing stage. From the figure, it is observed that the diameter of the crater is approximately 120 km. Fig. 3 shows the accumulation of different materials that are splashed into the atmosphere as a function of time. Between 10 and 15 km altitude, the material originated from the crust amounts the most and those from the seawater and the mantle follow. Between 20 and 30 km altitude, the material from the crust also amounts to the largest value and material from the seawater and the mantle is relatively less compared with the lower altitude. In both cases, practically no material from the asteroid is observed. Fig. 4 shows the material distribution with velocity arrows at different times for Case 2. Although the calculation stopped at 20 s after the impact, it looks as if the depth of the crater is almost at its maximum of about 60 km, which is less than Case 1. The final diameter is about 80 km and it is also smaller compared with that in Case 1. Considering the fact that the impact energy of Case 2 is about 1.7 times larger than that of Case 1, the result demonstrates that the size or the volume of a crater due to normal impact of an asteroid or a comet is a function of the impactor size as well as the impact kinetic energy. It is noted that the ejecta seems more violent compared with that in Case 1, and more melted impactor material recoils into the atmosphere.
ARTICLE IN PRESS T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
717
Table 1 Numerical models adopted and physical properties
Atmosphere
Material
Eq. of state
Constitutive Eq.
Failure model
Void
—
—
—
a
Ocean
Water
Polynomial Ref. density r0 ¼ 1000 Kg/m3
—
Hydro Tensile limit 3.0 MPa
Crust
Cordieriteb (Dolomite)
Tillotson Ref. density r0 ¼ 2800 Kg/m3
VonMises Shear Modulus: 53 GPa Yield Stress: 344 GPa
Hydro Tensile limit 27.5 MPa
Mantle
Steatiteb (Dolomite)
Tillotson Ref. density r0 ¼ 3300 Kg/m3
VonMises Shear Modulus: 45 GPa Yield Stress: 619 GPa
Hydro Tensile limit 68.8 MPa
Asteroid
Limestoneb (Limestone)
Tillotson Ref. density r0 ¼ 4000 Kg/m3
VonMises Shear Modulus: 53 GPa Yield Stress: 344 GPa
Hydro Tensile limit 27.5 MPa
Comet
Limestone+iceb (Limestone)
Tillotson Ref. density r0 ¼ 3250 Kg/m3
VonMises Shear Modulus: 53 GPa Yield Stress: 344 GPa
Hydro Tensile limit 27.5 MPa
a
Bakken and Anderson 1967 [19]. Tillotson EOS parameters for ‘dolomite’ and ‘limestone’ [20] except densities are used. The densities of asterod and comet are determined from astronomical considerations. b
Fig. 2. Axisymmetric numerical results of Case1: diameter; 21 km, impact velocity; 20 km/s.
ARTICLE IN PRESS 718
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
Fig. 3. Accumulation of ejected mass from different sources at different altitude as function of time for Case1: left; between altitudes of 10 and 15 km, right; between altitudes of 20 and 30 km.
Fig. 4. Axisymmetric numerical results of Case 2: diameter; 13 km, impact velocity; 60 km/s.
ARTICLE IN PRESS T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
719
Fig. 5. Accumulation of ejected mass from different sources at different altitude as function of time for Case 2: left; between altitudes of 10 and 15 km, right; altitudes of 20 and 30 km.
Fig. 5 shows the accumulation of different materials that are splashed into the atmosphere at the specified altitude. Although it is still very little, the material from the impactor is visible in Case 2 reflecting the more violent features of impact in the early stage. According to the geological observations, more material from the impactor (asteroid or comet) is expected. The main reason for the current result, in that a very little amount of material from the asteroid/comet is ejected into the atmosphere, is most likely due the assumption of normal impact. In general, the most probable impact angle is 451 and the situation can be quite different from that of normal impact. The effect of impact angle is investigated in the next section using two-dimensional (2D) simulations. Among the material parameters, the density is considered to be the most important one for hypervelocity impact. In the present work, actual densities of material are used. However, other Tillotson EOS parameters are taken from those either for Dolomite or Limestone as indicated in Table 1. Although we cannot quantitatively estimate how much error would be induced with this approximation, we can still hope that the numerical results give us useful information. Investigations of the Tillotson parameters that can control the shock-induced vaporization are needed in future in order to accurately estimate the amount of mass splashed into the atmosphere. 3.3. 2D simulations of oblique collisions Hypervelocity collisions of asteroid or comet with the Earth at different impact angles were carried out in order to investigate the effect of impact angle to the appearance of the ejecta and the melted impactor. Simulations were carried out with the 2D plane strain model due to restricted computational resources. The diameter of the impactor (comet) is 13 km and the impact velocity is set to 60 km/s. Note that the shape of the comet now is essentially cylindrical while it is spherical in the previous section. The impact angles measured from the vertical direction are 01, 301, 451 and 601 and are named as Case A, Case B, Case C and Case D, respectively. Fig. 6 shows the domain of numerical simulation. Other numerical parameters such as the grid resolution are the same as in the previous axisymmetric calculations. 3.4. 2D results and discussions Fig. 7 compares the results of axisymmetric calculation (Case 2) and 2D calculation for a 13 km diameter comet impacting normally to the earth with the impact velocity of 60 km (Case A). The figures displayed are at
ARTICLE IN PRESS 720
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
Fig. 6. Schematic of computational domain and asteroid impact angles: Atmosphere 50 km, Ocean 5 km, Crust 6 km, Mantle 300 km.
Fig. 7. Normal impact simulations with axisymmetric (left: Case 2) and 2D geometries (right: Case A) 4 s after impact.
4 s after the impact event. The kinetic energy is 6.3 1024 J for the axisymmetric case and 7.8 1020 J for the 2D case assuming the width of comet to be a unit length. Then, the kinetic impact energies per unit colliding area are obtained by dividing the kinetic energies with the projection areas of the comet as 5.1 1016 and 6.0 1016 J/m2 for the axisymmetric and 2D cases, respectively. The impact energy density is about 18% higher in the 2D case. Although the general appearances of the two are more or less the same, the 2D case ejects much more material into the atmosphere. The difference between the two is more than expected and it can be attributed to the different ways of how impact energy is diffused after the impact. It is also noted that the impact material (comet) is completely disintegrated in the axisymmetric case while it is only partially so in the 2D simulation. It is natural to consider the axisymmetric case to be more realistic and we have to keep this observation in mind.
ARTICLE IN PRESS T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
721
Fig. 8. Oblique impact simulations with 2D geometry: diameter; 13 km, impact velocity; 60 km/s, 4 s after impact, impact angles: Case A: 01, Case B: 301, Case C: 451, Case D: 601.
Fig. 9. Oblique impact simulations with 2D geometry: 13 km, impact velocity; 60 km/s, 8 s after impact, impact angles: CaseA: 01, Case B: 301, Case C: 451, Case D: 601.
Namely, although the simulations for Cases A–D may provide us with some general data of the effect of impact angle on the current impact problem, we still need three-dimensional simulations to get some quantitative conclusions. Fig. 8 compares 2D results (Cases A–D) at 4 s after the impact. As expected, the material is ejected almost into the direction of mirror reflection to the incident impact angle. It is also noticed that, for a larger impact
ARTICLE IN PRESS 722
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
angle, the depth of crater is smaller and its diameter is larger. Fig. 9 shows the results at 8 s after the impact. Comparing it with Fig. 8, it is noticed that the angle dependency of ejecta flow gets weaker in later stages. The difference between Cases C and D at 4 s after the impact (Fig. 8) is quite notable but Cases C and D in Fig. 9 are not much different with each other. It is known that the shape of crater is circular regardless of the impact angle. The tendency observed here in Figs. 8 and 9 may be reflection of this known fact. 4. Conclusions Two cases of normal impact of an asteroid/comet with the earth are numerically simulated. It is demonstrated that the materials from different origin are splashed into the atmosphere. The size, i.e. the depth and the diameter of the crater, is a function of the diameter as well as the kinetic energy of the impactor. It is demonstrated that the crater size of a hypervelocity impact with higher impact energy can be smaller if the diameter of the impactor is smaller. The effect of impact angle to the flow of ejected materials is investigated using 2D numerical analysis. Although it is only qualitative, a tendency that the dependency of the appearance of mass ejection into the atmosphere on the impact angle is more significant in the early stage and becomes less notable at later time is observed. References [1] Kaiho K, Lamolda MA. Catastrophic extinction of planktonic foraminifera at the Cretaceous–Tertiary boundary evidenced by stable isotopes and foraminiferal abundance at Caravaca, Spain. Geology 1999;27:355–8. [2] Rampino MR, Adler AC. Evidence for abrupt latest Permian mass extinction of foraminifera: results of tests for the Signor–Lipps effect. Geology 1998;26:415–8. [3] Bowring SA, Erwin DH, Jin YG, Martin MW, David EK, Wang W. U/Pb zircon geochronology and tempo of the end-Permian mass extinction. Science 1998;280:1039–45. [4] Kaiho K, Kajiwara Y, Tazaki K, Ueshima M, Takeda N, Kawahata H, et al. Oceanic primary productivity and dissolved oxygen levels at the Cretaceous/Tertiary boundary: their decrease, subsequent warming, and recovery. Paleoceanography 1999;14:511–24. [5] Alvarez W. Toward a theory of impact crises. EOS (Transactions, American Geophysical Union), 1986; 67: 649, 653–655, 658. [6] Sheehan PM, Fastovsky DE, Barreto C, Hoffmann RG. Dinosaur abundance was not declining in a ‘‘3 m gap’’ at the top of the Hell Creek Formation, Montana and North Dakota. Geology 2000;28:523–6. [7] Erwin DH. The great Paleozoic crisis. New York: Columbia University Press; 1993 327p. [8] Retallack GJ. Permian–Triassic extinction on land. Science 1995;267:77–80. [9] Eshet Y, Rampino MR, Visscher H. Fungal event and palynological record of ecological crisis and recovery across the Permian–Triassic boundary. Geology 1995;23:967–70. [10] Jin YG, Wang Y, Wang W, Shang QH, Cao CQ, Erwin DH. Pattern of marine mass extinction near the Permian–Triassic boundary in south China. Science 2000;289:432–6. [11] Becker L, et al. Impact event at the Permian–Triassic boundary: evidence from extraterrestrial noble gases in fullerenes. Science 2001;291:1530–3. [12] Kaiho K, Kajiwara Y, Miura Y. End-Permian catastrophe by a bolide impact: evidence of a gigantic release of sulfur from the mantle: reply. Geology 2002;September:856. [13] Pierazzo E, Melosh HJ. Hydrocode modelling of Chicxulub as an oblique impact event. Earth Planet Sci Lett 1999;165:163–76. [14] Jones AP, Price GD, Price NJ, DeCarli PS, Clegg RA. Impact induced melting and the development of large igneous provinces. Earth Planet Sci Lett 2002;202:551–61. [15] Stoeffler D, Artemieva NA, Pierazzo E. Modeling the Ries–Steinheim impact event and the formation of the Moldavite strewn field. Meteoritics Planet Sci 2002;37:1893–908. [16] Wilkins ML. Calculation of elastic-plastic flow. Lawrence Livermore Laboratory Report, UCRL-7322, revision 1, 1973. [17] Robertson NJ, Hayhurst CJ, Fairlie GE. Numerical simulation of explosion phenomena. Int J Comput Appl Technol 1994;7(3–6):316–29. [18] Tillotson JH. Metallic equations of state for hypervelocity impact. GA-3216, General Atomic, CA, July, 1962. [19] Bakken LH, Anderson PD. The complete equation of state handbook. SANDIA SCLTM-67-118, November, 1967. [20] Hageman LJ, Walsh JM. HELP: a multi-material Eulerian program for compressible fluid and elastic-plastic flows in two space dimensions and time. BRL Contract Report, vol. 1(39), 1971. p. 26.
International Journal of Impact Engineering 33 (2006) 713–722 www.elsevier.com/locate/ijimpeng
Numerical simulations of hypervelocity impact of asteroid/comet on the Earth T. Saitoa,, K. Kaihob, A. Abec, M. Katayamac, K. Takayamad a
Department of Mechanical Systems Engineering, Muroran Institute of Technology, 27-1 Mizumoto-cho, Muroran 050-8585, Japan b Institute of Geology and Paleontology, Tohoku University, Aoba, Aramaki, Sendai 980-8578, Japan c CRC Solutions Corp., 2-7-5 Minamisuna, Tokyo 136-8581, Japan d Biomedical Engineering Research Organization, Tohoku University, 2-1-1 Katahira, Sendai 980-8577, Japan Available online 15 November 2006
Abstract In this study, numerical simulations of asteroid/comet impact on the earth were carried out by using the hydrocode AUTODYN-2D (Century Dynamics Inc.). The ejected mass into the atmosphere originated from different sources was recorded at different altitudes as a function of time. It is claimed that the mass extinction at the Permian–Triassic (P/T) boundary was caused by an extraterrestrial impact and the following abrupt climate/environment change. The objective of this work is to substantiate the claim with numerical analysis. It is demonstrated that the materials from different origins are splashed into the atmosphere and that the crater size is a function of the diameter as well as the kinetic energy of the impactor. The effect of impact angle to the appearance of ejected materials is also investigated using two-dimensional plane strain model. r 2006 Elsevier Ltd. All rights reserved. Keywords: Simulation; Impact; Hypervelocity; Asteroid; Comet
1. Introduction In the history of the Earth, occurrences of species mass extinction have been recorded with approximate cycles of 26.6 million years. Among them, two devastating ones took place 65 and 251 million years ago. The former is the so-called Cretaceous–Tertiary (K/T) boundary separating the mesozoic era from the cenozoic one and the latter is the Permian–Triassic (P/T) boundary in between the paleozoic era and the mesozoic era [1–3]. It is widely accepted that an extraterrestrial impact was the cause of the extensive mortality at the K/T boundary [4–6]. However, the cause of the most severe biotic crisis in the Earth’s history at the P/T boundary leading to the subsequent origination of the Holocene biota on the Earth has not been determined, yet [7]. Until recently, biological and environmental changes at the end of Permian were thought to be more gradual than those recorded at the K–T boundary. Accordingly, many earth scientists have considered Corresponding author. Tel.: +81 143 46 5354; fax: +81 143 46 5353.
E-mail address: [email protected] (T. Saito). 0734-743X/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijimpeng.2006.09.012
ARTICLE IN PRESS 714
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
different explanations such as gradual sea-level fall, climate change, oceanic anoxia, and volcanism, etc. for the end-Permian event [7]. Then some evidences indicating an abrupt end-Permian extinction of organisms in both marine and terrestrial environments were presented [2,3,7–10]. Specifically, Becker et al. [11] reported that some end-Permian sediment contains fullerenes with trapped noble gases that are indicative of an extraterrestrial source. Based on the geochemical and paleontological analyses of well-preserved sediment at southern China, Kaiho et al. [12] claim that the mass extinction at the P/T boundary was caused by an extraterrestrial impact and the following abrupt climate/environment change. They explain the severe climate change as a result of a gigantic release of sulfur from the Earth’s mantle due to the impact. As a result of chemical reactions of the sulfur with atmosphere, significant oxygen consumption and the production of strong acid rain catastrophically changed the global environment. Large-scale forest fires, shielding effect of the sunshine by dust particles floating in stratosphere also occurred. All these events attributed in the break down of food chain and the global extinction resulted. A 130-km diameter impact crater was discovered in May 2000 at Woodleigh in Australia and scientists are investigating the possibility of this crater being the evidence of the impact at the P/T boundary. Kaiho also indicated a hypothesis that the huge crater recently discovered in Siberia might be a spalling hole created with this event. The objective of this paper is to substantiate the claim presented by Kaiho et al. by showing the possibility that the Earth’s mantle is released into the atmospheric–oceanic system by an extraterrestrial impact. Numerical simulation is one of the most powerful means of investigating dynamics of such nonstationary phenomena and the impact event has been numerically investigated by many researchers [13–15]. In this study, numerical simulations of asteroid/comet impact on the Earth were carried out by using the hydrocode AUTODYN-2D (Century Dynamics Inc.). Accumulations of ejected mass from different sources into the atmosphere were recorded at different altitudes. 2. Numerical simulations 2.1. Hydrocode For the purpose of analyzing hypervelocity impact phenomena, numerical tools that are categorized as hydrocode are generally used. Wilkins originally formulated it as HEMP code in the 1950s [16]. The fundamental governing equations are based upon the classical continuum mechanics. They describe the dynamic behavior of a continuous medium with a set of partial differential equations established through the conservations of mass, linear momentum and energy. The equation of state and the constitutive equation describing material properties are solved simultaneously together with the initial and boundary conditions. Many hydrocodes such as the HEMP code have adopted finite difference schemes, while others have been formulated in the finite element method. Almost all of them, employ explicit time integration schemes in order to solve extremely time-dependent transient phenomena. 2.2. Numerical modeling AUTODYNs-2D has been used extensively for hypervelocity impact problems [17]. The numerical code has several different solvers such as Eulerian, Lagrangian, arbitrary Lagrangian Eulerian (ALE), smoothed particle hydrodynamics (SPH) and others. In this study, the Eulerian scheme was used. The code automatically handles problems with multiple materials. In order to set up simulations with a hydrocode, appropriate descriptions for the equation of state (EOS), the constitutive relations including the material strength and the failure model together with physical/ mechanical properties of each material must be specified. The Tillotson EOS is suitable for the current study since it is capable of taking the effect of shock-induced vaporization into account [18]. It becomes equivalent to the Mie–Gru¨neisen–type shock-Hugoniot EOS in the lower-pressure region (p1 TPa). The Tillotson EOS takes four different forms depending on the degree of material compression m defined as m ¼ (r/r01) and the internal energy E, where r and r0 are the current and the reference densities,
ARTICLE IN PRESS T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
715
respectively. – Region I (mX0): the pressure, p1, is calculated by bTil p1 ¼ aTil þ Zr0 E þ ATil m þ BTil m2 , 1 þ E=E 0 Z2
(1)
– Region II (mo0, EpEs): the pressure, p2, is calculated by Eq. (1) by putting BTil ¼ 0. – Region III (mo0, EsoEo E0 s): the pressure, p3 , is calculated by p3 ¼ p2 þ ðp4 p2 Þ
E Es , E 0s E s
– Region IV (mo0, E0 spE): the pressure, p4, is calculated by ( ) bTil Zr0 E 1 1 2 p4 ¼ aTil Zr0 E þ þ ATil m exp b 1 exp a 1 , 1 þ E=E 0 Z2 Z Z
(2)
(3)
where E is the current specific internal energy (energy per unit mass) and Es is the one at the sublimation point. The symbol E0 s is defined as E0 s ¼ Es+eEv with the vaporization energy Ev, and e is the fraction of the vaporization energy to be added to Es to guarantee gas-like behavior as r approaches zero. The symbol E0 is the specific internal energy at 0 1C. Note here that the values p2 and p4 in Eq. (2) are evaluated by substituting Es and E0 s into the formula for the Region II and Region IV, respectively. The variables indicated by aTil, bTil, ATil, BTil, a, b are the characteristic properties of each material. Region I represents the compressed phase of the material and extends to shock pressures of about 1000 Mbar. Region II describes material that is compressed to an energy less than the sublimation energy Es. Therefore, the energy will be released adiabatically and the pressure returns to zero. Presumably, the material keeps the state as solid. There is no provision in this EOS to describe the material at pressures less than zero. Region IV is the expansion phase of material which has been shocked to an energy E0 s sufficiently large to ensure that it will expand as a gas at very large expansions. Since details of constitutive models for the solid materials are not expected to have significant effect in the current study of extremely high-velocity impact, a simple elastic-perfectly-plastic model (VonMises model) was employed. The spall failure criterion was used for solid materials. The material instantaneously fails when the hydro static pressure reaches its negative limiting value.
3. Hypervelocity impact simulations 3.1. Axisymmetric simulations: normal collisions Two cases of the normal impact are numerically simulated. In the first case (Case 1), a 21 km diameter asteroid (specific gravity: 4) impacts on the Earth with the velocity of 20 km/s, while in the second case (Case 2), a 13 km diameter comet (specific gravity: 1) collides with the Earth with the velocity of 60 km/s. Fig. 1 shows the simulation domain. Since calculations were carried out with axisymmetic geometry, only half of the region is displayed. The space resolution is 1 km/grid point and 300 and 361 cells were used in the radial and axial directions. The simulation includes the 5 km deep ocean and the 6 km thick crust in between the atmosphere and the mantle. Water is used as the oceanic material. The oceanic sediment layer is neglected since the thickness is approximately 1 km and no significant influence is expected with the current grid resolution. Major elements of the crust are SiO2 (60%) and Al2O3 (15%) and specific gravity of 2.8–2.9 while for the mantle they are SiO2 (45%) and MgO (35%) with specific gravity of 3.3. It is assumed that cordierite (2MgO–2Al2O3–5SiO2) and steatite (MgO–SiO2) have similar properties with the crust and mantle,
ARTICLE IN PRESS 716
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
Fig. 1. Schematic of the computational domain for axisymmetric calculations.
respectively. Their mechanical properties used in the current study are taken from the matweb database (http://www.matweb.com). In this study, judging from the impact velocity, it is implicitly assumed that the large slower impactor (Case 1) is an asteroid and the smaller one with much faster velocity (Case 2 and Case A–D in the next section) is a comet. The material of the asteroid and comet are not known and they are assumed rather arbitrarily: the asteroid is pure limestone and the comet consists of 75% limestone and 25% ice. Their densities are 4000 and 3250 kg/m3, respectively. Numerical models and the physical parameters adopted in this study are summarized in Table 1. 3.2. Axisymmetric numerical results and discussions Fig. 2 shows the material distribution with velocity arrows at different times for Case 1. The times after impact are indicated in the figure. The shock wave propagating in the mantle is clearly seen in the figure at 10 s after the impact. The depth of the crater asymptotically reaches its maximum of about 70 km in about 30 s, while the diameter still increases 60 s after the impact. It is noted that the interface between the atmosphere and the mantle inside the crater becomes unstable at this late stage. Although the calculation stopped at this moment, it is expected that this is the beginning of the collapsing stage. From the figure, it is observed that the diameter of the crater is approximately 120 km. Fig. 3 shows the accumulation of different materials that are splashed into the atmosphere as a function of time. Between 10 and 15 km altitude, the material originated from the crust amounts the most and those from the seawater and the mantle follow. Between 20 and 30 km altitude, the material from the crust also amounts to the largest value and material from the seawater and the mantle is relatively less compared with the lower altitude. In both cases, practically no material from the asteroid is observed. Fig. 4 shows the material distribution with velocity arrows at different times for Case 2. Although the calculation stopped at 20 s after the impact, it looks as if the depth of the crater is almost at its maximum of about 60 km, which is less than Case 1. The final diameter is about 80 km and it is also smaller compared with that in Case 1. Considering the fact that the impact energy of Case 2 is about 1.7 times larger than that of Case 1, the result demonstrates that the size or the volume of a crater due to normal impact of an asteroid or a comet is a function of the impactor size as well as the impact kinetic energy. It is noted that the ejecta seems more violent compared with that in Case 1, and more melted impactor material recoils into the atmosphere.
ARTICLE IN PRESS T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
717
Table 1 Numerical models adopted and physical properties
Atmosphere
Material
Eq. of state
Constitutive Eq.
Failure model
Void
—
—
—
a
Ocean
Water
Polynomial Ref. density r0 ¼ 1000 Kg/m3
—
Hydro Tensile limit 3.0 MPa
Crust
Cordieriteb (Dolomite)
Tillotson Ref. density r0 ¼ 2800 Kg/m3
VonMises Shear Modulus: 53 GPa Yield Stress: 344 GPa
Hydro Tensile limit 27.5 MPa
Mantle
Steatiteb (Dolomite)
Tillotson Ref. density r0 ¼ 3300 Kg/m3
VonMises Shear Modulus: 45 GPa Yield Stress: 619 GPa
Hydro Tensile limit 68.8 MPa
Asteroid
Limestoneb (Limestone)
Tillotson Ref. density r0 ¼ 4000 Kg/m3
VonMises Shear Modulus: 53 GPa Yield Stress: 344 GPa
Hydro Tensile limit 27.5 MPa
Comet
Limestone+iceb (Limestone)
Tillotson Ref. density r0 ¼ 3250 Kg/m3
VonMises Shear Modulus: 53 GPa Yield Stress: 344 GPa
Hydro Tensile limit 27.5 MPa
a
Bakken and Anderson 1967 [19]. Tillotson EOS parameters for ‘dolomite’ and ‘limestone’ [20] except densities are used. The densities of asterod and comet are determined from astronomical considerations. b
Fig. 2. Axisymmetric numerical results of Case1: diameter; 21 km, impact velocity; 20 km/s.
ARTICLE IN PRESS 718
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
Fig. 3. Accumulation of ejected mass from different sources at different altitude as function of time for Case1: left; between altitudes of 10 and 15 km, right; between altitudes of 20 and 30 km.
Fig. 4. Axisymmetric numerical results of Case 2: diameter; 13 km, impact velocity; 60 km/s.
ARTICLE IN PRESS T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
719
Fig. 5. Accumulation of ejected mass from different sources at different altitude as function of time for Case 2: left; between altitudes of 10 and 15 km, right; altitudes of 20 and 30 km.
Fig. 5 shows the accumulation of different materials that are splashed into the atmosphere at the specified altitude. Although it is still very little, the material from the impactor is visible in Case 2 reflecting the more violent features of impact in the early stage. According to the geological observations, more material from the impactor (asteroid or comet) is expected. The main reason for the current result, in that a very little amount of material from the asteroid/comet is ejected into the atmosphere, is most likely due the assumption of normal impact. In general, the most probable impact angle is 451 and the situation can be quite different from that of normal impact. The effect of impact angle is investigated in the next section using two-dimensional (2D) simulations. Among the material parameters, the density is considered to be the most important one for hypervelocity impact. In the present work, actual densities of material are used. However, other Tillotson EOS parameters are taken from those either for Dolomite or Limestone as indicated in Table 1. Although we cannot quantitatively estimate how much error would be induced with this approximation, we can still hope that the numerical results give us useful information. Investigations of the Tillotson parameters that can control the shock-induced vaporization are needed in future in order to accurately estimate the amount of mass splashed into the atmosphere. 3.3. 2D simulations of oblique collisions Hypervelocity collisions of asteroid or comet with the Earth at different impact angles were carried out in order to investigate the effect of impact angle to the appearance of the ejecta and the melted impactor. Simulations were carried out with the 2D plane strain model due to restricted computational resources. The diameter of the impactor (comet) is 13 km and the impact velocity is set to 60 km/s. Note that the shape of the comet now is essentially cylindrical while it is spherical in the previous section. The impact angles measured from the vertical direction are 01, 301, 451 and 601 and are named as Case A, Case B, Case C and Case D, respectively. Fig. 6 shows the domain of numerical simulation. Other numerical parameters such as the grid resolution are the same as in the previous axisymmetric calculations. 3.4. 2D results and discussions Fig. 7 compares the results of axisymmetric calculation (Case 2) and 2D calculation for a 13 km diameter comet impacting normally to the earth with the impact velocity of 60 km (Case A). The figures displayed are at
ARTICLE IN PRESS 720
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
Fig. 6. Schematic of computational domain and asteroid impact angles: Atmosphere 50 km, Ocean 5 km, Crust 6 km, Mantle 300 km.
Fig. 7. Normal impact simulations with axisymmetric (left: Case 2) and 2D geometries (right: Case A) 4 s after impact.
4 s after the impact event. The kinetic energy is 6.3 1024 J for the axisymmetric case and 7.8 1020 J for the 2D case assuming the width of comet to be a unit length. Then, the kinetic impact energies per unit colliding area are obtained by dividing the kinetic energies with the projection areas of the comet as 5.1 1016 and 6.0 1016 J/m2 for the axisymmetric and 2D cases, respectively. The impact energy density is about 18% higher in the 2D case. Although the general appearances of the two are more or less the same, the 2D case ejects much more material into the atmosphere. The difference between the two is more than expected and it can be attributed to the different ways of how impact energy is diffused after the impact. It is also noted that the impact material (comet) is completely disintegrated in the axisymmetric case while it is only partially so in the 2D simulation. It is natural to consider the axisymmetric case to be more realistic and we have to keep this observation in mind.
ARTICLE IN PRESS T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
721
Fig. 8. Oblique impact simulations with 2D geometry: diameter; 13 km, impact velocity; 60 km/s, 4 s after impact, impact angles: Case A: 01, Case B: 301, Case C: 451, Case D: 601.
Fig. 9. Oblique impact simulations with 2D geometry: 13 km, impact velocity; 60 km/s, 8 s after impact, impact angles: CaseA: 01, Case B: 301, Case C: 451, Case D: 601.
Namely, although the simulations for Cases A–D may provide us with some general data of the effect of impact angle on the current impact problem, we still need three-dimensional simulations to get some quantitative conclusions. Fig. 8 compares 2D results (Cases A–D) at 4 s after the impact. As expected, the material is ejected almost into the direction of mirror reflection to the incident impact angle. It is also noticed that, for a larger impact
ARTICLE IN PRESS 722
T. Saito et al. / International Journal of Impact Engineering 33 (2006) 713–722
angle, the depth of crater is smaller and its diameter is larger. Fig. 9 shows the results at 8 s after the impact. Comparing it with Fig. 8, it is noticed that the angle dependency of ejecta flow gets weaker in later stages. The difference between Cases C and D at 4 s after the impact (Fig. 8) is quite notable but Cases C and D in Fig. 9 are not much different with each other. It is known that the shape of crater is circular regardless of the impact angle. The tendency observed here in Figs. 8 and 9 may be reflection of this known fact. 4. Conclusions Two cases of normal impact of an asteroid/comet with the earth are numerically simulated. It is demonstrated that the materials from different origin are splashed into the atmosphere. The size, i.e. the depth and the diameter of the crater, is a function of the diameter as well as the kinetic energy of the impactor. It is demonstrated that the crater size of a hypervelocity impact with higher impact energy can be smaller if the diameter of the impactor is smaller. The effect of impact angle to the flow of ejected materials is investigated using 2D numerical analysis. Although it is only qualitative, a tendency that the dependency of the appearance of mass ejection into the atmosphere on the impact angle is more significant in the early stage and becomes less notable at later time is observed. References [1] Kaiho K, Lamolda MA. Catastrophic extinction of planktonic foraminifera at the Cretaceous–Tertiary boundary evidenced by stable isotopes and foraminiferal abundance at Caravaca, Spain. Geology 1999;27:355–8. [2] Rampino MR, Adler AC. Evidence for abrupt latest Permian mass extinction of foraminifera: results of tests for the Signor–Lipps effect. Geology 1998;26:415–8. [3] Bowring SA, Erwin DH, Jin YG, Martin MW, David EK, Wang W. U/Pb zircon geochronology and tempo of the end-Permian mass extinction. Science 1998;280:1039–45. [4] Kaiho K, Kajiwara Y, Tazaki K, Ueshima M, Takeda N, Kawahata H, et al. Oceanic primary productivity and dissolved oxygen levels at the Cretaceous/Tertiary boundary: their decrease, subsequent warming, and recovery. Paleoceanography 1999;14:511–24. [5] Alvarez W. Toward a theory of impact crises. EOS (Transactions, American Geophysical Union), 1986; 67: 649, 653–655, 658. [6] Sheehan PM, Fastovsky DE, Barreto C, Hoffmann RG. Dinosaur abundance was not declining in a ‘‘3 m gap’’ at the top of the Hell Creek Formation, Montana and North Dakota. Geology 2000;28:523–6. [7] Erwin DH. The great Paleozoic crisis. New York: Columbia University Press; 1993 327p. [8] Retallack GJ. Permian–Triassic extinction on land. Science 1995;267:77–80. [9] Eshet Y, Rampino MR, Visscher H. Fungal event and palynological record of ecological crisis and recovery across the Permian–Triassic boundary. Geology 1995;23:967–70. [10] Jin YG, Wang Y, Wang W, Shang QH, Cao CQ, Erwin DH. Pattern of marine mass extinction near the Permian–Triassic boundary in south China. Science 2000;289:432–6. [11] Becker L, et al. Impact event at the Permian–Triassic boundary: evidence from extraterrestrial noble gases in fullerenes. Science 2001;291:1530–3. [12] Kaiho K, Kajiwara Y, Miura Y. End-Permian catastrophe by a bolide impact: evidence of a gigantic release of sulfur from the mantle: reply. Geology 2002;September:856. [13] Pierazzo E, Melosh HJ. Hydrocode modelling of Chicxulub as an oblique impact event. Earth Planet Sci Lett 1999;165:163–76. [14] Jones AP, Price GD, Price NJ, DeCarli PS, Clegg RA. Impact induced melting and the development of large igneous provinces. Earth Planet Sci Lett 2002;202:551–61. [15] Stoeffler D, Artemieva NA, Pierazzo E. Modeling the Ries–Steinheim impact event and the formation of the Moldavite strewn field. Meteoritics Planet Sci 2002;37:1893–908. [16] Wilkins ML. Calculation of elastic-plastic flow. Lawrence Livermore Laboratory Report, UCRL-7322, revision 1, 1973. [17] Robertson NJ, Hayhurst CJ, Fairlie GE. Numerical simulation of explosion phenomena. Int J Comput Appl Technol 1994;7(3–6):316–29. [18] Tillotson JH. Metallic equations of state for hypervelocity impact. GA-3216, General Atomic, CA, July, 1962. [19] Bakken LH, Anderson PD. The complete equation of state handbook. SANDIA SCLTM-67-118, November, 1967. [20] Hageman LJ, Walsh JM. HELP: a multi-material Eulerian program for compressible fluid and elastic-plastic flows in two space dimensions and time. BRL Contract Report, vol. 1(39), 1971. p. 26.