Numerical modeling of the 2013 meteorite entry in Lake Chebarkul, Russia

Abstract

The results of the numerical simulation of possible hydrodynamic perturbations in Lake Chebarkul (Russia) as a consequence of the meteorite fall of 2013 (15 February) are presented. The numerical modeling is based on the Navier–Stokes equations for a two-phase fluid. The results of the simulation of a meteorite entering the water at an angle of 20${}^{\circ }$ are given. Numerical experiments are carried out both when the lake is covered with ice and when it is not. The estimation of size of the destructed ice cover is made. It is shown that the size of the observed ice hole at the place of the meteorite fall is in good agreement with the theoretical predictions, as well as with other estimates. The heights of tsunami waves generated by a small meteorite entering the lake are small enough (a few centimeters) according to the estimations. However, the danger of a tsunami of meteorite or asteroid origin should not be underestimated.

Introduction

On 15 February 2013 at 09:20 LT (local time), in the vicinity of the city of Chelyabinsk, Russia (Fig. 1), a meteorite exploded and collapsed in the earth's atmosphere as a result of inhibition. Its fall was accompanied by a series of atmospheric explosions and propagation of shock waves over the territory of the Chelyabinsk region. Small fragments of the meteorite came down on the region. It is the largest of the known celestial bodies which have fallen to the ground after the Tunguska meteorite in 1908. The following characteristics of the Chelyabinsk meteorite are given in recent publications (Zetser, 2014; Ionov, 2013; Kopeikin et al., 2013; Emel'yanenko et al., 2013; Popova et al., 2013; Berngardt et al., 2013; Gokhberg et al., 2013; Krasnov et al., 2014; Seleznev et al., 2013; De Groot-Hedlin and Hedlin, 2014):

• the meteorite with a diameter of 16–19 m flew into the earth's atmosphere at about 20${}^{\circ }$ to the horizon at a velocity of $\sim$ 17–22 km s${}^{-\mathrm{1}}$.

• The meteorite was destroyed in several stages, and the main explosion occurred at an altitude of about 23 km. The analysis of the meteorite fall led to the estimation of the explosion energy from 380 to 1000 kt of TNT.

• Surviving fragments of the Chelyabinsk meteorite supposedly flew at velocities up to 150–300 m s${}^{-\mathrm{1}}$, and according to the updated data the velocity was 156 m s${}^{-\mathrm{1}}$.

• The largest meteorite fragment, weighing about 550 kg, was recovered from Lake Chebarkul near the Krutik Peninsula, and it had an irregular oval shape with an average outer diameter of about 1 m.

• The celestial body exploded when it hit the ice and water. Its fragments flew more than 100 m and formed an ice hole of a round shape with a diameter of about 8 m (Fig. 1).

The ground-penetrating radar (GPR) survey of the crash site showed that the crater, which was formed by the meteorite impact on the bottom, was observed 30 m away from the ice hole (Kopeikin et al., 2013).

The ice hole that was formed on the surface of Lake Chebarkul.

[Figure omitted. See PDF]

The results of the numerical modeling of the processes that began in the water after the meteorite entered Lake Chebarkul are presented in this paper. We consider two cases – the occurrence of a meteorite in the ice-covered lake (as it was in February 2013) and in the lake without ice to estimate the amplitude of the wave in the case of a possible fall in the summer. First of all, in Sect. 2 a well-known parametric model of the tsunami source of meteoric origin (Ward and Asphaug, 2000) which estimates the disturbances in the water at the site of the meteorite entry is used. It is shown that the size of the source is twice the size of the observed ice hole; that is why this model cannot explain the observed phenomenon. In Sect. 3 we describe the hydrodynamic model of wave generation based on the Navier–Stokes equations as well as the parameters of its numerical discretization. In Sect. 4 we present the results of numerical experiments on the tsunami wave generation in open water. The results of numerical experiments on disturbance generation in the lake covered with ice are presented in Sect. 5. The zones of maximum and minimum pressure which could potentially lead to the destruction of the ice cover are detected. In Sect. 6 the results of estimations of the ice cover destruction, based on the calculation of the stresses, are shown. The area of ice destruction is estimated according to the famous semi-empirical formulas and pressure values obtained in the numerical experiments. This estimation is in good agreement with the observed data, as well as with the estimates made previously in Ivanov (2014). The results are summarized in Sect. 7.

Preliminary estimates of wave heights in the lake with a free surface

To estimate the possible parameters of the water surface displacement when a meteorite falls, we use a simplified model in which the generation of waves by a meteorite entering the water is parameterized by certain initial conditions (Ward and Asphaug, 2000; Mirchina and Pelinovsky, 1988; Kharif and Pelinovsky, 2005; Levin and Nosov, 2009; Torsvick et al, 2010). It is suggested that in the initial stage of a crater formation a meteorite, vertically entering the water, creates a radially symmetrical cavity on the water surface, which can be described by a simple function: $$\begin{array}{}& {\displaystyle }& {\displaystyle }{\mathit{\eta }}^{\mathrm{imp}}\left(r\right)=D_{\mathrm{C}}\left(\mathrm{1},-,r^{\mathrm{2}},/,R_{\mathrm{C}}^{\mathrm{2}}\right),r\le R_{\mathrm{D}},& {\displaystyle }& {\displaystyle }{\mathit{\eta }}^{\mathrm{imp}}\left(r\right)=\mathrm{0},r>R_{\mathrm{D}},\end{array}$$ where $D_{\mathrm{C}}$ is the depth of the cavity, and $R_{\mathrm{C}}$ and $R_{\mathrm{D}}$ are inner and outer radii of the cavity, respectively. In the case of $R_{\mathrm{D}}$ $=$ $\sqrt{\mathrm{2}}{R}_{\mathrm{C}}$ the water ejected from the cavity forms the outer splash–ring structure typical of the fall of the object in the water and the volume of which corresponds exactly to the volume of water discharged from the cavity.

Considering that the meteorite kinetic energy is converted into the potential energy of the water level displacement, the following simple analytical formulas are derived to calculate the radius and depth of the cavity (see, for instance, Ward and Asphaug, 2000): $$\begin{array}{cc}{\displaystyle }{D}_{\mathrm{C}}=& {\displaystyle }\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\epsilon }{\mathit{\rho }}_{\mathrm{I}}{R}_{\mathrm{I}}^{\mathrm{3}}{V}_{\mathrm{I}}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{w}}g{R}_{\mathrm{C}}^{\mathrm{2}}}}},R_{\mathrm{C}}=R_{\mathrm{I}}{\left(\mathrm{2},\mathit{\epsilon },{\displaystyle \frac{V_{\mathrm{I}}^{\mathrm{2}}}{g{R}_{\mathrm{I}}}}\right)}^{\mathit{\delta }}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{I}}}{{\mathit{\rho }}_{\mathrm{w}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\\ & {\displaystyle }& {\displaystyle }\left({\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{I}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}-\mathit{\delta }},{\left({\displaystyle \frac{\mathrm{1}}{q{R}_{\mathrm{I}}^{\mathit{\alpha }-\mathrm{1}}}}\right)}^{\mathrm{2}\mathit{\delta }}\right),\end{array}$$ where ${\mathit{\rho }}_{\mathrm{w}}$ is density of water; $g$ is gravitational acceleration; $\mathit{\epsilon }$ is the proportion of the kinetic energy of the meteorite, converted into the tsunami energy; ${\mathit{\rho }}_{\mathrm{I}}$, $R_{\mathrm{I}}$ and $V_{\mathrm{I}}$ are density, radius and velocity of the meteorite $q$; and $\mathit{\alpha }$ are the coefficients associated with the properties of the meteorite and the water layer. According to Levin and Nosov (2009) about 16% of the kinetic energy of the falling body is converted into the energy of tsunami waves. The parameter $\mathit{\alpha }$ is 1.27, and the value of $q$ and $\mathit{\delta }$ is calculated as follows: $${\displaystyle }{\displaystyle }\mathit{\delta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\alpha }+\mathrm{2}}},q\approx \mathrm{0.39}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{I}}}}\right)}^{\mathrm{0.26}}{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{I}}^{\mathrm{0.27}}}}.$$ If the diameter of a meteorite is 1 m, the density is 3.3 g cm${}^{-\mathrm{3}}$ and the water entering velocity is 156 m s${}^{-\mathrm{1}}$, the inner radius of the cavity is equal to $R_{\mathrm{C}}$ $=$ 5.6 m, the outer radius is $R_{\mathrm{D}}$ $=$ 7.8 m and the depth of the cavity formed on the surface of Lake Chebarkul is $D_{\mathrm{C}}$ $=$ 3.2 m approximately. The resulting cavity is schematically depicted in Fig. 2.

The obtained estimation of the cavity diameter (15.5 m) exceeds twice the size of the ice hole (7–8 m) formed by the fall of a Chelyabinsk meteorite fragment (Fig. 1). The estimation of the cavity size is made for the conditions of open water. If we take into account the ice cover, which prevents the water from forming vertical upward discharges at the edges of the resulting cavities, and consider only part of the underwater cavity (indicated in Fig. 2 by dashed lines), we will obtain a diameter of a little more than 10 m, which is closer to the observed size of the ice hole. Another reason for the differences in the estimation of the source size, in our opinion, might be the wrong extrapolation of empirical formulas for the small size meteorites, especially since it does not enter the water vertically.

The initial disturbance of the water surface at the source.

[Figure omitted. See PDF]

The hydrodynamic model of the tsunami source, based on the Navier–Stokes equations

The models based on the numerical solution of the Navier–Stokes equations (Ferziger and Peric, 2002), which have become popular in solving the problems of asteroid tsunami, allow simulating a real meteorite entering the water (Kozelkov et al., 2015, 2016a). In this paper we solve the Navier–Stokes equations for a two-component (water and air) incompressible fluid in a gravitational field, and both air and water are considered to be incompressible. The incompressibility/compressibility of air will mainly influence the rate of change of the wave parameters in the source during the collapse of the cavity. Qualitatively, the process will be described correctly, but the quantitative difference can be significant. This is directly related to the Mach number, determined by the characteristic velocity of the meteorite, which in our calculations is 156 m s${}^{-\mathrm{1}}$, which is approximately a Mach number of 0.5. It is known that the compressibility of the medium must be taken into account when the Mach number exceeds 0.3. In our case, the Mach number is slightly higher than the recommended value to take into account the compressibility, and therefore, in the first stage of calculations, the compressibility can be neglected. Of course, this will introduce a certain error in the collapse of the cavern, but it is estimated that it will not exceed 15–20 % for speed and pressure. After the wave has emerged from the source, the effect of compressibility on the process of its propagation is negligible.

The Navier–Stokes equations are used to describe each component of the fluid. Their integration into a single system is done using the “volume-of-fluid” (VOF) method (Hirt and Nichols, 1981). The numerical solution of the resulting system is computed using the methods SIMPLE/PISO (Ferziger and Peric, 2002; Issa, 1986). The movement of the meteorite is modeled using the immersed boundary method (IBM; Mittal and Iaccarino, 2005). This method involves the allocation of cells completely or partially occupied by a rigid body in the computational field and resistance force, proposed in Mohd-Yusof (1997). This approach to modeling the moving rigid body is quite simple. It does not require dynamic rebuilding of the computational grid, thus giving good results for the practical tasks that do not require a detailed description of the boundary layer near the surface of a solid body (Posa et al., 2011).

The numerical solution of coupled water wave–ice deformation differential equations is computationally intensive. Due to the absence of observations related to waves induced by this event, there is no possibility to refine the calculations, and for the first step the simplified problem can be considered with “rigid-wall” boundary conditions at all spatial boundaries of the computational domain. In the present calculations zero initial conditions for the velocity field in the water, the water level and the standard hydrostatic distribution of the pressure are used. It assumes a zero velocity and zero value of the pressure gradient and the volume fraction of the components, as well as zero shear stress. In fact, the calculation stops when the perturbation reaches the lateral boundaries so that these boundaries do not actually affect the results. The upper part of the computational domain is an open border with a given static pressure of 1 atm and a zero velocity gradient. The bottom of the basin is regarded as being non-deformable, and when a rigid body reaches the bottom, its velocity is artificially set to 0. The value of the component volume fraction at the upper boundary is different and depends on the direction of the flow: the inlet flow volume fraction of water is equal to 0, and that of the air is equal to 1. At the output flow the gradient of the volume fraction of each component is assumed to be 0. Physically, this condition means that all the components are free to leave the settlement area through the upper boundary, while only the air enters it. For future studies, of course, the “ice” block of hydroelastic equations with adequate breaking conditions should be added into the model.

The described model is implemented in the software package LOGOS, which is used for the numerical simulations presented below. The LOGOS software package has been tested for the given class of problems and has shown sufficiently reliable results (Kozelkov et al., 2013, 2014, 2016b, Kozelkov and Kurulin, 2015; Volkov et al., 2013; Betelin et al., 2014; Deryugin et al., 2015).

The bathymetric data of Lake Chebarkul are not available, and the depth of a model region relative to the surface is taken as a constant of 10 m, corresponding to the average depth of Lake Chebarkul. The height above the zero level surface of the water is 40 m. So the discrete area of the model was constructed with a size of 160 $\times$ 160 $\times$ 50 m. It is a non-structured three-dimensional (3-D) grid of truncated polyhedrons of arbitrary shape (Fig. 3). This type of grid is the only possible one for the areas of complex geometric configurations. It is built by the preprocessor of the software package LOGOS.

Model discrete area (a) and a fragment of the computational grid at the site of the meteorite entering (b). Zones of mesh refinement near the free surface and the area of the movement of the meteorite are shown in black.

[Figure omitted. See PDF]

Perturbation of water at different points in time when the meteorite enters the water at an angle of 20${}^{\circ }$ (the meteorite comes from left to right).

[Figure omitted. See PDF]

The perturbation of the volume fraction (blue – water; red – air) in the cross section of the computational domain taken along the centerline of the meteorite fall (the body enters from right to left at an angle of 20${}^{\circ }$) at different times.

[Figure omitted. See PDF]

The form of the water surface at different times.

[Figure omitted. See PDF]

Numerical experiments on the generation of tsunami waves in open water

In the first numerical experiment within the hydrodynamic model described above, we examine the case of wave generation in the water without ice cover. Figure 4 illustrates a three-dimensional distribution of the water when the meteorite enters the water at an angle of 20${}^{\circ }$ to the vertical. The diameter of the meteorite is 1 m, its density is 3.3 g cm${}^{-\mathrm{3}}$ and the velocity of the entry is 156 m s${}^{-\mathrm{1}}$. The visualization is carried out by the post-processing software package LOGOS.

A powerful upward surge of water is formed at the first moments of the collision (the formation of the splash), and then waves are generated. The head wave has a positive polarity (the hump). Under the water a strong air bubble is generated by the movement of the body in the water. At first this bubble has a cylindrical structure, and then, during the destruction and surfacing, it becomes oval. The destruction of the underwater bubble near the surface generates an additional, “secondary” wave. This process is in good agreement with the classic description of a body falling into water (Aristoff, 2009).

The anisotropy of wave propagation is visible – the waves propagate faster in the direction of a meteorite fall than in the opposite direction, and their amplitudes decrease rapidly. The details of the cross section of perturbations are shown in Figs. 5 and 6. Figure 5 shows the surface of the water in the cross section taken along the centerline of the meteorite fall on a large scale, and Fig. 6 shows it in the whole computational domain.

The meteorite passes through about a quarter of the water column in 0.05 s and almost reaches the bottom in 0.15 s. At 0.25 s the meteorite is already lying on the bottom, so it has passed a distance of about 27 m from the ice hole on the surface. This result is in good agreement with the detailed analysis of a geo-radar cross-sectional survey (Kopeikin et al., 2013). When the meteorite falls, it forms a cavity with a diameter of about 4 m near the surface in 0.05 s. The diameter becomes 7–8 m before the start of the collapse, which begins when the meteorite reaches the bottom. This value corresponds to the diameter of the ice hole formed as a result of a real meteorite fall. The collapse of the cavity begins 0.25 s later and ends within 0.5 s. At the time of the cavity collapse a “primary wave” with a height of about 3.5 m and a large air pocket in the water column are generated. The height of the “core” wave corresponds to the theoretical estimation (Eq. 2) in the framework of the parametric source. The height of “splash” lifting at the time of the cavity collapse is close to the mark of 10 m, but these splashes do not participate in the surface wave generation (see time 1–2 s).

The transformation of the “air” pocket starts at 0.5 s. At the time of 2 s it reaches the water surface and generates a secondary wave about 1 m high. Two seconds after a meteorite enters the water, the height of the primary wave is about 1 m; it completely collapses and begins to move away from the source. At this time, the secondary wave of almost the same height as the primary one is generated. The process of the collapse ends 3 s after the meteorite fall. The maximum height of the splash, formed after the fall, is observed approximately 0.75 s later. It is about 10 m, which corresponds to the depth of the basin (Fig. 5). This height is visible at the withering “cap” of the splash, which breaks away from the main water discharge of about 4 m.

At the time of 1 s the first wave, which reaches a height of about 1 m at 1.5 s, begins to propagate from the source. At 2 s this wave collapses under the air pocket impact, which reaches the water surface and generates the so-called secondary wave. The primary and secondary waves collide with each other in our numerical experiment. The maximum height of the generated wave is observed at the time of 3 s, and it is about 3 m high. During the next second the water near the site of the meteorite fall oscillates, and 5 s after the meteorite enters the water two waves go out from the source. They can be clearly seen in Fig. 6 at the time of 6 s – their amplitude is only about 20 cm. These two leading waves go out from the source and then propagate on Lake Chebarkul. Generally speaking, there should be two waves formed by the meteorite – a primary wave from the intermediate cavity and a secondary wave generated by the collapse of the cavity and by the underwater bubble (Aristoff, 2009; Kozelkov et al., 2015). However, here the waves had an impact on each other and created a single wave, coming from the source. Subsequently, this wave was attenuated due to the cylindrical divergence and dispersion by law $r^{-(\mathrm{0.5}-\mathrm{1})}$; therefore, its height at a distance greater than 100 m will not exceed a few centimeters, which is why such a wave is not dangerous outside the area of the meteorite entering the water.

In conclusion we would like to point out that the difference between the sizes of the ice hole, obtained by using the parametric model and within the Navier–Stokes equations, may be related to the fact that the model does not take into account the angle of the fall. According to Fig. 6 (at 0.2 s – the meteorite immersion is total) the major discharge of water falls on the side of the cavity, which is located along the moving body. Its diameter is just 7–8 m. If the fall had been vertical, the cavity would have looked different. Figure 7 illustrates the cavities resulting from the different falls when the meteorite hits the bottom.

One can see that the geometry of the cavities varies and two intense waves are generated if the fall is vertical. In the fall at an angle a wave is formed along the meteorite's movement, while behind it the disturbances are minimal, and they will go in this direction just after the collapse of the “left” side.

The cavity formed by the meteorite fall at an angle (a) and in a vertical fall (b).

[Figure omitted. See PDF]

The variants of the grid model for the arbitrary shoreline with the marked meteorite fall zone (a – the full shoreline of the lake; b – the marked fragment of the meteorite fall site).

[Figure omitted. See PDF]

The tsunami propagation along Lake Chebarkul.

[Figure omitted. See PDF]

As was already noted, the bathymetric map of Lake Chebarkul is not available. It is possible to examine tsunami wave propagation on the lake using digitized coastline and constant average depth (for Lake Chebarkul it is 10 m). The map of Lake Chebarkul's shoreline and its digitized version, built with a model grid (right), are shown in Fig. 8. As we can see, the coastline has a complicated configuration, and it is advisable to construct a model grid with the help of an automatic generator of arbitrary unstructured grids.

The grid constructed by this generator is shown in Fig. 8 (left panel). In the area free from the shoreline (open water) the grid has mainly a hexagonal structure. Lake dimensions in length and width are approximately the same and are about 2 km. With the meteorite diameter of 1 m it is reasonable to allocate a fall area to grind the model grid to the desired level (Fig. 8, right panel) – at least 10 cells along the meteorite diameter (typical cell size of 10 cm). It is obvious that in the case of building a grid for the entire lake with the characteristic size, its size will be very large, which will significantly affect the account time and the CPU box used.

A typical cell size outside the fall area can be assumed optional, here it ranges from 10 to 50 m in the horizontal plane. To simulate the wave propagation accurately, it is necessary to thicken the grid near the interface of air and the water surface (Fig. 8, top right panel). Usually, this thickening is carried out by the given law of geometric progression, so the typical cell sizes near the surface can be selected for the desired size in order to track the desired wave height. The grid consists of approximately 15 million cells. The results of tsunami propagation simulation on Lake Chebarkul without ice cover are shown in Fig. 9.

Already at the ninth second we can see that the wave has sufficiently small amplitude of about 10 cm. About 30 s after the fall the wave reaches the nearest shore, and its height is also about 10 cm. At the 45th second we can see the bounce of the wave from the shore. This, in its turn, suggests a fairly good grid resolution which allows for reproduction of centimeter amplitude waves. At the 60th second we can see how the bounced wave follows the primary one, which arose after the collapse of the water crater formed in the course of a meteorite fall. Ninety seconds after the fall, the primary wave and the “catch-up way” are practically damped. The wave does not reach the opposite bank.

Numerical experiments taking into account ice

The Chelyabinsk meteorite fell in wintertime when the ice thickness on Lake Chebarkul was 70–80 cm. Modeling the meteorite fall into the water with punching of the ice is quite a challenge. We have to take into account the relationship of processes of hydrodynamics and destruction. Therefore, here to estimate the destruction of the ice cover and the possible splash zone of the water on the ice surface, we simulate an idealized case of a meteorite fall. Let a meteorite enter an already “punched” hole in the ice surface. In this case the ice is simulated as the rigid-wall boundary condition, and the zero water level is placed in the hole. The hole in the computational experiment is rectangular with a size of 160 $\times$ 130 cm${}^{\mathrm{2}}$, close to the size of the meteorite. The meteorite entering parameters are exactly the same as in the previous case. Figure 10 shows the domain of the computational experiment.

The perturbation of the volume fraction (blue – water; white – air) in the cross section of the computational domain taken along the centerline of the meteorite fall at different times in the presence of ice (the body enters from right to left at an angle of 20${}^{\circ }$).

[Figure omitted. See PDF]

The pressure field in the computational domain (the horizontal line – ice; the short black strip – entrance area).

[Figure omitted. See PDF]

The height of the splash formed by the fall of the meteorite into the hole is over 10 m at the time of 0.2 s. Thus the splash rises almost vertically upwards, without any deviation in any direction. The main part of the splash consists of splashes, and the height of the discharge of the water bulk is about 2.5 m.

About 1 s after the fall the splash enters the final stage of the collapse on the ice surface. When the splash collapses, it begins to deviate in the direction towards the entering body (left), without generating strong disturbances on the ice surface and under it. The maximum height of the wave falling on the ice surface is 1 m at the time of 1 s. The collapsing splash spans a distance of about 10 m on the ice surface in the direction of the fall (Fig. 10 at time 2 s).

The resulting size of the ice hole during the Chelyabinsk event was several times larger than the diameter of the meteorite; i.e., the observed destruction of the ice cover around the area of the fall was much larger than the area of the collision. To analyze this phenomenon it is necessary to examine the pressure distribution in the studied area. Figure 11 shows the distribution of pressure in the entire computational domain at different times during the meteorite moving into the water column.

As shown in Fig. 11, the maximum pressure zones are observed in the frontal point of the body contact with water (in the so-called breaking point – point 1 in the figure) and under the ice surface at the start of the body immersion (point 2 in the figure).

A shock wave is observed in the fluid immediately after the meteorite enters the water. It generates a zone of high pressure around the ice surface (point 4).

Around the body the fluid accelerates from point 1 (breaking point) to point 2 (and a symmetrical point on the other side of the body). At the breaking point the pressure is at its maximum. As the distance from the middle of the body increases (point 2) the movement is slowed down and the pressure is reduced. In the central part of the back side of the body (near point 3) the pressure rises again, but in a very narrow area. This picture corresponds to the classic description of the boundary layer separation flow around the body with a blunt stern (Schlichting, 1960).

The picture of the pressure distribution shows the area of the possible ice destruction. It consists of two domains – a zone of high pressure in the shock wave (point 4) and a low-pressure zone, which is formed after the body has passed (point 5). In the zone of increased pressure, the possible destruction of the ice is observed when it is pushed upwards, and in the zone of reduced pressure the possible destruction occurs when it moves downwards.

Estimation of the ice cover destruction

An estimation of the ice cover destruction on Lake Chebarkul after the fall of the Chelyabinsk meteorite fragments was recently made in Ivanov (2014). The authors used both well-known empirical data of the ice destruction by explosions in the water and direct numerical simulation of the body impact on the ice. Both estimates, made for the meteorite vertical fall with a velocity of 100 m s${}^{-\mathrm{1}}$, predict that the size of the ice hole is 6–8 m by explosive or impact action with 3–10 MJ of energy, which corresponds to the energy of the Chelyabinsk meteorite.

Here we use a different approach based on the method of calculating the strength of a material (Barber, 2010). This method is built on the assumption that the determining parameter of the structure reliability is stress or, more precisely, the state of tension at a point. The estimated value of the stress is compared with the maximum permissible stress value for the material obtained from experimental studies. The conclusion about the strength of the material is made from the comparison of its calculated and limit stresses. The strength property has the form $${\displaystyle }{\displaystyle }{\mathit{\sigma }}_{\max }\le \left[\mathit{\sigma }\right],$$ where ${\mathit{\sigma }}_{\max }$ is the maximum calculated stress that occurs in the material and [$\mathit{\sigma }$] is the ultimate strength of the material.

Upon the impact of a meteorite, the fluid under the influence of excess pressure rushes straight up. We can consider the ice cover to be a thin elastic plate of infinite size, resting on an elastic foundation (Peschanskii, 1967) – in this case it is water. A bend occurs under this loading. In the presence of bending due to the lateral load, Eq. (4) is transformed as follows: $${\displaystyle }{\displaystyle }M/W\le \left[\mathit{\sigma }\right],$$ where $M$ is the bending moment due to lateral load and $W$ is the moment of resistance to bending. If the load, uniformly distributed over the area of a circle of some radius, affects infinite-size ice, the maximum bending moment can be calculated to determine the load. The expression for the stress produced by the action of the maximum bending moment for winter ice at temperatures of $-$25 ${}^{\circ }$C is contained in (Peschanskii, 1967) $${\displaystyle }{\displaystyle }q={\displaystyle \frac{{\mathit{\sigma }}_{\max }{h}^{\mathrm{5}/\mathrm{4}}}{\mathrm{31}\left(\mathrm{0.76},r_{\mathrm{0}},-,h^{\mathrm{34}}\right)}},$$ where $h$ is the thickness of the ice layer and $r_{\mathrm{0}}$ is the radius of the area where the load is distributed. It should be noted that this formula is obtained to determine the strength of the ice cover to allow its use as a crossing. These expressions allow us to estimate the area of lake ice cover destruction that resulted from the under-ice pressure which was induced by underwater motion of the body. Considering the spatial distribution of the shock pressure, one can assume that this formula gives an adequate estimation.

According to Fig. 11 the maximum under-ice pressure is 2 MPa, and the radius of the zone is 2 m. For an ice thickness of 1 m and an ambient temperature of about $-$25${}^{\circ }$, the stress calculated by Eq. (6) has a value of 35 $\times$ 103 kPa. The data of the values of river ice breaking bending stress are given in Peschanskii (1967). The maximum available value is shown for a temperature of $-$15 ${}^{\circ }$C and amounts to 50 kPa, which is 3 orders smaller than the obtained value. Therefore, it can be argued with high probability that at this point the ice cover has been destroyed. That is why the size of the ice hole is bigger than the radius of the meteorite.

The pressure around the ice is also about 2 MPa in the reduced-pressure zone. This fact means the destruction of ice in this zone. Thus, according to Fig. 9 the zone of destruction is equal to approximately 6 m from the point of the direct fall; i.e., the overall ice hole size obtained in the calculation is about 7 m, which agrees well with the observed data. This estimate is also in good agreement with Ivanov (2014), where it is noted that an ice hole with a diameter of about 6–7 m is to be formed in the lake if the velocity of the body fall is about 100 m s${}^{-\mathrm{1}}$.

A similar method of estimating the ice cover destruction is given in Fransson (2009). It suggests counting the capacity of ice using the results of Peschanskii (1967) and Barber (2010). Here, the force that causes the cracks in the ice is determined on the basis of the ice cover characteristic length $L$ in accordance with (Peschanskii, 1967) $$\begin{array}{}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{cr}}={\displaystyle \frac{\mathit{\pi }\mathit{\beta }}{\mathrm{3}(\mathrm{1}+n)f\left(\mathit{\beta }\right)}}{\mathit{\sigma }}_f{h}^{\mathrm{2}},& {\displaystyle }& {\displaystyle }f\left(\mathit{\beta }\right)=(\mathrm{0.6195}-\ln \mathit{\beta }){\displaystyle \frac{\mathit{\beta }}{\mathrm{2}}}+{\displaystyle \frac{\mathit{\pi }{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{64}}}+\mathrm{\dots },\end{array}$$ where $\mathit{\beta }$ $=$ $r_{\mathrm{0}}/L$ is the relative radius of the load; $f\left(\mathit{\beta }\right)$ is a function of the stress intensity; and ${\mathit{\sigma }}_f$ is the bending strength of the ice cover.

The strength leading to the ice cover destruction is determined by the semi-empirical dependence from $${\displaystyle }{\displaystyle }{P}_u=\mathrm{1.25}(\mathrm{1}+\mathrm{1.68}\mathit{\beta }){\mathit{\sigma }}_f{h}^{\mathrm{2}}.$$ Using the same data as before for the estimation, according to Eq. (9) we find that the value of the ice cover stress is 41.3 $\times$ 103 kPa, which almost corresponds to the previous estimation. It is 3 orders greater than the ice breaking bending stress (Peschanskii, 1967). Thus, both estimations allow explaining the observed size of the ice hole formed after the meteorite entering the water, as well as the absence of ice destruction over large areas by weak tsunami waves.

Conclusions

The results of the numerical modeling of perturbations formed on a meteorite entering Lake Chebarkul, Russia, are shown. The initial data correspond to the event that occurred on 15 February 2013. The characteristics of the waves both for the water with the ice cover and without it are found. In the second case it is shown that at a distance of about 100 m from the source the wave height is a few centimeters and ceases to be dangerous. In the case of a meteorite entering the lake covered with ice, the splash which has a height of about 10 m is formed and collapses quickly in the direction of the body motion. The main part of the splash consists of splashes, and the height of the water bulk discharge is about 2.5 m. The collapse zone of the splash is about 10 m. To estimate the ice cover destruction, approximate formulas are presented. They confirm the breaking of the ice at the site of a meteorite entering the lake. Numerical computations and estimations correctly predict the diameter of the ice hole observed on the lake after a meteorite entering it, and they are in good agreement with the estimates made in Ivanov (2014).

So, from the viewpoint of tsunami formation, the energy of the Chelyabinsk meteorite fragment that fell into the water is not enough to generate large waves. The whole effect is manifested only at the point of a meteorite entering the water. In the case of a larger meteorite fragment entering the water, tsunami waves can be significant, and estimations of the danger of these events should be made for inland lakes and seas.

All data sets of our computations are given in the manuscript.

The authors declare that they have no conflict of interest.

Acknowledgements

The presented results were obtained with the financial support of the grant of the President of the Russian Federation for state support of leading scientific schools of the Russian Federation (NSH-6637.2016.5) and RFBR grants (15-45-02061 and 16-01-00267). Edited by: I. Didenkulova Reviewed by: C. L. Mader and two anonymous referees